1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Routines for the Minima Hopping global optimization scheme
8!> \author Ole Schuett
9! **************************************************************************************************
10MODULE glbopt_minhop
11   USE bibliography,                    ONLY: Goedecker2004,&
12                                              cite_reference
13   USE glbopt_history,                  ONLY: history_add,&
14                                              history_finalize,&
15                                              history_fingerprint,&
16                                              history_fingerprint_match,&
17                                              history_fingerprint_type,&
18                                              history_init,&
19                                              history_lookup,&
20                                              history_type
21   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
22                                              section_vals_type,&
23                                              section_vals_val_get
24   USE kinds,                           ONLY: default_string_length,&
25                                              dp
26   USE physcon,                         ONLY: kelvin
27   USE swarm_message,                   ONLY: swarm_message_add,&
28                                              swarm_message_get,&
29                                              swarm_message_type
30#include "../base/base_uses.f90"
31
32   IMPLICIT NONE
33   PRIVATE
34
35   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'glbopt_minhop'
36
37   PUBLIC :: minhop_type
38   PUBLIC :: minhop_init, minhop_finalize
39   PUBLIC :: minhop_steer
40
41   TYPE worker_state_type
42      REAL(KIND=dp)                                       :: Eaccept = -1.0
43      REAL(KIND=dp)                                       :: temp = -1.0
44      REAL(KIND=dp), DIMENSION(:), ALLOCATABLE            :: pos
45      REAL(KIND=dp)                                       :: Epot = -1.0
46      TYPE(history_fingerprint_type)                      :: fp
47      REAL(KIND=dp), DIMENSION(:), ALLOCATABLE            :: pos_hop
48      REAL(KIND=dp)                                       :: Epot_hop = HUGE(1.0)
49      TYPE(history_fingerprint_type)                      :: fp_hop
50      INTEGER                                             :: minima_id = -1
51      INTEGER                                             :: iframe = 1
52   END TYPE worker_state_type
53
54   TYPE minima_state_type
55      REAL(KIND=dp)                                       :: Eaccept = -1.0
56      REAL(KIND=dp)                                       :: temp = -1.0
57      REAL(KIND=dp), DIMENSION(:), ALLOCATABLE            :: pos
58      REAL(KIND=dp)                                       :: Epot = -1.0
59      TYPE(history_fingerprint_type)                      :: fp
60      LOGICAL                                             :: disabled = .FALSE.
61      INTEGER                                             :: n_active = 0
62      INTEGER                                             :: n_sampled = 0
63   END TYPE minima_state_type
64
65   TYPE minhop_type
66      PRIVATE
67      TYPE(history_type), DIMENSION(:), ALLOCATABLE        :: history
68      TYPE(worker_state_type), DIMENSION(:), ALLOCATABLE  :: worker_state
69      TYPE(minima_state_type), DIMENSION(:), ALLOCATABLE  :: minima_state
70      INTEGER                                             :: n_minima = 0
71      REAL(KIND=dp)                                       :: beta1 = 0
72      REAL(KIND=dp)                                       :: beta2 = 0
73      REAL(KIND=dp)                                       :: beta3 = 0
74      REAL(KIND=dp)                                       :: Eaccept0 = 0
75      REAL(KIND=dp)                                       :: temp_init = 0
76      REAL(KIND=dp)                                       :: temp_max = 0
77      REAL(KIND=dp)                                       :: temp_min = 0
78      REAL(KIND=dp)                                       :: alpha1 = 0
79      REAL(KIND=dp)                                       :: alpha2 = 0
80      INTEGER                                             :: n_accepted = 0
81      INTEGER                                             :: n_rejected = 0
82      INTEGER                                             :: iw = 0
83      INTEGER                                             :: n_workers = 0
84      LOGICAL                                             :: share_history = .FALSE.
85   END TYPE minhop_type
86
87CONTAINS
88
89! **************************************************************************************************
90!> \brief Initializes master for Minima Hopping
91!> \param this ...
92!> \param glbopt_section ...
93!> \param n_workers ...
94!> \param iw ...
95!> \author Ole Schuett
96! **************************************************************************************************
97   SUBROUTINE minhop_init(this, glbopt_section, n_workers, iw)
98      TYPE(minhop_type)                                  :: this
99      TYPE(section_vals_type), POINTER                   :: glbopt_section
100      INTEGER, INTENT(IN)                                :: n_workers, iw
101
102      INTEGER                                            :: i, n_histories
103      REAL(kind=dp)                                      :: temp_in_kelvin
104      TYPE(section_vals_type), POINTER                   :: history_section, minhop_section
105
106      CALL cite_reference(Goedecker2004)
107
108      ! read input
109      minhop_section => section_vals_get_subs_vals(glbopt_section, "MINIMA_HOPPING")
110      CALL section_vals_val_get(minhop_section, "BETA_1", r_val=this%beta1)
111      CALL section_vals_val_get(minhop_section, "BETA_2", r_val=this%beta2)
112      CALL section_vals_val_get(minhop_section, "BETA_3", r_val=this%beta3)
113      CALL section_vals_val_get(minhop_section, "ALPHA_1", r_val=this%alpha1)
114      CALL section_vals_val_get(minhop_section, "ALPHA_2", r_val=this%alpha2)
115      CALL section_vals_val_get(minhop_section, "E_ACCEPT_INIT", r_val=this%Eaccept0)
116      CALL section_vals_val_get(minhop_section, "TEMPERATURE_INIT", r_val=temp_in_kelvin)
117      this%temp_init = temp_in_kelvin/kelvin
118      CALL section_vals_val_get(minhop_section, "SHARE_HISTORY", l_val=this%share_history)
119
120      ! allocate history / histories
121      history_section => section_vals_get_subs_vals(glbopt_section, "HISTORY")
122      n_histories = n_workers
123      IF (this%share_history) n_histories = 1
124      ALLOCATE (this%history(n_histories))
125
126      !only the first history shall write to iw
127      CALL history_init(this%history(1), history_section, iw=iw)
128      DO i = 2, n_histories
129         CALL history_init(this%history(i), history_section, iw=-1)
130      END DO
131
132      ALLOCATE (this%worker_state(n_workers))
133      this%n_workers = n_workers
134      this%iw = iw
135
136      IF (this%iw > 0) THEN
137         WRITE (this%iw, '(A,T71,F10.3)') " MINHOP| beta_1", this%beta1
138         WRITE (this%iw, '(A,T71,F10.3)') " MINHOP| beta_2", this%beta2
139         WRITE (this%iw, '(A,T71,F10.3)') " MINHOP| beta_3", this%beta3
140         WRITE (this%iw, '(A,T71,F10.3)') " MINHOP| alpha_1", this%alpha1
141         WRITE (this%iw, '(A,T71,F10.3)') " MINHOP| alpha_2", this%alpha2
142         WRITE (this%iw, '(A,T71,F10.3)') " MINHOP| Initial acceptance energy [Hartree]", this%Eaccept0
143         WRITE (this%iw, '(A,T71,F10.3)') " MINHOP| Initial temperature [Kelvin]", this%temp_init*kelvin
144         WRITE (this%iw, '(A,T71,L10)') " MINHOP| All workers share a single history", this%share_history
145      ENDIF
146   END SUBROUTINE minhop_init
147
148! **************************************************************************************************
149!> \brief Central steering routine of Minima Hopping
150!> \param this ...
151!> \param report ...
152!> \param cmd ...
153!> \author Ole Schuett
154! **************************************************************************************************
155   SUBROUTINE minhop_steer(this, report, cmd)
156      TYPE(minhop_type)                                  :: this
157      TYPE(swarm_message_type)                           :: report, cmd
158
159      CHARACTER(len=default_string_length)               :: status
160      INTEGER                                            :: hid, iframe, wid
161      LOGICAL                                            :: minima_known
162      REAL(KIND=dp)                                      :: report_Epot
163      REAL(KIND=dp), DIMENSION(:), POINTER               :: report_positions
164      TYPE(history_fingerprint_type)                     :: report_fp
165
166      NULLIFY (report_positions)
167      CALL swarm_message_get(report, "worker_id", wid)
168      CALL swarm_message_get(report, "status", status)
169
170      IF (TRIM(status) == "initial_hello") THEN
171         this%worker_state(wid)%temp = this%temp_init
172         this%worker_state(wid)%Eaccept = this%Eaccept0
173         CALL swarm_message_add(cmd, "command", "md_and_gopt")
174         CALL swarm_message_add(cmd, "iframe", 1)
175         CALL swarm_message_add(cmd, "temperature", this%worker_state(wid)%temp)
176         IF (this%iw > 0) WRITE (this%iw, '(1X,A,1X,I10,1X,A,7X,F10.3)') &
177            "MINHOP| Sending worker", wid, &
178            "initial temperature [Kelvin]", this%worker_state(wid)%temp*kelvin
179         RETURN
180      ENDIF
181
182      hid = wid ! history_id = worker_id unless ....
183      IF (this%share_history) hid = 1 !...there is ONE shared history.
184
185      CALL swarm_message_get(report, "Epot", report_Epot)
186      CALL swarm_message_get(report, "positions", report_positions)
187
188      report_fp = history_fingerprint(report_Epot, report_positions)
189
190      IF (.NOT. ALLOCATED(this%worker_state(wid)%pos)) THEN
191         !init (first real report)
192         this%worker_state(wid)%Epot = report_Epot
193         ALLOCATE (this%worker_state(wid)%pos(SIZE(report_positions)))
194         this%worker_state(wid)%pos(:) = report_positions
195         this%worker_state(wid)%fp = report_fp
196      END IF
197
198      IF (history_fingerprint_match(this%history(hid), this%worker_state(wid)%fp, report_fp)) THEN
199         ! not escaped
200         IF (this%iw > 0) WRITE (this%iw, '(A)') " MINHOP| Not escaped"
201         this%worker_state(wid)%temp = this%worker_state(wid)%temp*this%beta1 !increasing temperature
202      ELSE
203         ! escaped
204         CALL history_lookup(this%history(hid), report_fp, minima_known)
205         IF (minima_known) THEN
206            IF (this%iw > 0) WRITE (this%iw, '(A)') " MINHOP| Escaped, old minima"
207            this%worker_state(wid)%temp = this%worker_state(wid)%temp*this%beta2 !increasing temperature
208         ELSE
209            IF (this%iw > 0) WRITE (this%iw, '(A)') " MINHOP| Escaped, new minima"
210            this%worker_state(wid)%temp = this%worker_state(wid)%temp*this%beta3 !decreasing temperature
211            CALL history_add(this%history(hid), report_fp)
212         ENDIF
213
214         IF (report_Epot < this%worker_state(wid)%Epot_hop) THEN
215            ! new locally lowest
216            IF (this%iw > 0) WRITE (this%iw, '(A)') " MINHOP| New locally lowest"
217            this%worker_state(wid)%Epot_hop = report_Epot
218            IF (.NOT. ALLOCATED(this%worker_state(wid)%pos_hop)) &
219               ALLOCATE (this%worker_state(wid)%pos_hop(SIZE(report_positions)))
220            this%worker_state(wid)%pos_hop(:) = report_positions
221            this%worker_state(wid)%fp_hop = report_fp
222         ENDIF
223
224         IF (this%worker_state(wid)%Epot_hop - this%worker_state(wid)%Epot < this%worker_state(wid)%Eaccept) THEN
225            ! accept
226            IF (this%iw > 0) WRITE (this%iw, '(A)') " MINHOP| Accept"
227            this%worker_state(wid)%Epot = this%worker_state(wid)%Epot_hop
228            this%worker_state(wid)%pos(:) = this%worker_state(wid)%pos_hop
229            this%worker_state(wid)%fp = this%worker_state(wid)%fp_hop
230            this%worker_state(wid)%Epot_hop = HUGE(1.0)
231
232            this%worker_state(wid)%Eaccept = this%worker_state(wid)%Eaccept*this%alpha1 !decreasing Eaccept
233            this%n_accepted = this%n_accepted + 1
234         ELSE
235            ! not accept
236            IF (this%iw > 0) WRITE (this%iw, '(A)') " MINHOP| Reject"
237            this%worker_state(wid)%Eaccept = this%worker_state(wid)%Eaccept*this%alpha2 !increasing Eaccept
238            this%n_rejected = this%n_rejected + 1
239         ENDIF
240      END IF
241
242      IF (this%iw > 0) THEN
243         WRITE (this%iw, '(A,15X,E20.10)') &
244            " MINHOP| Worker's acceptance Energy [Hartree]", this%worker_state(wid)%Eaccept
245         WRITE (this%iw, '(A,22X,F20.3)') &
246            " MINHOP| Worker's temperature [Kelvin]", this%worker_state(wid)%temp*kelvin
247      END IF
248
249      CALL swarm_message_get(report, "iframe", iframe)
250      CALL swarm_message_add(cmd, "iframe", iframe)
251      CALL swarm_message_add(cmd, "command", "md_and_gopt")
252      CALL swarm_message_add(cmd, "positions", this%worker_state(wid)%pos)
253      CALL swarm_message_add(cmd, "temperature", this%worker_state(wid)%temp)
254
255      IF (this%iw > 0) THEN
256         WRITE (this%iw, '(A,30X,I10)') &
257            " MINHOP| Total number of accepted minima", this%n_accepted
258         WRITE (this%iw, '(A,30X,I10)') &
259            " MINHOP| Total number of rejected minima", this%n_rejected
260      ENDIF
261
262      DEALLOCATE (report_positions)
263   END SUBROUTINE minhop_steer
264
265! **************************************************************************************************
266!> \brief Finalizes master for Minima Hopping
267!> \param this ...
268!> \author Ole Schuett
269! **************************************************************************************************
270   SUBROUTINE minhop_finalize(this)
271      TYPE(minhop_type)                                  :: this
272
273      INTEGER                                            :: i
274
275      DO i = 1, SIZE(this%history)
276         CALL history_finalize(this%history(i))
277      END DO
278   END SUBROUTINE minhop_finalize
279
280END MODULE glbopt_minhop
281
282