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