Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 : END TYPE glbopt_worker_type
70 :
71 : CONTAINS
72 :
73 : ! **************************************************************************************************
74 : !> \brief Initializes worker for global optimization
75 : !> \param worker ...
76 : !> \param input_declaration ...
77 : !> \param para_env ...
78 : !> \param root_section ...
79 : !> \param input_path ...
80 : !> \param worker_id ...
81 : !> \param iw ...
82 : !> \author Ole Schuett
83 : ! **************************************************************************************************
84 6 : SUBROUTINE glbopt_worker_init(worker, input_declaration, para_env, root_section, &
85 : input_path, worker_id, iw)
86 : TYPE(glbopt_worker_type), INTENT(INOUT) :: worker
87 : TYPE(section_type), POINTER :: input_declaration
88 : TYPE(mp_para_env_type), POINTER :: para_env
89 : TYPE(section_vals_type), POINTER :: root_section
90 : CHARACTER(LEN=*), INTENT(IN) :: input_path
91 : INTEGER, INTENT(in) :: worker_id, iw
92 :
93 : INTEGER :: i
94 : REAL(kind=dp) :: dist_in_angstrom
95 : TYPE(section_vals_type), POINTER :: glbopt_section
96 :
97 3 : worker%root_section => root_section
98 3 : worker%id = worker_id
99 3 : worker%iw = iw
100 :
101 : ! ======= Create f_env =======
102 : CALL create_force_env(worker%f_env_id, &
103 : input_declaration=input_declaration, &
104 : input_path=input_path, &
105 : input=root_section, &
106 : output_unit=worker%iw, &
107 3 : mpi_comm=para_env)
108 :
109 : ! ======= More setup stuff =======
110 3 : CALL f_env_add_defaults(worker%f_env_id, worker%f_env)
111 3 : worker%force_env => worker%f_env%force_env
112 3 : CALL force_env_get(worker%force_env, globenv=worker%globenv, subsys=worker%subsys)
113 :
114 : ! We want different random-number-streams for each worker
115 6 : DO i = 1, worker_id
116 6 : CALL worker%globenv%gaussian_rng_stream%reset_to_next_substream()
117 : END DO
118 :
119 3 : CALL cp_subsys_get(worker%subsys, natom=worker%n_atoms)
120 :
121 : ! fetch original value from input
122 3 : CALL section_vals_val_get(root_section, "MOTION%GEO_OPT%MAX_ITER", i_val=worker%gopt_max_iter)
123 3 : glbopt_section => section_vals_get_subs_vals(root_section, "SWARM%GLOBAL_OPT")
124 :
125 3 : CALL section_vals_val_get(glbopt_section, "BUMP_STEPS_UPWARDS", i_val=worker%bump_steps_upwards)
126 3 : CALL section_vals_val_get(glbopt_section, "BUMP_STEPS_DOWNWARDS", i_val=worker%bump_steps_downwards)
127 3 : CALL section_vals_val_get(glbopt_section, "MD_BUMPS_MAX", i_val=worker%md_bumps_max)
128 3 : CALL section_vals_val_get(glbopt_section, "FRAGMENTATION_THRESHOLD", r_val=dist_in_angstrom)
129 3 : worker%fragmentation_threshold = dist_in_angstrom/angstrom
130 3 : END SUBROUTINE glbopt_worker_init
131 :
132 : ! **************************************************************************************************
133 : !> \brief Central execute routine of global optimization worker
134 : !> \param worker ...
135 : !> \param cmd ...
136 : !> \param report ...
137 : !> \author Ole Schuett
138 : ! **************************************************************************************************
139 13 : SUBROUTINE glbopt_worker_execute(worker, cmd, report)
140 : TYPE(glbopt_worker_type), INTENT(INOUT) :: worker
141 : TYPE(swarm_message_type), INTENT(IN) :: cmd
142 : TYPE(swarm_message_type), INTENT(INOUT) :: report
143 :
144 : CHARACTER(len=default_string_length) :: command
145 :
146 13 : CALL swarm_message_get(cmd, "command", command)
147 13 : IF (TRIM(command) == "md_and_gopt") THEN
148 13 : CALL run_mdgopt(worker, cmd, report)
149 : ELSE
150 0 : CPABORT("Worker: received unknown command")
151 : END IF
152 :
153 13 : END SUBROUTINE glbopt_worker_execute
154 :
155 : ! **************************************************************************************************
156 : !> \brief Performs an escape attempt as need by e.g. Minima Hopping
157 : !> \param worker ...
158 : !> \param cmd ...
159 : !> \param report ...
160 : !> \author Ole Schuett
161 : ! **************************************************************************************************
162 13 : SUBROUTINE run_mdgopt(worker, cmd, report)
163 : TYPE(glbopt_worker_type), INTENT(INOUT) :: worker
164 : TYPE(swarm_message_type), INTENT(IN) :: cmd
165 : TYPE(swarm_message_type), INTENT(INOUT) :: report
166 :
167 : INTEGER :: gopt_steps, iframe, md_steps, &
168 : n_fragments, prev_iframe
169 : REAL(kind=dp) :: Epot, temperature
170 13 : REAL(KIND=dp), DIMENSION(:), POINTER :: positions
171 13 : TYPE(glbopt_mdctrl_data_type), TARGET :: mdctrl_data
172 : TYPE(mdctrl_type), POINTER :: mdctrl_p
173 : TYPE(mdctrl_type), TARGET :: mdctrl
174 :
175 13 : NULLIFY (positions)
176 :
177 13 : CALL swarm_message_get(cmd, "temperature", temperature)
178 13 : CALL swarm_message_get(cmd, "iframe", iframe)
179 13 : IF (iframe > 1) THEN
180 11 : CALL swarm_message_get(cmd, "positions", positions)
181 11 : CALL unpack_subsys_particles(worker%subsys, r=positions)
182 : END IF
183 :
184 : ! setup mdctrl callback
185 39 : ALLOCATE (mdctrl_data%epot_history(worker%bump_steps_downwards + worker%bump_steps_upwards + 1))
186 78 : mdctrl_data%epot_history = 0.0_dp
187 13 : mdctrl_data%md_bump_counter = 0
188 13 : mdctrl_data%bump_steps_upwards = worker%bump_steps_upwards
189 13 : mdctrl_data%bump_steps_downwards = worker%bump_steps_downwards
190 13 : mdctrl_data%md_bumps_max = worker%md_bumps_max
191 13 : mdctrl_data%output_unit = worker%iw
192 13 : mdctrl%glbopt => mdctrl_data
193 13 : mdctrl_p => mdctrl
194 :
195 13 : prev_iframe = iframe
196 13 : IF (iframe == 0) iframe = 1 ! qs_mol_dyn behaves differently for STEP_START_VAL=0
197 13 : CALL section_vals_val_set(worker%root_section, "MOTION%MD%STEP_START_VAL", i_val=iframe - 1)
198 13 : CALL section_vals_val_set(worker%root_section, "MOTION%MD%TEMPERATURE", r_val=temperature)
199 :
200 13 : IF (worker%iw > 0) THEN
201 13 : WRITE (worker%iw, '(A,33X,F20.3)') ' GLBOPT| MD temperature [K]', temperature*kelvin
202 13 : WRITE (worker%iw, '(A,29X,I10)') " GLBOPT| Starting MD at trajectory frame ", iframe
203 : END IF
204 :
205 : ! run MD
206 13 : CALL qs_mol_dyn(worker%force_env, worker%globenv, mdctrl=mdctrl_p)
207 :
208 13 : iframe = mdctrl_data%itimes + 1
209 13 : md_steps = iframe - prev_iframe
210 13 : IF (worker%iw > 0) WRITE (worker%iw, '(A,I4,A)') " GLBOPT| md ended after ", md_steps, " steps."
211 :
212 : ! fix fragmentation
213 17 : IF (.NOT. ASSOCIATED(positions)) ALLOCATE (positions(3*worker%n_atoms))
214 13 : CALL pack_subsys_particles(worker%subsys, r=positions)
215 13 : n_fragments = 0
216 : DO
217 13 : n_fragments = n_fragments + 1
218 13 : IF (fix_fragmentation(positions, worker%fragmentation_threshold)) EXIT
219 : END DO
220 13 : CALL unpack_subsys_particles(worker%subsys, r=positions)
221 :
222 13 : IF (n_fragments > 0 .AND. worker%iw > 0) THEN
223 13 : WRITE (worker%iw, '(A,13X,I10)') " GLBOPT| Ran fix_fragmentation times:", n_fragments
224 : END IF
225 :
226 : ! setup geometry optimization
227 13 : IF (worker%iw > 0) WRITE (worker%iw, '(A,13X,I10)') " GLBOPT| Starting local optimisation at trajectory frame ", iframe
228 13 : CALL section_vals_val_set(worker%root_section, "MOTION%GEO_OPT%STEP_START_VAL", i_val=iframe - 1)
229 : CALL section_vals_val_set(worker%root_section, "MOTION%GEO_OPT%MAX_ITER", &
230 13 : i_val=iframe + worker%gopt_max_iter)
231 :
232 : ! run geometry optimization
233 13 : CALL cp_geo_opt(worker%force_env, worker%globenv, rm_restart_info=.FALSE.)
234 :
235 13 : prev_iframe = iframe
236 13 : CALL section_vals_val_get(worker%root_section, "MOTION%GEO_OPT%STEP_START_VAL", i_val=iframe)
237 13 : iframe = iframe + 2 ! Compensates for different START_VAL interpretation.
238 13 : gopt_steps = iframe - prev_iframe - 1
239 13 : IF (worker%iw > 0) WRITE (worker%iw, '(A,I4,A)') " GLBOPT| gopt ended after ", gopt_steps, " steps."
240 13 : CALL force_env_get(worker%force_env, potential_energy=Epot)
241 13 : IF (worker%iw > 0) WRITE (worker%iw, '(A,25X,E20.10)') ' GLBOPT| Potential Energy [Hartree]', Epot
242 :
243 : ! assemble report
244 13 : CALL swarm_message_add(report, "Epot", Epot)
245 13 : CALL swarm_message_add(report, "iframe", iframe)
246 13 : CALL swarm_message_add(report, "md_steps", md_steps)
247 13 : CALL swarm_message_add(report, "gopt_steps", gopt_steps)
248 13 : CALL pack_subsys_particles(worker%subsys, r=positions)
249 13 : CALL swarm_message_add(report, "positions", positions)
250 :
251 13 : DEALLOCATE (positions)
252 52 : END SUBROUTINE run_mdgopt
253 :
254 : ! **************************************************************************************************
255 : !> \brief Helper routine for run_mdgopt, fixes a fragmented atomic cluster.
256 : !> \param positions ...
257 : !> \param bondlength ...
258 : !> \return ...
259 : !> \author Stefan Goedecker
260 : ! **************************************************************************************************
261 13 : FUNCTION fix_fragmentation(positions, bondlength) RESULT(all_connected)
262 : REAL(KIND=dp), DIMENSION(:) :: positions
263 : REAL(KIND=dp) :: bondlength
264 : LOGICAL :: all_connected
265 :
266 : REAL(KIND=dp), PARAMETER :: bond_thre = 1.25_dp
267 :
268 : INTEGER :: cluster_edge, fragment_edge, i, j, &
269 : n_particles, stack_size
270 13 : INTEGER, ALLOCATABLE, DIMENSION(:) :: stack
271 13 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: marked
272 : REAL(KIND=dp) :: d, dr(3), min_dist, s
273 :
274 13 : n_particles = SIZE(positions)/3
275 52 : ALLOCATE (stack(n_particles), marked(n_particles))
276 :
277 13 : marked = .FALSE.; stack_size = 0
278 :
279 : ! First particle taken as root of flooding, mark it and push to stack
280 13 : marked(1) = .TRUE.; stack(1) = 1; stack_size = 1
281 :
282 143 : DO WHILE (stack_size > 0)
283 130 : i = stack(stack_size); stack_size = stack_size - 1 !pop
284 1443 : DO j = 1, n_particles
285 5330 : IF (NORM2(diff(positions, i, j)) < bond_thre*bondlength) THEN ! they are close = they are connected
286 2590 : IF (.NOT. marked(j)) THEN
287 117 : marked(j) = .TRUE.
288 117 : stack(stack_size + 1) = j; stack_size = stack_size + 1 !push
289 : END IF
290 : END IF
291 : END DO
292 : END DO
293 :
294 143 : all_connected = ALL(marked) !did we visit every particle?
295 13 : IF (all_connected) RETURN
296 :
297 : ! make sure we keep the larger chunk
298 0 : IF (COUNT(marked) < n_particles/2) marked(:) = .NOT. (marked(:))
299 :
300 0 : min_dist = HUGE(1.0_dp)
301 0 : cluster_edge = -1
302 0 : fragment_edge = -1
303 0 : DO i = 1, n_particles
304 0 : IF (marked(i)) CYCLE
305 0 : DO j = 1, n_particles
306 0 : IF (.NOT. marked(j)) CYCLE
307 0 : d = NORM2(diff(positions, i, j))
308 0 : IF (d < min_dist) THEN
309 0 : min_dist = d
310 0 : cluster_edge = i
311 0 : fragment_edge = j
312 : END IF
313 : END DO
314 : END DO
315 :
316 0 : dr = diff(positions, cluster_edge, fragment_edge)
317 0 : s = 1.0_dp - bondlength/NORM2(dr)
318 0 : DO i = 1, n_particles
319 0 : IF (marked(i)) CYCLE
320 0 : positions(3*i - 2:3*i) = positions(3*i - 2:3*i) - s*dr
321 : END DO
322 :
323 13 : END FUNCTION fix_fragmentation
324 :
325 : ! **************************************************************************************************
326 : !> \brief Helper routine for fix_fragmentation, calculates atomic distance
327 : !> \param positions ...
328 : !> \param i ...
329 : !> \param j ...
330 : !> \return ...
331 : !> \author Ole Schuett
332 : ! **************************************************************************************************
333 1300 : PURE FUNCTION diff(positions, i, j) RESULT(dr)
334 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: positions
335 : INTEGER, INTENT(IN) :: i, j
336 : REAL(KIND=dp), DIMENSION(3) :: dr
337 :
338 5200 : dr = positions(3*i - 2:3*i) - positions(3*j - 2:3*j)
339 1300 : END FUNCTION diff
340 :
341 : ! **************************************************************************************************
342 : !> \brief Finalizes worker for global optimization
343 : !> \param worker ...
344 : !> \author Ole Schuett
345 : ! **************************************************************************************************
346 6 : SUBROUTINE glbopt_worker_finalize(worker)
347 : TYPE(glbopt_worker_type), INTENT(INOUT) :: worker
348 :
349 : INTEGER :: ierr
350 :
351 3 : CALL f_env_rm_defaults(worker%f_env)
352 3 : CALL destroy_force_env(worker%f_env_id, ierr)
353 3 : IF (ierr /= 0) CPABORT("destroy_force_env failed")
354 3 : END SUBROUTINE glbopt_worker_finalize
355 :
356 0 : END MODULE glbopt_worker
|
