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 Worker routines used by global optimization schemes
10 : !> \author Ole Schuett
11 : ! **************************************************************************************************
12 : MODULE glbopt_worker
13 : USE cp_subsys_types, ONLY: cp_subsys_get,&
14 : cp_subsys_type,&
15 : pack_subsys_particles,&
16 : unpack_subsys_particles
17 : USE f77_interface, ONLY: create_force_env,&
18 : destroy_force_env,&
19 : f_env_add_defaults,&
20 : f_env_rm_defaults,&
21 : f_env_type
22 : USE force_env_types, ONLY: force_env_get,&
23 : force_env_type
24 : USE geo_opt, ONLY: cp_geo_opt
25 : USE global_types, ONLY: global_environment_type
26 : USE input_section_types, ONLY: section_type,&
27 : section_vals_get_subs_vals,&
28 : section_vals_type,&
29 : section_vals_val_get,&
30 : section_vals_val_set
31 : USE kinds, ONLY: default_string_length,&
32 : dp
33 : USE md_run, ONLY: qs_mol_dyn
34 : USE mdctrl_types, ONLY: glbopt_mdctrl_data_type,&
35 : mdctrl_type
36 : USE message_passing, ONLY: mp_para_env_type
37 : USE physcon, ONLY: angstrom,&
38 : kelvin
39 : USE swarm_message, ONLY: swarm_message_add,&
40 : swarm_message_get,&
41 : swarm_message_type
42 : #include "../base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 : PRIVATE
46 :
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'glbopt_worker'
48 :
49 : PUBLIC :: glbopt_worker_init, glbopt_worker_finalize
50 : PUBLIC :: glbopt_worker_execute
51 : PUBLIC :: glbopt_worker_type
52 :
53 : TYPE glbopt_worker_type
54 : PRIVATE
55 : INTEGER :: id = -1
56 : INTEGER :: iw = -1
57 : INTEGER :: f_env_id = -1
58 : TYPE(f_env_type), POINTER :: f_env => NULL()
59 : TYPE(force_env_type), POINTER :: force_env => NULL()
60 : TYPE(cp_subsys_type), POINTER :: subsys => NULL()
61 : TYPE(section_vals_type), POINTER :: root_section => NULL()
62 : TYPE(global_environment_type), POINTER :: globenv => NULL()
63 : INTEGER :: gopt_max_iter = 0
64 : INTEGER :: bump_steps_downwards = 0
65 : INTEGER :: bump_steps_upwards = 0
66 : INTEGER :: md_bumps_max = 0
67 : REAL(KIND=dp) :: fragmentation_threshold = 0.0_dp
68 : INTEGER :: n_atoms = -1
69 : !REAL(KIND=dp) :: adaptive_timestep = 0.0
70 : END TYPE glbopt_worker_type
71 :
72 : CONTAINS
73 :
74 : ! **************************************************************************************************
75 : !> \brief Initializes worker for global optimization
76 : !> \param worker ...
77 : !> \param input_declaration ...
78 : !> \param para_env ...
79 : !> \param root_section ...
80 : !> \param input_path ...
81 : !> \param worker_id ...
82 : !> \param iw ...
83 : !> \author Ole Schuett
84 : ! **************************************************************************************************
85 6 : SUBROUTINE glbopt_worker_init(worker, input_declaration, para_env, root_section, &
86 : input_path, worker_id, iw)
87 : TYPE(glbopt_worker_type), INTENT(INOUT) :: worker
88 : TYPE(section_type), POINTER :: input_declaration
89 : TYPE(mp_para_env_type), POINTER :: para_env
90 : TYPE(section_vals_type), POINTER :: root_section
91 : CHARACTER(LEN=*), INTENT(IN) :: input_path
92 : INTEGER, INTENT(in) :: worker_id, iw
93 :
94 : INTEGER :: i
95 : REAL(kind=dp) :: dist_in_angstrom
96 : TYPE(section_vals_type), POINTER :: glbopt_section
97 :
98 3 : worker%root_section => root_section
99 3 : worker%id = worker_id
100 3 : worker%iw = iw
101 :
102 : ! ======= Create f_env =======
103 : CALL create_force_env(worker%f_env_id, &
104 : input_declaration=input_declaration, &
105 : input_path=input_path, &
106 : input=root_section, &
107 : output_unit=worker%iw, &
108 3 : mpi_comm=para_env)
109 :
110 : ! ======= More setup stuff =======
111 3 : CALL f_env_add_defaults(worker%f_env_id, worker%f_env)
112 3 : worker%force_env => worker%f_env%force_env
113 3 : CALL force_env_get(worker%force_env, globenv=worker%globenv, subsys=worker%subsys)
114 :
115 : ! We want different random-number-streams for each worker
116 6 : DO i = 1, worker_id
117 6 : CALL worker%globenv%gaussian_rng_stream%reset_to_next_substream()
118 : END DO
119 :
120 3 : CALL cp_subsys_get(worker%subsys, natom=worker%n_atoms)
121 :
122 : ! fetch original value from input
123 3 : CALL section_vals_val_get(root_section, "MOTION%GEO_OPT%MAX_ITER", i_val=worker%gopt_max_iter)
124 3 : glbopt_section => section_vals_get_subs_vals(root_section, "SWARM%GLOBAL_OPT")
125 :
126 3 : CALL section_vals_val_get(glbopt_section, "BUMP_STEPS_UPWARDS", i_val=worker%bump_steps_upwards)
127 3 : CALL section_vals_val_get(glbopt_section, "BUMP_STEPS_DOWNWARDS", i_val=worker%bump_steps_downwards)
128 3 : CALL section_vals_val_get(glbopt_section, "MD_BUMPS_MAX", i_val=worker%md_bumps_max)
129 3 : CALL section_vals_val_get(glbopt_section, "FRAGMENTATION_THRESHOLD", r_val=dist_in_angstrom)
130 : !CALL section_vals_val_get(glbopt_section,"MD_ADAPTIVE_TIMESTEP", r_val=worker%adaptive_timestep)
131 3 : worker%fragmentation_threshold = dist_in_angstrom/angstrom
132 3 : END SUBROUTINE glbopt_worker_init
133 :
134 : ! **************************************************************************************************
135 : !> \brief Central execute routine of global optimization worker
136 : !> \param worker ...
137 : !> \param cmd ...
138 : !> \param report ...
139 : !> \author Ole Schuett
140 : ! **************************************************************************************************
141 25 : SUBROUTINE glbopt_worker_execute(worker, cmd, report)
142 : TYPE(glbopt_worker_type), INTENT(INOUT) :: worker
143 : TYPE(swarm_message_type), INTENT(IN) :: cmd
144 : TYPE(swarm_message_type), INTENT(INOUT) :: report
145 :
146 : CHARACTER(len=default_string_length) :: command
147 :
148 25 : CALL swarm_message_get(cmd, "command", command)
149 25 : IF (TRIM(command) == "md_and_gopt") THEN
150 25 : CALL run_mdgopt(worker, cmd, report)
151 : ELSE
152 0 : CPABORT("Worker: received unknown command")
153 : END IF
154 :
155 25 : END SUBROUTINE glbopt_worker_execute
156 :
157 : ! **************************************************************************************************
158 : !> \brief Performs an escape attempt as need by e.g. Minima Hopping
159 : !> \param worker ...
160 : !> \param cmd ...
161 : !> \param report ...
162 : !> \author Ole Schuett
163 : ! **************************************************************************************************
164 25 : SUBROUTINE run_mdgopt(worker, cmd, report)
165 : TYPE(glbopt_worker_type), INTENT(INOUT) :: worker
166 : TYPE(swarm_message_type), INTENT(IN) :: cmd
167 : TYPE(swarm_message_type), INTENT(INOUT) :: report
168 :
169 : INTEGER :: gopt_steps, iframe, md_steps, &
170 : n_fragments, prev_iframe
171 : REAL(kind=dp) :: Epot, temperature
172 25 : REAL(KIND=dp), DIMENSION(:), POINTER :: positions
173 25 : TYPE(glbopt_mdctrl_data_type), TARGET :: mdctrl_data
174 : TYPE(mdctrl_type), POINTER :: mdctrl_p
175 : TYPE(mdctrl_type), TARGET :: mdctrl
176 :
177 25 : NULLIFY (positions)
178 :
179 25 : CALL swarm_message_get(cmd, "temperature", temperature)
180 25 : CALL swarm_message_get(cmd, "iframe", iframe)
181 25 : IF (iframe > 1) THEN
182 23 : CALL swarm_message_get(cmd, "positions", positions)
183 23 : CALL unpack_subsys_particles(worker%subsys, r=positions)
184 : END IF
185 :
186 : ! setup mdctrl callback
187 75 : ALLOCATE (mdctrl_data%epot_history(worker%bump_steps_downwards + worker%bump_steps_upwards + 1))
188 150 : mdctrl_data%epot_history = 0.0
189 25 : mdctrl_data%md_bump_counter = 0
190 25 : mdctrl_data%bump_steps_upwards = worker%bump_steps_upwards
191 25 : mdctrl_data%bump_steps_downwards = worker%bump_steps_downwards
192 25 : mdctrl_data%md_bumps_max = worker%md_bumps_max
193 25 : mdctrl_data%output_unit = worker%iw
194 25 : mdctrl%glbopt => mdctrl_data
195 25 : mdctrl_p => mdctrl
196 :
197 : !IF(worker%adaptive_timestep > 0.0) THEN
198 : ! !TODO: 300K is hard encoded
199 : ! boltz = 1.0 + exp( -temperature * kelvin / 150.0 )
200 : ! timestep = 4.0 * ( boltz - 1.0 ) / boltz / femtoseconds
201 : ! !timestep = 0.01_dp / femtoseconds
202 : ! !timestep = SQRT(MIN(0.5, 2.0/(1+exp(-300.0/(temperature*kelvin))))) / femtoseconds
203 : ! CALL section_vals_val_set(worker%root_section, "MOTION%MD%TIMESTEP", r_val=timestep)
204 : ! IF(worker%iw>0)&
205 : ! WRITE (worker%iw,'(A,35X,F20.3)') ' GLBOPT| MD timestep [fs]',timestep*femtoseconds
206 : !ENDIF
207 :
208 25 : prev_iframe = iframe
209 25 : IF (iframe == 0) iframe = 1 ! qs_mol_dyn behaves differently for STEP_START_VAL=0
210 25 : CALL section_vals_val_set(worker%root_section, "MOTION%MD%STEP_START_VAL", i_val=iframe - 1)
211 25 : CALL section_vals_val_set(worker%root_section, "MOTION%MD%TEMPERATURE", r_val=temperature)
212 :
213 25 : IF (worker%iw > 0) THEN
214 25 : WRITE (worker%iw, '(A,33X,F20.3)') ' GLBOPT| MD temperature [K]', temperature*kelvin
215 25 : WRITE (worker%iw, '(A,29X,I10)') " GLBOPT| Starting MD at trajectory frame ", iframe
216 : END IF
217 :
218 : ! run MD
219 25 : CALL qs_mol_dyn(worker%force_env, worker%globenv, mdctrl=mdctrl_p)
220 :
221 25 : iframe = mdctrl_data%itimes + 1
222 25 : md_steps = iframe - prev_iframe
223 25 : IF (worker%iw > 0) WRITE (worker%iw, '(A,I4,A)') " GLBOPT| md ended after ", md_steps, " steps."
224 :
225 : ! fix fragmentation
226 29 : IF (.NOT. ASSOCIATED(positions)) ALLOCATE (positions(3*worker%n_atoms))
227 25 : CALL pack_subsys_particles(worker%subsys, r=positions)
228 25 : n_fragments = 0
229 : DO
230 25 : n_fragments = n_fragments + 1
231 25 : IF (fix_fragmentation(positions, worker%fragmentation_threshold)) EXIT
232 : END DO
233 25 : CALL unpack_subsys_particles(worker%subsys, r=positions)
234 :
235 25 : IF (n_fragments > 0 .AND. worker%iw > 0) &
236 25 : WRITE (worker%iw, '(A,13X,I10)') " GLBOPT| Ran fix_fragmentation times:", n_fragments
237 :
238 : ! setup geometry optimization
239 25 : IF (worker%iw > 0) WRITE (worker%iw, '(A,13X,I10)') " GLBOPT| Starting local optimisation at trajectory frame ", iframe
240 25 : CALL section_vals_val_set(worker%root_section, "MOTION%GEO_OPT%STEP_START_VAL", i_val=iframe - 1)
241 : CALL section_vals_val_set(worker%root_section, "MOTION%GEO_OPT%MAX_ITER", &
242 25 : i_val=iframe + worker%gopt_max_iter)
243 :
244 : ! run geometry optimization
245 25 : CALL cp_geo_opt(worker%force_env, worker%globenv, rm_restart_info=.FALSE.)
246 :
247 25 : prev_iframe = iframe
248 25 : CALL section_vals_val_get(worker%root_section, "MOTION%GEO_OPT%STEP_START_VAL", i_val=iframe)
249 25 : iframe = iframe + 2 ! Compensates for different START_VAL interpretation.
250 25 : gopt_steps = iframe - prev_iframe - 1
251 25 : IF (worker%iw > 0) WRITE (worker%iw, '(A,I4,A)') " GLBOPT| gopt ended after ", gopt_steps, " steps."
252 25 : CALL force_env_get(worker%force_env, potential_energy=Epot)
253 25 : IF (worker%iw > 0) WRITE (worker%iw, '(A,25X,E20.10)') ' GLBOPT| Potential Energy [Hartree]', Epot
254 :
255 : ! assemble report
256 25 : CALL swarm_message_add(report, "Epot", Epot)
257 25 : CALL swarm_message_add(report, "iframe", iframe)
258 25 : CALL swarm_message_add(report, "md_steps", md_steps)
259 25 : CALL swarm_message_add(report, "gopt_steps", gopt_steps)
260 25 : CALL pack_subsys_particles(worker%subsys, r=positions)
261 25 : CALL swarm_message_add(report, "positions", positions)
262 :
263 25 : DEALLOCATE (positions)
264 100 : END SUBROUTINE run_mdgopt
265 :
266 : ! **************************************************************************************************
267 : !> \brief Helper routine for run_mdgopt, fixes a fragmented atomic cluster.
268 : !> \param positions ...
269 : !> \param bondlength ...
270 : !> \return ...
271 : !> \author Stefan Goedecker
272 : ! **************************************************************************************************
273 25 : FUNCTION fix_fragmentation(positions, bondlength) RESULT(all_connected)
274 : REAL(KIND=dp), DIMENSION(:) :: positions
275 : REAL(KIND=dp) :: bondlength
276 : LOGICAL :: all_connected
277 :
278 : INTEGER :: cluster_edge, fragment_edge, i, j, &
279 : n_particles, stack_size
280 25 : INTEGER, ALLOCATABLE, DIMENSION(:) :: stack
281 25 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: marked
282 : REAL(KIND=dp) :: d, dr(3), min_dist, s
283 :
284 25 : n_particles = SIZE(positions)/3
285 100 : ALLOCATE (stack(n_particles), marked(n_particles))
286 :
287 275 : marked = .FALSE.; stack_size = 0
288 :
289 : ! First particle taken as root of flooding, mark it and push to stack
290 25 : marked(1) = .TRUE.; stack(1) = 1; stack_size = 1
291 :
292 275 : DO WHILE (stack_size > 0)
293 250 : i = stack(stack_size); stack_size = stack_size - 1 !pop
294 2775 : DO j = 1, n_particles
295 2750 : IF (norm(diff(positions, i, j)) < 1.25*bondlength) THEN ! they are close = they are connected
296 4960 : IF (.NOT. marked(j)) THEN
297 225 : marked(j) = .TRUE.
298 225 : stack(stack_size + 1) = j; stack_size = stack_size + 1; !push
299 : END IF
300 : END IF
301 : END DO
302 : END DO
303 :
304 275 : all_connected = ALL(marked) !did we visit every particle?
305 25 : IF (all_connected) RETURN
306 :
307 : ! make sure we keep the larger chunk
308 0 : IF (COUNT(marked) < n_particles/2) marked(:) = .NOT. (marked(:))
309 :
310 0 : min_dist = HUGE(1.0)
311 0 : cluster_edge = -1
312 0 : fragment_edge = -1
313 0 : DO i = 1, n_particles
314 0 : IF (marked(i)) CYCLE
315 0 : DO j = 1, n_particles
316 0 : IF (.NOT. marked(j)) CYCLE
317 0 : d = norm(diff(positions, i, j))
318 0 : IF (d < min_dist) THEN
319 0 : min_dist = d
320 0 : cluster_edge = i
321 0 : fragment_edge = j
322 : END IF
323 : END DO
324 : END DO
325 :
326 0 : dr = diff(positions, cluster_edge, fragment_edge)
327 0 : s = 1.0 - bondlength/norm(dr)
328 0 : DO i = 1, n_particles
329 0 : IF (marked(i)) CYCLE
330 0 : positions(3*i - 2:3*i) = positions(3*i - 2:3*i) - s*dr
331 : END DO
332 :
333 25 : END FUNCTION fix_fragmentation
334 :
335 : ! **************************************************************************************************
336 : !> \brief Helper routine for fix_fragmentation, calculates atomic distance
337 : !> \param positions ...
338 : !> \param i ...
339 : !> \param j ...
340 : !> \return ...
341 : !> \author Ole Schuett
342 : ! **************************************************************************************************
343 2500 : PURE FUNCTION diff(positions, i, j) RESULT(dr)
344 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: positions
345 : INTEGER, INTENT(IN) :: i, j
346 : REAL(KIND=dp), DIMENSION(3) :: dr
347 :
348 10000 : dr = positions(3*i - 2:3*i) - positions(3*j - 2:3*j)
349 2500 : END FUNCTION diff
350 :
351 : ! **************************************************************************************************
352 : !> \brief Helper routine for fix_fragmentation, calculates vector norm
353 : !> \param vec ...
354 : !> \return ...
355 : !> \author Ole Schuett
356 : ! **************************************************************************************************
357 2500 : PURE FUNCTION norm(vec) RESULT(res)
358 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vec
359 : REAL(KIND=dp) :: res
360 :
361 10000 : res = SQRT(DOT_PRODUCT(vec, vec))
362 2500 : END FUNCTION norm
363 :
364 : ! **************************************************************************************************
365 : !> \brief Finalizes worker for global optimization
366 : !> \param worker ...
367 : !> \author Ole Schuett
368 : ! **************************************************************************************************
369 6 : SUBROUTINE glbopt_worker_finalize(worker)
370 : TYPE(glbopt_worker_type), INTENT(INOUT) :: worker
371 :
372 : INTEGER :: ierr
373 :
374 3 : CALL f_env_rm_defaults(worker%f_env)
375 3 : CALL destroy_force_env(worker%f_env_id, ierr)
376 3 : IF (ierr /= 0) CPABORT("destroy_force_env failed")
377 3 : END SUBROUTINE glbopt_worker_finalize
378 :
379 0 : END MODULE glbopt_worker
|
