1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Master's routines for global optimization
8!> \author Ole Schuett
9! **************************************************************************************************
10MODULE glbopt_master
11   USE atomic_kind_types,               ONLY: atomic_kind_type,&
12                                              deallocate_atomic_kind_set
13   USE colvar_types,                    ONLY: colvar_p_release,&
14                                              colvar_p_type
15   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
16                                              cp_logger_type
17   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
18                                              cp_print_key_unit_nr
19   USE cp_para_types,                   ONLY: cp_para_env_type
20   USE cp_units,                        ONLY: cp_unit_from_cp2k
21   USE exclusion_types,                 ONLY: exclusion_release,&
22                                              exclusion_type
23   USE glbopt_mincrawl,                 ONLY: mincrawl_finalize,&
24                                              mincrawl_init,&
25                                              mincrawl_steer,&
26                                              mincrawl_type
27   USE glbopt_minhop,                   ONLY: minhop_finalize,&
28                                              minhop_init,&
29                                              minhop_steer,&
30                                              minhop_type
31   USE input_constants,                 ONLY: dump_xmol,&
32                                              glbopt_do_mincrawl,&
33                                              glbopt_do_minhop
34   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
35                                              section_vals_release,&
36                                              section_vals_retain,&
37                                              section_vals_type,&
38                                              section_vals_val_get
39   USE kinds,                           ONLY: default_string_length,&
40                                              dp,&
41                                              int_8
42   USE molecule_kind_types,             ONLY: deallocate_molecule_kind_set,&
43                                              molecule_kind_type
44   USE molecule_types,                  ONLY: deallocate_global_constraint,&
45                                              deallocate_molecule_set,&
46                                              global_constraint_type,&
47                                              molecule_type
48   USE particle_methods,                ONLY: write_particle_coordinates
49   USE particle_types,                  ONLY: deallocate_particle_set,&
50                                              particle_type
51   USE swarm_message,                   ONLY: swarm_message_get,&
52                                              swarm_message_type
53   USE topology,                        ONLY: topology_control
54#include "../base/base_uses.f90"
55
56   IMPLICIT NONE
57   PRIVATE
58
59   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'glbopt_master'
60
61   PUBLIC :: glbopt_master_type
62   PUBLIC :: glbopt_master_init, glbopt_master_finalize
63   PUBLIC :: glbopt_master_steer
64
65   TYPE glbopt_master_type
66      PRIVATE
67      REAL(KIND=dp)                                       :: E_lowest = HUGE(1.0_dp)
68      REAL(KIND=dp)                                       :: E_target = TINY(1.0_dp)
69      INTEGER                                             :: iw = 0
70      INTEGER                                             :: progress_traj_unit = 0
71      INTEGER(int_8)                                      :: total_md_steps = 0
72      INTEGER(int_8)                                      :: total_gopt_steps = 0
73      INTEGER(int_8)                                      :: count_reports = 0
74      INTEGER                                             :: method
75      TYPE(minhop_type), POINTER                           :: minhop
76      TYPE(mincrawl_type), POINTER                        :: mincrawl
77      INTEGER                                             :: i_iteration = 0
78      TYPE(atomic_kind_type), DIMENSION(:), POINTER       :: atomic_kind_set => Null()
79      TYPE(particle_type), DIMENSION(:), POINTER          :: particle_set => Null()
80      TYPE(section_vals_type), POINTER                    :: glbopt_section => Null()
81   END TYPE glbopt_master_type
82
83CONTAINS
84
85! **************************************************************************************************
86!> \brief Initializes the master of the global optimizer
87!> \param this ...
88!> \param para_env ...
89!> \param root_section ...
90!> \param n_walkers ...
91!> \param iw ...
92!> \author Ole Schuett
93! **************************************************************************************************
94   SUBROUTINE glbopt_master_init(this, para_env, root_section, n_walkers, iw)
95      TYPE(glbopt_master_type), INTENT(INOUT)            :: this
96      TYPE(cp_para_env_type), POINTER                    :: para_env
97      TYPE(section_vals_type), POINTER                   :: root_section
98      INTEGER, INTENT(IN)                                :: n_walkers, iw
99
100      CHARACTER(len=*), PARAMETER :: routineN = 'glbopt_master_init', &
101         routineP = moduleN//':'//routineN
102
103      TYPE(cp_logger_type), POINTER                      :: logger
104
105      NULLIFY (logger)
106
107      this%iw = iw
108
109      this%glbopt_section => section_vals_get_subs_vals(root_section, "SWARM%GLOBAL_OPT")
110      CALL section_vals_retain(this%glbopt_section)
111
112      CALL section_vals_val_get(this%glbopt_section, "E_TARGET", r_val=this%E_target)
113      CALL section_vals_val_get(this%glbopt_section, "METHOD", i_val=this%method)
114
115      CALL glbopt_read_particle_set(this, para_env, root_section)
116
117      logger => cp_get_default_logger()
118      this%progress_traj_unit = cp_print_key_unit_nr(logger, &
119                                                     this%glbopt_section, "PROGRESS_TRAJECTORY", &
120                                                     middle_name="progress", extension=".xyz", &
121                                                     file_action="WRITE", file_position="REWIND")
122
123      SELECT CASE (this%method)
124      CASE (glbopt_do_minhop)
125         ALLOCATE (this%minhop)
126         CALL minhop_init(this%minhop, this%glbopt_section, n_walkers, iw)
127      CASE (glbopt_do_mincrawl)
128         ALLOCATE (this%mincrawl)
129         CALL mincrawl_init(this%mincrawl, this%glbopt_section, n_walkers, iw, this%particle_set)
130      CASE DEFAULT
131         CPABORT("Unknown glbopt_method")
132      END SELECT
133   END SUBROUTINE glbopt_master_init
134
135! **************************************************************************************************
136!> \brief Helper-routine for glbopt_master_init, reads part of SUBSYS-section
137!> \param this ...
138!> \param para_env ...
139!> \param root_section ...
140!> \author Ole Schuett
141! **************************************************************************************************
142   SUBROUTINE glbopt_read_particle_set(this, para_env, root_section)
143      TYPE(glbopt_master_type), INTENT(INOUT)            :: this
144      TYPE(cp_para_env_type), POINTER                    :: para_env
145      TYPE(section_vals_type), POINTER                   :: root_section
146
147      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
148      TYPE(colvar_p_type), DIMENSION(:), POINTER         :: colvar_p
149      TYPE(exclusion_type), DIMENSION(:), POINTER        :: exclusions
150      TYPE(global_constraint_type), POINTER              :: gci
151      TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
152      TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
153      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
154      TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
155
156      NULLIFY (atomic_kind_set, particle_set, molecule_kind_set, molecule_set)
157      NULLIFY (colvar_p, gci, exclusions, force_env_section, subsys_section)
158
159      force_env_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
160      subsys_section => section_vals_get_subs_vals(root_section, "FORCE_EVAL%SUBSYS")
161
162      CALL topology_control(atomic_kind_set, &
163                            particle_set, &
164                            molecule_kind_set, &
165                            molecule_set, &
166                            colvar_p, &
167                            gci, &
168                            root_section=root_section, &
169                            para_env=para_env, &
170                            force_env_section=force_env_section, &
171                            subsys_section=subsys_section, &
172                            use_motion_section=.FALSE., &
173                            exclusions=exclusions)
174
175      !This is the only thing we need to write trajectories.
176      this%atomic_kind_set => atomic_kind_set
177      this%particle_set => particle_set
178      CALL exclusion_release(exclusions)
179      CALL deallocate_molecule_set(molecule_set)
180      CALL deallocate_molecule_kind_set(molecule_kind_set)
181      CALL deallocate_global_constraint(gci)
182      CALL colvar_p_release(colvar_p)
183
184   END SUBROUTINE glbopt_read_particle_set
185
186! **************************************************************************************************
187!> \brief Central steering routine of global optimizer
188!> \param this ...
189!> \param report ...
190!> \param cmd ...
191!> \param should_stop ...
192!> \author Ole Schuett
193! **************************************************************************************************
194   SUBROUTINE glbopt_master_steer(this, report, cmd, should_stop)
195      TYPE(glbopt_master_type), INTENT(INOUT)            :: this
196      TYPE(swarm_message_type)                           :: report, cmd
197      LOGICAL, INTENT(INOUT)                             :: should_stop
198
199      CHARACTER(len=*), PARAMETER :: routineN = 'glbopt_master_steer', &
200         routineP = moduleN//':'//routineN
201
202      CALL progress_report(this, report)
203
204      SELECT CASE (this%method)
205      CASE (glbopt_do_minhop)
206         CALL minhop_steer(this%minhop, report, cmd)
207      CASE (glbopt_do_mincrawl)
208         CALL mincrawl_steer(this%mincrawl, report, cmd)
209      CASE DEFAULT
210         CPABORT("Unknown glbopt_method")
211      END SELECT
212
213      IF (this%E_lowest < this%E_target) THEN
214         IF (this%iw > 0) WRITE (this%iw, "(A,I8,A)") &
215            " GLBOPT| Reached E_pot < E_target after ", this%i_iteration, " iterations. Quitting."
216         should_stop = .TRUE.
217      END IF
218   END SUBROUTINE glbopt_master_steer
219
220! **************************************************************************************************
221!> \brief Helper routine for glbopt_master_steer(), updates stats, etc.
222!> \param this ...
223!> \param report ...
224!> \author Ole Schuett
225! **************************************************************************************************
226   SUBROUTINE progress_report(this, report)
227      TYPE(glbopt_master_type), INTENT(INOUT)            :: this
228      TYPE(swarm_message_type)                           :: report
229
230      CHARACTER(len=default_string_length)               :: status
231      INTEGER                                            :: gopt_steps, md_steps, report_worker_id
232      REAL(KIND=dp)                                      :: report_Epot
233
234      this%i_iteration = this%i_iteration + 1
235
236      CALL swarm_message_get(report, "worker_id", report_worker_id)
237      CALL swarm_message_get(report, "status", status)
238
239      IF (TRIM(status) == "ok") THEN
240         CALL swarm_message_get(report, "Epot", report_Epot)
241         CALL swarm_message_get(report, "md_steps", md_steps)
242         CALL swarm_message_get(report, "gopt_steps", gopt_steps)
243         this%total_md_steps = this%total_md_steps + md_steps
244         this%total_gopt_steps = this%total_gopt_steps + gopt_steps
245         this%count_reports = this%count_reports + 1
246
247         IF (report_Epot < this%E_lowest) THEN
248            this%E_lowest = report_Epot
249            CALL write_progress_traj(this, report)
250         ENDIF
251
252         IF (this%iw > 0) THEN
253            WRITE (this%iw, '(A,46X,I8)') &
254               " GLBOPT| Reporting worker ", report_worker_id
255            WRITE (this%iw, '(A,20X,E15.8)') &
256               " GLBOPT| Reported potential energy [Hartree] ", report_Epot
257            WRITE (this%iw, '(A,13X,E15.8)') &
258               " GLBOPT| Lowest reported potential energy [Hartree] ", this%E_lowest
259            WRITE (this%iw, '(A,T71,F10.1)') &
260               " GLBOPT| Average number of MD steps", REAL(this%total_md_steps, KIND=dp)/this%count_reports
261            WRITE (this%iw, '(A,T71,F10.1)') &
262               " GLBOPT| Average number of Geo-Opt steps", REAL(this%total_gopt_steps, KIND=dp)/this%count_reports
263         END IF
264      END IF
265   END SUBROUTINE progress_report
266
267! **************************************************************************************************
268!> \brief Helper routine for progress_report(), write frames to trajectory.
269!> \param this ...
270!> \param report ...
271!> \author Ole Schuett
272! **************************************************************************************************
273   SUBROUTINE write_progress_traj(this, report)
274      TYPE(glbopt_master_type), INTENT(INOUT)            :: this
275      TYPE(swarm_message_type), INTENT(IN)               :: report
276
277      CHARACTER(len=default_string_length)               :: title, unit_str
278      INTEGER                                            :: report_worker_id
279      REAL(KIND=dp)                                      :: report_Epot, unit_conv
280      REAL(KIND=dp), DIMENSION(:), POINTER               :: report_positions
281
282      NULLIFY (report_positions)
283
284      IF (this%progress_traj_unit <= 0) RETURN
285
286      CALL swarm_message_get(report, "worker_id", report_worker_id)
287      CALL swarm_message_get(report, "positions", report_positions)
288      CALL swarm_message_get(report, "Epot", report_Epot)
289
290      WRITE (title, '(A,I8,A,I5,A,F20.10)') 'i = ', this%i_iteration, &
291         ' worker_id = ', report_worker_id, ' Epot = ', report_Epot
292
293      !get the conversion factor for the length unit
294      CALL section_vals_val_get(this%glbopt_section, "PROGRESS_TRAJECTORY%UNIT", &
295                                c_val=unit_str)
296      unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
297      CALL write_particle_coordinates(this%particle_set, &
298                                      iunit=this%progress_traj_unit, &
299                                      output_format=dump_xmol, &
300                                      content="POS", &
301                                      title=TRIM(title), &
302                                      array=report_positions, &
303                                      unit_conv=unit_conv)
304      DEALLOCATE (report_positions)
305   END SUBROUTINE write_progress_traj
306
307! **************************************************************************************************
308!> \brief Finalized the master of the global optimizer
309!> \param this ...
310!> \author Ole
311! **************************************************************************************************
312   SUBROUTINE glbopt_master_finalize(this)
313      TYPE(glbopt_master_type), INTENT(INOUT)            :: this
314
315      CHARACTER(len=*), PARAMETER :: routineN = 'glbopt_master_finalize', &
316         routineP = moduleN//':'//routineN
317
318      TYPE(cp_logger_type), POINTER                      :: logger
319
320      NULLIFY (logger)
321
322      SELECT CASE (this%method)
323      CASE (glbopt_do_minhop)
324         CALL minhop_finalize(this%minhop)
325         DEALLOCATE (this%minhop)
326      CASE (glbopt_do_mincrawl)
327         CALL mincrawl_finalize(this%mincrawl)
328         DEALLOCATE (this%mincrawl)
329      CASE DEFAULT
330         CPABORT("Unknown glbopt_method")
331      END SELECT
332
333      logger => cp_get_default_logger()
334      CALL cp_print_key_finished_output(this%progress_traj_unit, logger, &
335                                        this%glbopt_section, "PROGRESS_TRAJECTORY")
336
337      CALL section_vals_release(this%glbopt_section)
338      CALL deallocate_particle_set(this%particle_set)
339      CALL deallocate_atomic_kind_set(this%atomic_kind_set)
340
341   END SUBROUTINE glbopt_master_finalize
342
343END MODULE glbopt_master
344
345