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