Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Contains types used for a Dimer Method calculations
10 : !> \par History
11 : !> -Luca Bellucci 11.2017 - CNR-NANO, Pisa
12 : !> add K-DIMER from doi:10.1063/1.4898664
13 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
14 : ! **************************************************************************************************
15 : MODULE dimer_methods
16 : USE bibliography, ONLY: Henkelman2014,&
17 : cite_reference
18 : USE cp_log_handling, ONLY: cp_get_default_logger,&
19 : cp_logger_type
20 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
21 : cp_print_key_unit_nr
22 : USE cp_subsys_types, ONLY: cp_subsys_get,&
23 : cp_subsys_type
24 : USE dimer_types, ONLY: dimer_env_type,&
25 : dimer_fixed_atom_control
26 : USE dimer_utils, ONLY: get_theta
27 : USE force_env_methods, ONLY: force_env_calc_energy_force
28 : USE force_env_types, ONLY: force_env_get
29 : USE geo_opt, ONLY: cp_rot_opt
30 : USE gopt_f_types, ONLY: gopt_f_type
31 : USE input_constants, ONLY: do_first_rotation_step,&
32 : do_second_rotation_step,&
33 : do_third_rotation_step
34 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
35 : section_vals_type
36 : USE kinds, ONLY: dp
37 : USE motion_utils, ONLY: rot_ana,&
38 : thrs_motion
39 : USE particle_list_types, ONLY: particle_list_type
40 : #include "../base/base_uses.f90"
41 :
42 : IMPLICIT NONE
43 : PRIVATE
44 :
45 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dimer_methods'
47 :
48 : PUBLIC :: cp_eval_at_ts
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief Computes the dimer energy/gradients (including the rotation of the dimer)
54 : !> \param gopt_env ...
55 : !> \param x ...
56 : !> \param f ...
57 : !> \param gradient ...
58 : !> \param calc_force ...
59 : !> \date 01.2008
60 : !> \par History
61 : !> none
62 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino]
63 : ! **************************************************************************************************
64 1816 : RECURSIVE SUBROUTINE cp_eval_at_ts(gopt_env, x, f, gradient, calc_force)
65 : TYPE(gopt_f_type), POINTER :: gopt_env
66 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
67 : REAL(KIND=dp), INTENT(OUT) :: f
68 : REAL(KIND=dp), DIMENSION(:), POINTER :: gradient
69 : LOGICAL, INTENT(IN) :: calc_force
70 :
71 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_eval_at_ts'
72 :
73 : INTEGER :: handle, iw
74 : LOGICAL :: eval_analytical
75 : REAL(KIND=dp) :: angle1, angle2, f1, gm1, gm2, norm, swf
76 : TYPE(cp_logger_type), POINTER :: logger
77 : TYPE(dimer_env_type), POINTER :: dimer_env
78 : TYPE(section_vals_type), POINTER :: print_section
79 :
80 1816 : NULLIFY (dimer_env, logger, print_section)
81 1816 : CALL timeset(routineN, handle)
82 1816 : logger => cp_get_default_logger()
83 :
84 1816 : CPASSERT(ASSOCIATED(gopt_env))
85 1816 : dimer_env => gopt_env%dimer_env
86 1816 : CPASSERT(ASSOCIATED(dimer_env))
87 1816 : iw = cp_print_key_unit_nr(logger, gopt_env%geo_section, "PRINT%PROGRAM_RUN_INFO", extension=".log")
88 : ! Possibly rotate Dimer or just compute Gradients of point 0 for Translation
89 1816 : IF (gopt_env%dimer_rotation) THEN
90 : IF (debug_this_module .AND. (iw > 0)) THEN
91 : WRITE (iw, '(A)') "NVEC:"
92 : WRITE (iw, '(3F15.9)') dimer_env%nvec
93 : END IF
94 2154 : SELECT CASE (dimer_env%rot%rotation_step)
95 : CASE (do_first_rotation_step, do_third_rotation_step)
96 726 : eval_analytical = .TRUE.
97 726 : IF ((dimer_env%rot%rotation_step == do_third_rotation_step) .AND. (dimer_env%rot%interpolate_gradient)) THEN
98 : eval_analytical = .FALSE.
99 : END IF
100 : IF (eval_analytical) THEN
101 : ! Compute energy, gradients and rotation vector for R1
102 132 : CALL cp_eval_at_ts_low(gopt_env, x, 1, dimer_env, calc_force, f1, dimer_env%rot%g1)
103 : ELSE
104 594 : angle1 = dimer_env%rot%angle1
105 594 : angle2 = dimer_env%rot%angle2
106 : dimer_env%rot%g1 = SIN(angle1 - angle2)/SIN(angle1)*dimer_env%rot%g1 + &
107 : SIN(angle2)/SIN(angle1)*dimer_env%rot%g1p + &
108 10932 : (1.0_dp - COS(angle2) - SIN(angle2)*TAN(angle1/2.0_dp))*dimer_env%rot%g0
109 : END IF
110 :
111 : ! Determine the theta vector (i.e. the search direction for line minimization)
112 24954 : gradient = -2.0_dp*(dimer_env%rot%g1 - dimer_env%rot%g0)
113 : IF (debug_this_module .AND. (iw > 0)) THEN
114 : WRITE (iw, '(A)') "G1 vector:"
115 : WRITE (iw, '(3F15.9)') dimer_env%rot%g1
116 : WRITE (iw, '(A)') "-2*(G1-G0) vector:"
117 : WRITE (iw, '(3F15.9)') gradient
118 : END IF
119 726 : CALL get_theta(gradient, dimer_env, norm)
120 726 : f = norm
121 726 : dimer_env%cg_rot%norm_theta_old = dimer_env%cg_rot%norm_theta
122 726 : dimer_env%cg_rot%norm_theta = norm
123 :
124 : IF (debug_this_module .AND. (iw > 0)) THEN
125 : WRITE (iw, '(A,F20.10)') "Rotational Force step (1,3): module:", norm
126 : WRITE (iw, '(3F15.9)') gradient
127 : END IF
128 :
129 : ! Compute curvature and derivative of the curvature w.r.t. the rotational angle
130 12840 : dimer_env%rot%curvature = DOT_PRODUCT(dimer_env%rot%g1 - dimer_env%rot%g0, dimer_env%nvec)/dimer_env%dr
131 12840 : dimer_env%rot%dCdp = 2.0_dp*DOT_PRODUCT(dimer_env%rot%g1 - dimer_env%rot%g0, gradient)/dimer_env%dr
132 :
133 726 : dimer_env%rot%rotation_step = do_second_rotation_step
134 12840 : gradient = -gradient
135 : CASE (do_second_rotation_step)
136 : ! Compute energy, gradients and rotation vector for R1
137 702 : CALL cp_eval_at_ts_low(gopt_env, x, 1, dimer_env, calc_force, f1, dimer_env%rot%g1p)
138 12708 : dimer_env%rot%curvature = DOT_PRODUCT(dimer_env%rot%g1p - dimer_env%rot%g0, dimer_env%nvec)/dimer_env%dr
139 702 : dimer_env%rot%rotation_step = do_third_rotation_step
140 :
141 : ! Determine the theta vector (i.e. the search direction for line minimization)
142 : ! This is never used for getting a new theta but is consistent in order to
143 : ! give back the right value of f
144 24714 : gradient = -2.0_dp*(dimer_env%rot%g1p - dimer_env%rot%g0)
145 702 : CALL get_theta(gradient, dimer_env, norm)
146 702 : f = norm
147 :
148 2856 : IF (debug_this_module .AND. (iw > 0)) THEN
149 : WRITE (iw, '(A)') "Rotational Force step (1,3):"
150 : WRITE (iw, '(3F15.9)') gradient
151 : END IF
152 : END SELECT
153 : ELSE
154 : ! Compute energy, gradients and rotation vector for R0
155 388 : CALL cp_eval_at_ts_low(gopt_env, x, 0, dimer_env, calc_force, f, dimer_env%rot%g0)
156 :
157 : ! The dimer is rotated only when we are out of the translation line search
158 388 : IF (.NOT. gopt_env%do_line_search) THEN
159 : IF (debug_this_module .AND. (iw > 0)) THEN
160 : WRITE (iw, '(A)') "Entering the rotation module"
161 : WRITE (iw, '(A)') "G0 vector:"
162 : WRITE (iw, '(3F15.9)') dimer_env%rot%g0
163 : END IF
164 : CALL cp_rot_opt(gopt_env%gopt_dimer_env, x, gopt_env%gopt_dimer_param, &
165 132 : gopt_env%gopt_dimer_env%geo_section)
166 132 : dimer_env%rot%rotation_step = do_first_rotation_step
167 : END IF
168 :
169 388 : print_section => section_vals_get_subs_vals(gopt_env%gopt_dimer_env%geo_section, "PRINT")
170 :
171 : ! Correcting gradients for Translation K-DIMER or STANDARD
172 388 : IF (dimer_env%kdimer) THEN
173 6 : CALL cite_reference(Henkelman2014)
174 : ! K-DIMER
175 6 : IF (iw > 0) THEN
176 3 : WRITE (iw, '(T2,A)') "DIMER| Correcting gradients for Translation with K-DIMER method"
177 : END IF
178 6 : swf = 1.0_dp + EXP(dimer_env%beta*dimer_env%rot%curvature)
179 6 : gm2 = 1.0_dp - (1.0_dp/swf)
180 6 : gm1 = (2.0_dp/swf) - 1.0_dp
181 : gradient = gm2*(dimer_env%rot%g0 - 2.0_dp*DOT_PRODUCT(dimer_env%rot%g0, dimer_env%nvec)*dimer_env%nvec) &
182 168 : - gm1*(DOT_PRODUCT(dimer_env%rot%g0, dimer_env%nvec)*dimer_env%nvec)
183 6 : CALL remove_rot_transl_component(gopt_env, gradient, print_section)
184 : IF (debug_this_module .AND. (iw > 0)) WRITE (iw, *) "K-DIMER", dimer_env%beta, dimer_env%rot%curvature, &
185 : dimer_env%rot%dCdp, gm1, gm2, swf
186 : ELSE
187 382 : IF (iw > 0) THEN
188 191 : WRITE (iw, '(T2,A)') "DIMER| Correcting gradients for Translation with standard method"
189 : END IF
190 382 : IF (dimer_env%rot%curvature > 0) THEN
191 4518 : gradient = -DOT_PRODUCT(dimer_env%rot%g0, dimer_env%nvec)*dimer_env%nvec
192 90 : CALL remove_rot_transl_component(gopt_env, gradient, print_section)
193 : ELSE
194 11578 : gradient = dimer_env%rot%g0 - 2.0_dp*DOT_PRODUCT(dimer_env%rot%g0, dimer_env%nvec)*dimer_env%nvec
195 292 : CALL remove_rot_transl_component(gopt_env, gradient, print_section)
196 : END IF
197 : END IF
198 : IF (debug_this_module .AND. (iw > 0)) THEN
199 : WRITE (iw, *) "final gradient:", gradient
200 : WRITE (iw, '(A,F20.10)') "norm gradient:", SQRT(DOT_PRODUCT(gradient, gradient))
201 : END IF
202 388 : IF (.NOT. gopt_env%do_line_search) THEN
203 1908 : f = SQRT(DOT_PRODUCT(gradient, gradient))
204 : ELSE
205 3772 : f = -DOT_PRODUCT(gradient, dimer_env%tsl%tls_vec)
206 : END IF
207 : END IF
208 1816 : CALL cp_print_key_finished_output(iw, logger, gopt_env%geo_section, "PRINT%PROGRAM_RUN_INFO")
209 1816 : CALL timestop(handle)
210 1816 : END SUBROUTINE cp_eval_at_ts
211 :
212 : ! **************************************************************************************************
213 : !> \brief This function removes translational forces after project of the gradient
214 : !> \param gopt_env ...
215 : !> \param gradient ...
216 : !> \param print_section ...
217 : !> \par History
218 : !> 2016/03/02 [LTong] added fixed atom constraint for gradient
219 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
220 : ! **************************************************************************************************
221 388 : SUBROUTINE remove_rot_transl_component(gopt_env, gradient, print_section)
222 : TYPE(gopt_f_type), POINTER :: gopt_env
223 : REAL(KIND=dp), DIMENSION(:) :: gradient
224 : TYPE(section_vals_type), POINTER :: print_section
225 :
226 : CHARACTER(len=*), PARAMETER :: routineN = 'remove_rot_transl_component'
227 :
228 : INTEGER :: dof, handle, i, j, natoms
229 : REAL(KIND=dp) :: norm, norm_gradient_old
230 388 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: D, mat
231 : TYPE(cp_subsys_type), POINTER :: subsys
232 : TYPE(particle_list_type), POINTER :: particles
233 :
234 388 : CALL timeset(routineN, handle)
235 388 : NULLIFY (mat)
236 388 : CALL force_env_get(gopt_env%force_env, subsys=subsys)
237 388 : CALL cp_subsys_get(subsys, particles=particles)
238 :
239 388 : natoms = particles%n_els
240 5680 : norm_gradient_old = SQRT(DOT_PRODUCT(gradient, gradient))
241 388 : IF (norm_gradient_old > 0.0_dp) THEN
242 388 : IF (natoms > 1) THEN
243 : CALL rot_ana(particles%els, mat, dof, print_section, keep_rotations=.FALSE., &
244 350 : mass_weighted=.FALSE., natoms=natoms)
245 :
246 : ! Orthogonalize gradient with respect to the full set of Roto-Trasl vectors
247 1400 : ALLOCATE (D(3*natoms, dof))
248 : ! Check First orthogonality in the first element of the basis set
249 2406 : DO i = 1, dof
250 63664 : D(:, i) = mat(:, i)
251 7436 : DO j = i + 1, dof
252 81380 : norm = DOT_PRODUCT(mat(:, i), mat(:, j))
253 7086 : CPASSERT(ABS(norm) < thrs_motion)
254 : END DO
255 : END DO
256 2406 : DO i = 1, dof
257 32860 : norm = DOT_PRODUCT(gradient, D(:, i))
258 33210 : gradient = gradient - norm*D(:, i)
259 : END DO
260 350 : DEALLOCATE (D)
261 350 : DEALLOCATE (mat)
262 : END IF
263 : END IF
264 : ! apply constraint
265 388 : CALL dimer_fixed_atom_control(gradient, subsys)
266 388 : CALL timestop(handle)
267 388 : END SUBROUTINE remove_rot_transl_component
268 :
269 : ! **************************************************************************************************
270 : !> \brief ...
271 : !> \param gopt_env ...
272 : !> \param x ...
273 : !> \param dimer_index ...
274 : !> \param dimer_env ...
275 : !> \param calc_force ...
276 : !> \param f ...
277 : !> \param gradient ...
278 : !> \par History
279 : !> none
280 : !> \author Luca Bellucci and Teodoro Laino - created [tlaino] - 01.2008
281 : ! **************************************************************************************************
282 2444 : SUBROUTINE cp_eval_at_ts_low(gopt_env, x, dimer_index, dimer_env, calc_force, &
283 1222 : f, gradient)
284 : TYPE(gopt_f_type), POINTER :: gopt_env
285 : REAL(KIND=dp), DIMENSION(:), POINTER :: x
286 : INTEGER, INTENT(IN) :: dimer_index
287 : TYPE(dimer_env_type), POINTER :: dimer_env
288 : LOGICAL, INTENT(IN) :: calc_force
289 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: f
290 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: gradient
291 :
292 : CHARACTER(len=*), PARAMETER :: routineN = 'cp_eval_at_ts_low'
293 :
294 : INTEGER :: handle, idg, idir, ip
295 : TYPE(cp_subsys_type), POINTER :: subsys
296 : TYPE(particle_list_type), POINTER :: particles
297 :
298 1222 : CALL timeset(routineN, handle)
299 1222 : idg = 0
300 1222 : CALL force_env_get(gopt_env%force_env, subsys=subsys)
301 1222 : CALL cp_subsys_get(subsys, particles=particles)
302 7580 : DO ip = 1, particles%n_els
303 26654 : DO idir = 1, 3
304 19074 : idg = idg + 1
305 25432 : particles%els(ip)%r(idir) = x(idg) + REAL(dimer_index, KIND=dp)*dimer_env%nvec(idg)*dimer_env%dr
306 : END DO
307 : END DO
308 :
309 : ! Compute energy and forces
310 1222 : CALL force_env_calc_energy_force(gopt_env%force_env, calc_force=calc_force)
311 :
312 : ! Possibly take the potential energy
313 1222 : IF (PRESENT(f)) THEN
314 1222 : CALL force_env_get(gopt_env%force_env, potential_energy=f)
315 : END IF
316 :
317 : ! Possibly take the gradients
318 1222 : IF (PRESENT(gradient)) THEN
319 1222 : idg = 0
320 1222 : CALL cp_subsys_get(subsys, particles=particles)
321 7580 : DO ip = 1, particles%n_els
322 26654 : DO idir = 1, 3
323 19074 : idg = idg + 1
324 19074 : CPASSERT(SIZE(gradient) >= idg)
325 25432 : gradient(idg) = -particles%els(ip)%f(idir)
326 : END DO
327 : END DO
328 : END IF
329 1222 : CALL timestop(handle)
330 1222 : END SUBROUTINE cp_eval_at_ts_low
331 :
332 : END MODULE dimer_methods
|
