(git:0de0cc2)
glbopt_worker.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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 ! **************************************************************************************************
13  USE cp_subsys_types, ONLY: cp_subsys_get,&
14  cp_subsys_type,&
17  USE f77_interface, ONLY: create_force_env,&
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,&
28  section_vals_type,&
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 
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  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  worker%root_section => root_section
99  worker%id = worker_id
100  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  mpi_comm=para_env)
109 
110  ! ======= More setup stuff =======
111  CALL f_env_add_defaults(worker%f_env_id, worker%f_env)
112  worker%force_env => worker%f_env%force_env
113  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  DO i = 1, worker_id
117  CALL worker%globenv%gaussian_rng_stream%reset_to_next_substream()
118  END DO
119 
120  CALL cp_subsys_get(worker%subsys, natom=worker%n_atoms)
121 
122  ! fetch original value from input
123  CALL section_vals_val_get(root_section, "MOTION%GEO_OPT%MAX_ITER", i_val=worker%gopt_max_iter)
124  glbopt_section => section_vals_get_subs_vals(root_section, "SWARM%GLOBAL_OPT")
125 
126  CALL section_vals_val_get(glbopt_section, "BUMP_STEPS_UPWARDS", i_val=worker%bump_steps_upwards)
127  CALL section_vals_val_get(glbopt_section, "BUMP_STEPS_DOWNWARDS", i_val=worker%bump_steps_downwards)
128  CALL section_vals_val_get(glbopt_section, "MD_BUMPS_MAX", i_val=worker%md_bumps_max)
129  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  worker%fragmentation_threshold = dist_in_angstrom/angstrom
132  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  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  CALL swarm_message_get(cmd, "command", command)
149  IF (trim(command) == "md_and_gopt") THEN
150  CALL run_mdgopt(worker, cmd, report)
151  ELSE
152  cpabort("Worker: received unknown command")
153  END IF
154 
155  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: positions
173  TYPE(glbopt_mdctrl_data_type), TARGET :: mdctrl_data
174  TYPE(mdctrl_type), POINTER :: mdctrl_p
175  TYPE(mdctrl_type), TARGET :: mdctrl
176 
177  NULLIFY (positions)
178 
179  CALL swarm_message_get(cmd, "temperature", temperature)
180  CALL swarm_message_get(cmd, "iframe", iframe)
181  IF (iframe > 1) THEN
182  CALL swarm_message_get(cmd, "positions", positions)
183  CALL unpack_subsys_particles(worker%subsys, r=positions)
184  END IF
185 
186  ! setup mdctrl callback
187  ALLOCATE (mdctrl_data%epot_history(worker%bump_steps_downwards + worker%bump_steps_upwards + 1))
188  mdctrl_data%epot_history = 0.0
189  mdctrl_data%md_bump_counter = 0
190  mdctrl_data%bump_steps_upwards = worker%bump_steps_upwards
191  mdctrl_data%bump_steps_downwards = worker%bump_steps_downwards
192  mdctrl_data%md_bumps_max = worker%md_bumps_max
193  mdctrl_data%output_unit = worker%iw
194  mdctrl%glbopt => mdctrl_data
195  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  prev_iframe = iframe
209  IF (iframe == 0) iframe = 1 ! qs_mol_dyn behaves differently for STEP_START_VAL=0
210  CALL section_vals_val_set(worker%root_section, "MOTION%MD%STEP_START_VAL", i_val=iframe - 1)
211  CALL section_vals_val_set(worker%root_section, "MOTION%MD%TEMPERATURE", r_val=temperature)
212 
213  IF (worker%iw > 0) THEN
214  WRITE (worker%iw, '(A,33X,F20.3)') ' GLBOPT| MD temperature [K]', temperature*kelvin
215  WRITE (worker%iw, '(A,29X,I10)') " GLBOPT| Starting MD at trajectory frame ", iframe
216  END IF
217 
218  ! run MD
219  CALL qs_mol_dyn(worker%force_env, worker%globenv, mdctrl=mdctrl_p)
220 
221  iframe = mdctrl_data%itimes + 1
222  md_steps = iframe - prev_iframe
223  IF (worker%iw > 0) WRITE (worker%iw, '(A,I4,A)') " GLBOPT| md ended after ", md_steps, " steps."
224 
225  ! fix fragmentation
226  IF (.NOT. ASSOCIATED(positions)) ALLOCATE (positions(3*worker%n_atoms))
227  CALL pack_subsys_particles(worker%subsys, r=positions)
228  n_fragments = 0
229  DO
230  n_fragments = n_fragments + 1
231  IF (fix_fragmentation(positions, worker%fragmentation_threshold)) EXIT
232  END DO
233  CALL unpack_subsys_particles(worker%subsys, r=positions)
234 
235  IF (n_fragments > 0 .AND. worker%iw > 0) &
236  WRITE (worker%iw, '(A,13X,I10)') " GLBOPT| Ran fix_fragmentation times:", n_fragments
237 
238  ! setup geometry optimization
239  IF (worker%iw > 0) WRITE (worker%iw, '(A,13X,I10)') " GLBOPT| Starting local optimisation at trajectory frame ", iframe
240  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  i_val=iframe + worker%gopt_max_iter)
243 
244  ! run geometry optimization
245  CALL cp_geo_opt(worker%force_env, worker%globenv, rm_restart_info=.false.)
246 
247  prev_iframe = iframe
248  CALL section_vals_val_get(worker%root_section, "MOTION%GEO_OPT%STEP_START_VAL", i_val=iframe)
249  iframe = iframe + 2 ! Compensates for different START_VAL interpretation.
250  gopt_steps = iframe - prev_iframe - 1
251  IF (worker%iw > 0) WRITE (worker%iw, '(A,I4,A)') " GLBOPT| gopt ended after ", gopt_steps, " steps."
252  CALL force_env_get(worker%force_env, potential_energy=epot)
253  IF (worker%iw > 0) WRITE (worker%iw, '(A,25X,E20.10)') ' GLBOPT| Potential Energy [Hartree]', epot
254 
255  ! assemble report
256  CALL swarm_message_add(report, "Epot", epot)
257  CALL swarm_message_add(report, "iframe", iframe)
258  CALL swarm_message_add(report, "md_steps", md_steps)
259  CALL swarm_message_add(report, "gopt_steps", gopt_steps)
260  CALL pack_subsys_particles(worker%subsys, r=positions)
261  CALL swarm_message_add(report, "positions", positions)
262 
263  DEALLOCATE (positions)
264  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  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: stack
281  LOGICAL, ALLOCATABLE, DIMENSION(:) :: marked
282  REAL(kind=dp) :: d, dr(3), min_dist, s
283 
284  n_particles = SIZE(positions)/3
285  ALLOCATE (stack(n_particles), marked(n_particles))
286 
287  marked = .false.; stack_size = 0
288 
289  ! First particle taken as root of flooding, mark it and push to stack
290  marked(1) = .true.; stack(1) = 1; stack_size = 1
291 
292  DO WHILE (stack_size > 0)
293  i = stack(stack_size); stack_size = stack_size - 1 !pop
294  DO j = 1, n_particles
295  IF (norm(diff(positions, i, j)) < 1.25*bondlength) THEN ! they are close = they are connected
296  IF (.NOT. marked(j)) THEN
297  marked(j) = .true.
298  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  all_connected = all(marked) !did we visit every particle?
305  IF (all_connected) RETURN
306 
307  ! make sure we keep the larger chunk
308  IF (count(marked) < n_particles/2) marked(:) = .NOT. (marked(:))
309 
310  min_dist = huge(1.0)
311  cluster_edge = -1
312  fragment_edge = -1
313  DO i = 1, n_particles
314  IF (marked(i)) cycle
315  DO j = 1, n_particles
316  IF (.NOT. marked(j)) cycle
317  d = norm(diff(positions, i, j))
318  IF (d < min_dist) THEN
319  min_dist = d
320  cluster_edge = i
321  fragment_edge = j
322  END IF
323  END DO
324  END DO
325 
326  dr = diff(positions, cluster_edge, fragment_edge)
327  s = 1.0 - bondlength/norm(dr)
328  DO i = 1, n_particles
329  IF (marked(i)) cycle
330  positions(3*i - 2:3*i) = positions(3*i - 2:3*i) - s*dr
331  END DO
332 
333  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  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  dr = positions(3*i - 2:3*i) - positions(3*j - 2:3*j)
349  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  PURE FUNCTION norm(vec) RESULT(res)
358  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: vec
359  REAL(kind=dp) :: res
360 
361  res = sqrt(dot_product(vec, vec))
362  END FUNCTION norm
363 
364 ! **************************************************************************************************
365 !> \brief Finalizes worker for global optimization
366 !> \param worker ...
367 !> \author Ole Schuett
368 ! **************************************************************************************************
369  SUBROUTINE glbopt_worker_finalize(worker)
370  TYPE(glbopt_worker_type), INTENT(INOUT) :: worker
371 
372  INTEGER :: ierr
373 
374  CALL f_env_rm_defaults(worker%f_env)
375  CALL destroy_force_env(worker%f_env_id, ierr)
376  IF (ierr /= 0) cpabort("destroy_force_env failed")
377  END SUBROUTINE glbopt_worker_finalize
378 
379 END MODULE glbopt_worker
types that represent a subsys, i.e. a part of the system
subroutine, public unpack_subsys_particles(subsys, f, r, s, v, fscale, cell)
Unpack components of a subsystem particle sets into a single vector.
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell)
returns information about various attributes of the given subsys
subroutine, public pack_subsys_particles(subsys, f, r, s, v, fscale, cell)
Pack components of a subsystem particle sets into a single vector.
interface to use cp2k as library
Definition: f77_interface.F:20
recursive subroutine, public destroy_force_env(env_id, ierr, q_finalize)
deallocates the force_env with the given id
subroutine, public f_env_add_defaults(f_env_id, f_env, handle)
adds the default environments of the f_env to the stack of the defaults, and returns a new error and ...
recursive subroutine, public create_force_env(new_env_id, input_declaration, input_path, output_path, mpi_comm, output_unit, owns_out_unit, input, ierr, work_dir, initial_variables)
creates a new force environment using the given input, and writing the output to the given output uni...
subroutine, public f_env_rm_defaults(f_env, ierr, handle)
removes the default environments of the f_env to the stack of the defaults, and sets ierr accordingly...
Interface for the force calculations.
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env)
returns various attributes about the force environment
performs geometry optimization
Definition: geo_opt.F:13
subroutine, public cp_geo_opt(force_env, globenv, eval_opt_geo, rm_restart_info)
Main driver to perform geometry optimization.
Definition: geo_opt.F:58
Worker routines used by global optimization schemes.
Definition: glbopt_worker.F:12
subroutine, public glbopt_worker_init(worker, input_declaration, para_env, root_section, input_path, worker_id, iw)
Initializes worker for global optimization.
Definition: glbopt_worker.F:87
subroutine, public glbopt_worker_finalize(worker)
Finalizes worker for global optimization.
subroutine, public glbopt_worker_execute(worker, cmd, report)
Central execute routine of global optimization worker.
Define type storing the global information of a run. Keep the amount of stored data small....
Definition: global_types.F:21
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
Perform a molecular dynamics (MD) run using QUICKSTEP.
Definition: md_run.F:14
subroutine, public qs_mol_dyn(force_env, globenv, averages, rm_restart_info, hmc_e_initial, hmc_e_final, mdctrl)
Main driver module for Molecular Dynamics.
Definition: md_run.F:121
A common interface for passing a callback into the md_run loop.
Definition: mdctrl_types.F:13
Interface to the message passing library MPI.
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public kelvin
Definition: physcon.F:165
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144
Swarm-message, a convenient data-container for with build-in serialization.
Definition: swarm_message.F:12