Line data Source code
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 : MODULE optimize_input
8 : USE cell_types, ONLY: parse_cell_line
9 : USE cp2k_info, ONLY: write_restart_header
10 : USE cp_external_control, ONLY: external_control
11 : USE cp_log_handling, ONLY: cp_get_default_logger,&
12 : cp_logger_get_default_io_unit,&
13 : cp_logger_type
14 : USE cp_output_handling, ONLY: cp_add_iter_level,&
15 : cp_iterate,&
16 : cp_print_key_finished_output,&
17 : cp_print_key_unit_nr,&
18 : cp_rm_iter_level
19 : USE cp_parser_methods, ONLY: parser_read_line
20 : USE cp_parser_types, ONLY: cp_parser_type,&
21 : parser_create,&
22 : parser_release
23 : USE environment, ONLY: cp2k_get_walltime
24 : USE f77_interface, ONLY: calc_force,&
25 : create_force_env,&
26 : destroy_force_env,&
27 : set_cell
28 : USE input_constants, ONLY: opt_force_matching
29 : USE input_section_types, ONLY: section_type,&
30 : section_vals_get,&
31 : section_vals_get_subs_vals,&
32 : section_vals_type,&
33 : section_vals_val_get,&
34 : section_vals_val_set,&
35 : section_vals_write
36 : USE kinds, ONLY: default_path_length,&
37 : default_string_length,&
38 : dp
39 : USE machine, ONLY: m_flush,&
40 : m_walltime
41 : USE memory_utilities, ONLY: reallocate
42 : USE message_passing, ONLY: mp_comm_type,&
43 : mp_para_env_type
44 : USE parallel_rng_types, ONLY: UNIFORM,&
45 : rng_stream_type
46 : USE physcon, ONLY: bohr
47 : USE powell, ONLY: opt_state_type,&
48 : powell_optimize
49 : #include "./base/base_uses.f90"
50 :
51 : IMPLICIT NONE
52 : PRIVATE
53 :
54 : PUBLIC :: run_optimize_input
55 :
56 : TYPE fm_env_type
57 : CHARACTER(LEN=default_path_length) :: optimize_file_name = ""
58 :
59 : CHARACTER(LEN=default_path_length) :: ref_traj_file_name = ""
60 : CHARACTER(LEN=default_path_length) :: ref_force_file_name = ""
61 : CHARACTER(LEN=default_path_length) :: ref_cell_file_name = ""
62 :
63 : INTEGER :: group_size = -1
64 :
65 : REAL(KIND=dp) :: energy_weight = -1.0_dp
66 : REAL(KIND=dp) :: shift_mm = -1.0_dp
67 : REAL(KIND=dp) :: shift_qm = -1.0_dp
68 : LOGICAL :: shift_average = .FALSE.
69 : INTEGER :: frame_start = -1, frame_stop = -1, frame_stride = -1, frame_count = -1
70 : END TYPE
71 :
72 : TYPE variable_type
73 : CHARACTER(LEN=default_string_length) :: label = ""
74 : REAL(KIND=dp) :: value = -1.0_dp
75 : LOGICAL :: fixed = .FALSE.
76 : END TYPE
77 :
78 : TYPE oi_env_type
79 : INTEGER :: method = -1
80 : INTEGER :: seed = -1
81 : CHARACTER(LEN=default_path_length) :: project_name = ""
82 : TYPE(fm_env_type) :: fm_env = fm_env_type()
83 : TYPE(variable_type), DIMENSION(:), ALLOCATABLE :: variables
84 : REAL(KIND=dp) :: rhobeg = -1.0_dp, rhoend = -1.0_dp
85 : INTEGER :: maxfun = -1
86 : INTEGER :: iter_start_val = -1
87 : REAL(KIND=dp) :: randomize_variables = -1.0_dp
88 : REAL(KIND=dp) :: start_time = -1.0_dp, target_time = -1.0_dp
89 : END TYPE
90 :
91 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'optimize_input'
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief main entry point for methods aimed at optimizing parameters in a CP2K input file
97 : !> \param input_declaration ...
98 : !> \param root_section ...
99 : !> \param para_env ...
100 : !> \author Joost VandeVondele
101 : ! **************************************************************************************************
102 6 : SUBROUTINE run_optimize_input(input_declaration, root_section, para_env)
103 : TYPE(section_type), POINTER :: input_declaration
104 : TYPE(section_vals_type), POINTER :: root_section
105 : TYPE(mp_para_env_type), POINTER :: para_env
106 :
107 : CHARACTER(len=*), PARAMETER :: routineN = 'run_optimize_input'
108 :
109 : INTEGER :: handle, i_var
110 : REAL(KIND=dp) :: random_number, seed(3, 2)
111 6 : TYPE(oi_env_type) :: oi_env
112 6 : TYPE(rng_stream_type), ALLOCATABLE :: rng_stream
113 :
114 6 : CALL timeset(routineN, handle)
115 :
116 6 : oi_env%start_time = m_walltime()
117 :
118 6 : CALL parse_input(oi_env, root_section)
119 :
120 : ! if we have been asked to randomize the variables, we do this.
121 6 : IF (oi_env%randomize_variables .NE. 0.0_dp) THEN
122 36 : seed = REAL(oi_env%seed, KIND=dp)
123 4 : rng_stream = rng_stream_type("run_optimize_input", distribution_type=UNIFORM, seed=seed)
124 12 : DO i_var = 1, SIZE(oi_env%variables, 1)
125 12 : IF (.NOT. oi_env%variables(i_var)%fixed) THEN
126 : ! change with a random percentage the variable
127 8 : random_number = rng_stream%next()
128 : oi_env%variables(i_var)%value = oi_env%variables(i_var)%value* &
129 8 : (1.0_dp + (2*random_number - 1.0_dp)*oi_env%randomize_variables/100.0_dp)
130 : END IF
131 : END DO
132 : END IF
133 :
134 : ! proceed to actual methods
135 12 : SELECT CASE (oi_env%method)
136 : CASE (opt_force_matching)
137 6 : CALL force_matching(oi_env, input_declaration, root_section, para_env)
138 : CASE DEFAULT
139 6 : CPABORT("")
140 : END SELECT
141 :
142 6 : CALL timestop(handle)
143 :
144 6 : END SUBROUTINE run_optimize_input
145 :
146 : ! **************************************************************************************************
147 : !> \brief optimizes parameters by force/energy matching results against reference values
148 : !> \param oi_env ...
149 : !> \param input_declaration ...
150 : !> \param root_section ...
151 : !> \param para_env ...
152 : !> \author Joost VandeVondele
153 : ! **************************************************************************************************
154 6 : SUBROUTINE force_matching(oi_env, input_declaration, root_section, para_env)
155 : TYPE(oi_env_type) :: oi_env
156 : TYPE(section_type), POINTER :: input_declaration
157 : TYPE(section_vals_type), POINTER :: root_section
158 : TYPE(mp_para_env_type), POINTER :: para_env
159 :
160 : CHARACTER(len=*), PARAMETER :: routineN = 'force_matching'
161 :
162 : CHARACTER(len=default_path_length) :: input_path, output_path
163 : CHARACTER(len=default_string_length), &
164 6 : ALLOCATABLE, DIMENSION(:, :) :: initial_variables
165 : INTEGER :: color, energies_unit, handle, history_unit, i_atom, i_el, i_frame, i_free_var, &
166 : i_var, ierr, mepos_master, mepos_minion, n_atom, n_el, n_frames, n_free_var, n_groups, &
167 : n_var, new_env_id, num_pe_master, output_unit, restart_unit, state
168 6 : INTEGER, ALLOCATABLE, DIMENSION(:) :: free_var_index
169 6 : INTEGER, DIMENSION(:), POINTER :: group_distribution
170 : LOGICAL :: should_stop
171 : REAL(KIND=dp) :: e1, e2, e3, e4, e_pot, energy_weight, &
172 : re, rf, shift_mm, shift_qm, t1, t2, &
173 : t3, t4, t5
174 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: force, free_vars, pos
175 6 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy_traj, energy_traj_read, energy_var
176 6 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cell_traj, cell_traj_read, force_traj, &
177 6 : force_traj_read, force_var, pos_traj, &
178 6 : pos_traj_read
179 : TYPE(cp_logger_type), POINTER :: logger
180 : TYPE(mp_comm_type) :: mpi_comm_master, mpi_comm_minion, &
181 : mpi_comm_minion_primus
182 : TYPE(opt_state_type) :: ostate
183 : TYPE(section_vals_type), POINTER :: oi_section, variable_section
184 :
185 6 : CALL timeset(routineN, handle)
186 :
187 6 : logger => cp_get_default_logger()
188 6 : CALL cp_add_iter_level(logger%iter_info, "POWELL_OPT")
189 6 : output_unit = cp_logger_get_default_io_unit(logger)
190 :
191 6 : IF (output_unit > 0) THEN
192 3 : WRITE (output_unit, '(T2,A)') 'FORCE_MATCHING| good morning....'
193 : END IF
194 :
195 : ! do IO of ref traj / frc / cell
196 6 : NULLIFY (cell_traj)
197 6 : NULLIFY (cell_traj_read, force_traj_read, pos_traj_read, energy_traj_read)
198 6 : CALL read_reference_data(oi_env, para_env, force_traj_read, pos_traj_read, energy_traj_read, cell_traj_read)
199 6 : n_atom = SIZE(pos_traj_read, 2)
200 :
201 : ! adjust read data with respect to start/stop/stride
202 6 : IF (oi_env%fm_env%frame_stop < 0) oi_env%fm_env%frame_stop = SIZE(pos_traj_read, 3)
203 :
204 6 : IF (oi_env%fm_env%frame_count > 0) THEN
205 : oi_env%fm_env%frame_stride = (oi_env%fm_env%frame_stop - oi_env%fm_env%frame_start + 1 + &
206 0 : oi_env%fm_env%frame_count - 1)/(oi_env%fm_env%frame_count)
207 : END IF
208 6 : n_frames = (oi_env%fm_env%frame_stop - oi_env%fm_env%frame_start + oi_env%fm_env%frame_stride)/oi_env%fm_env%frame_stride
209 :
210 48 : ALLOCATE (force_traj(3, n_atom, n_frames), pos_traj(3, n_atom, n_frames), energy_traj(n_frames))
211 18 : IF (ASSOCIATED(cell_traj_read)) ALLOCATE (cell_traj(3, 3, n_frames))
212 :
213 6 : n_frames = 0
214 220 : DO i_frame = oi_env%fm_env%frame_start, oi_env%fm_env%frame_stop, oi_env%fm_env%frame_stride
215 214 : n_frames = n_frames + 1
216 13910 : force_traj(:, :, n_frames) = force_traj_read(:, :, i_frame)
217 13910 : pos_traj(:, :, n_frames) = pos_traj_read(:, :, i_frame)
218 214 : energy_traj(n_frames) = energy_traj_read(i_frame)
219 5356 : IF (ASSOCIATED(cell_traj)) cell_traj(:, :, n_frames) = cell_traj_read(:, :, i_frame)
220 : END DO
221 6 : DEALLOCATE (force_traj_read, pos_traj_read, energy_traj_read)
222 6 : IF (ASSOCIATED(cell_traj_read)) DEALLOCATE (cell_traj_read)
223 :
224 6 : n_el = 3*n_atom
225 24 : ALLOCATE (pos(n_el), force(n_el))
226 36 : ALLOCATE (energy_var(n_frames), force_var(3, n_atom, n_frames))
227 :
228 : ! split the para_env in a set of sub_para_envs that will do the force_env communications
229 6 : mpi_comm_master = para_env
230 6 : num_pe_master = para_env%num_pe
231 6 : mepos_master = para_env%mepos
232 18 : ALLOCATE (group_distribution(0:num_pe_master - 1))
233 6 : IF (oi_env%fm_env%group_size > para_env%num_pe) oi_env%fm_env%group_size = para_env%num_pe
234 :
235 6 : CALL mpi_comm_minion%from_split(mpi_comm_master, n_groups, group_distribution, subgroup_min_size=oi_env%fm_env%group_size)
236 6 : mepos_minion = mpi_comm_minion%mepos
237 6 : color = 0
238 6 : IF (mepos_minion == 0) color = 1
239 6 : CALL mpi_comm_minion_primus%from_split(mpi_comm_master, color)
240 :
241 : ! assign initial variables
242 6 : n_var = SIZE(oi_env%variables, 1)
243 18 : ALLOCATE (initial_variables(2, n_var))
244 6 : n_free_var = 0
245 18 : DO i_var = 1, n_var
246 12 : initial_variables(1, i_var) = oi_env%variables(i_var)%label
247 12 : WRITE (initial_variables(2, i_var), *) oi_env%variables(i_var)%value
248 18 : IF (.NOT. oi_env%variables(i_var)%fixed) n_free_var = n_free_var + 1
249 : END DO
250 30 : ALLOCATE (free_vars(n_free_var), free_var_index(n_free_var))
251 18 : i_free_var = 0
252 18 : DO i_var = 1, n_var
253 18 : IF (.NOT. oi_env%variables(i_var)%fixed) THEN
254 12 : i_free_var = i_free_var + 1
255 12 : free_var_index(i_free_var) = i_var
256 12 : free_vars(i_free_var) = oi_env%variables(free_var_index(i_free_var))%value
257 : END IF
258 : END DO
259 :
260 : ! create input and output file names.
261 6 : input_path = oi_env%fm_env%optimize_file_name
262 6 : WRITE (output_path, '(A,I0,A)') TRIM(oi_env%project_name)//"-worker-", group_distribution(mepos_master), ".out"
263 :
264 : ! initialize the powell optimizer
265 6 : energy_weight = oi_env%fm_env%energy_weight
266 6 : shift_mm = oi_env%fm_env%shift_mm
267 6 : shift_qm = oi_env%fm_env%shift_qm
268 :
269 6 : IF (para_env%is_source()) THEN
270 3 : ostate%nf = 0
271 3 : ostate%nvar = n_free_var
272 3 : ostate%rhoend = oi_env%rhoend
273 3 : ostate%rhobeg = oi_env%rhobeg
274 3 : ostate%maxfun = oi_env%maxfun
275 3 : ostate%iprint = 1
276 3 : ostate%unit = output_unit
277 3 : ostate%state = 0
278 : END IF
279 :
280 6 : IF (output_unit > 0) THEN
281 3 : WRITE (output_unit, '(T2,A,T60,I20)') 'FORCE_MATCHING| number of atoms per frame ', n_atom
282 3 : WRITE (output_unit, '(T2,A,T60,I20)') 'FORCE_MATCHING| number of frames ', n_frames
283 3 : WRITE (output_unit, '(T2,A,T60,I20)') 'FORCE_MATCHING| number of parallel groups ', n_groups
284 3 : WRITE (output_unit, '(T2,A,T60,I20)') 'FORCE_MATCHING| number of variables ', n_var
285 3 : WRITE (output_unit, '(T2,A,T60,I20)') 'FORCE_MATCHING| number of free variables ', n_free_var
286 3 : WRITE (output_unit, '(T2,A,A)') 'FORCE_MATCHING| optimize file name ', TRIM(input_path)
287 3 : WRITE (output_unit, '(T2,A,T60,F20.12)') 'FORCE_MATCHING| accuracy', ostate%rhoend
288 3 : WRITE (output_unit, '(T2,A,T60,F20.12)') 'FORCE_MATCHING| step size', ostate%rhobeg
289 3 : WRITE (output_unit, '(T2,A,T60,I20)') 'FORCE_MATCHING| max function evaluation', ostate%maxfun
290 3 : WRITE (output_unit, '(T2,A,T60,L20)') 'FORCE_MATCHING| shift average', oi_env%fm_env%shift_average
291 3 : WRITE (output_unit, '(T2,A)') 'FORCE_MATCHING| initial values:'
292 9 : DO i_var = 1, n_var
293 9 : WRITE (output_unit, '(T2,A,1X,E28.16)') TRIM(oi_env%variables(i_var)%label), oi_env%variables(i_var)%value
294 : END DO
295 3 : WRITE (output_unit, '(T2,A)') 'FORCE_MATCHING| switching to POWELL optimization of the free parameters'
296 3 : WRITE (output_unit, '()')
297 3 : WRITE (output_unit, '(T2,A20,A20,A11,A11)') 'iteration number', 'function value', 'time', 'time Force'
298 3 : CALL m_flush(output_unit)
299 : END IF
300 :
301 6 : t1 = m_walltime()
302 :
303 98 : DO
304 :
305 : ! globalize the state
306 104 : IF (para_env%is_source()) state = ostate%state
307 104 : CALL para_env%bcast(state)
308 :
309 : ! if required get the energy of this set of params
310 104 : IF (state == 2) THEN
311 :
312 86 : CALL cp_iterate(logger%iter_info, last=.FALSE.)
313 :
314 : ! create a new force env, updating the free vars as needed
315 258 : DO i_free_var = 1, n_free_var
316 172 : WRITE (initial_variables(2, free_var_index(i_free_var)), *) free_vars(i_free_var)
317 258 : oi_env%variables(free_var_index(i_free_var))%value = free_vars(i_free_var)
318 : END DO
319 :
320 : ierr = 0
321 : CALL create_force_env(new_env_id=new_env_id, input_declaration=input_declaration, &
322 : input_path=input_path, output_path=output_path, &
323 86 : mpi_comm=mpi_comm_minion, initial_variables=initial_variables, ierr=ierr)
324 :
325 : ! set to zero initialy, for easier mp_summing
326 3460 : energy_var = 0.0_dp
327 111428 : force_var = 0.0_dp
328 :
329 : ! compute energies and forces for all frames, doing the work on a minion sub group based on round robin
330 86 : t5 = 0.0_dp
331 3460 : DO i_frame = group_distribution(mepos_master) + 1, n_frames, n_groups
332 :
333 : ! set new cell if needed
334 3374 : IF (ASSOCIATED(cell_traj)) THEN
335 3374 : CALL set_cell(env_id=new_env_id, new_cell=cell_traj(:, :, i_frame), ierr=ierr)
336 : END IF
337 :
338 : ! copy pos from ref
339 : i_el = 0
340 30366 : DO i_atom = 1, n_atom
341 26992 : pos(i_el + 1) = pos_traj(1, i_atom, i_frame)
342 26992 : pos(i_el + 2) = pos_traj(2, i_atom, i_frame)
343 26992 : pos(i_el + 3) = pos_traj(3, i_atom, i_frame)
344 30366 : i_el = i_el + 3
345 : END DO
346 :
347 : ! evaluate energy/force with new pos
348 3374 : t3 = m_walltime()
349 3374 : CALL calc_force(env_id=new_env_id, pos=pos, n_el_pos=n_el, e_pot=e_pot, force=force, n_el_force=n_el, ierr=ierr)
350 3374 : t4 = m_walltime()
351 3374 : t5 = t5 + t4 - t3
352 :
353 : ! copy force and energy in place
354 3374 : energy_var(i_frame) = e_pot
355 3374 : i_el = 0
356 33826 : DO i_atom = 1, n_atom
357 26992 : force_var(1, i_atom, i_frame) = force(i_el + 1)
358 26992 : force_var(2, i_atom, i_frame) = force(i_el + 2)
359 26992 : force_var(3, i_atom, i_frame) = force(i_el + 3)
360 30366 : i_el = i_el + 3
361 : END DO
362 :
363 : END DO
364 :
365 : ! clean up force env, get ready for the next round
366 86 : CALL destroy_force_env(env_id=new_env_id, ierr=ierr)
367 :
368 : ! get data everywhere on the master group, we could reduce the amount of data by reducing to partial RMSD first
369 : ! furthermore, we should only do this operation among the masters of the minion group
370 86 : IF (mepos_minion == 0) THEN
371 3417 : CALL mpi_comm_minion_primus%sum(energy_var)
372 111385 : CALL mpi_comm_minion_primus%sum(force_var)
373 : END IF
374 :
375 : ! now evaluate the target function to be minimized (only valid on mepos_minion==0)
376 86 : IF (para_env%is_source()) THEN
377 55714 : rf = SQRT(SUM((force_var - force_traj)**2)/(REAL(n_frames, KIND=dp)*REAL(n_atom, KIND=dp)))
378 43 : IF (oi_env%fm_env%shift_average) THEN
379 0 : shift_mm = SUM(energy_var)/n_frames
380 0 : shift_qm = SUM(energy_traj)/n_frames
381 : END IF
382 1730 : re = SQRT(SUM(((energy_var - shift_mm) - (energy_traj - shift_qm))**2)/n_frames)
383 43 : ostate%f = energy_weight*re + rf
384 43 : t2 = m_walltime()
385 43 : WRITE (output_unit, '(T2,I20,F20.12,F11.3,F11.3)') oi_env%iter_start_val + ostate%nf, ostate%f, t2 - t1, t5
386 43 : CALL m_flush(output_unit)
387 43 : t1 = m_walltime()
388 : END IF
389 :
390 : ! the history file with the trajectory of the parameters
391 : history_unit = cp_print_key_unit_nr(logger, root_section, "OPTIMIZE_INPUT%HISTORY", &
392 86 : extension=".dat")
393 86 : IF (history_unit > 0) THEN
394 8 : WRITE (UNIT=history_unit, FMT="(I20,F20.12,1000F20.12)") oi_env%iter_start_val + ostate%nf, ostate%f, free_vars
395 : END IF
396 86 : CALL cp_print_key_finished_output(history_unit, logger, root_section, "OPTIMIZE_INPUT%HISTORY")
397 :
398 : ! the energy profile for all frames
399 : energies_unit = cp_print_key_unit_nr(logger, root_section, "OPTIMIZE_INPUT%FORCE_MATCHING%COMPARE_ENERGIES", &
400 86 : file_position="REWIND", extension=".dat")
401 86 : IF (energies_unit > 0) THEN
402 8 : WRITE (UNIT=energies_unit, FMT="(A20,A20,A20,A20)") "#frame", "ref", "fit", "diff"
403 324 : DO i_frame = 1, n_frames
404 316 : e1 = energy_traj(i_frame) - shift_qm
405 316 : e2 = energy_var(i_frame) - shift_mm
406 324 : WRITE (UNIT=energies_unit, FMT="(I20,F20.12,F20.12,F20.12)") i_frame, e1, e2, e1 - e2
407 : END DO
408 : END IF
409 86 : CALL cp_print_key_finished_output(energies_unit, logger, root_section, "OPTIMIZE_INPUT%FORCE_MATCHING%COMPARE_ENERGIES")
410 :
411 : ! the force profile for all frames
412 : energies_unit = cp_print_key_unit_nr(logger, root_section, "OPTIMIZE_INPUT%FORCE_MATCHING%COMPARE_FORCES", &
413 86 : file_position="REWIND", extension=".dat")
414 86 : IF (energies_unit > 0) THEN
415 43 : WRITE (UNIT=energies_unit, FMT="(A20,A20,A20,A20)") "#frame", "normalized diff", "diff", "ref", "ref sum"
416 1730 : DO i_frame = 1, n_frames
417 55671 : e1 = SQRT(SUM((force_var(:, :, i_frame) - force_traj(:, :, i_frame))**2))
418 55671 : e2 = SQRT(SUM((force_traj(:, :, i_frame))**2))
419 47236 : e3 = SQRT(SUM(SUM(force_traj(:, :, i_frame), DIM=2)**2))
420 47236 : e4 = SQRT(SUM(SUM(force_var(:, :, i_frame), DIM=2)**2))
421 1730 : WRITE (UNIT=energies_unit, FMT="(I20,F20.12,F20.12,F20.12,2F20.12)") i_frame, e1/e2, e1, e2, e3, e4
422 : END DO
423 : END IF
424 86 : CALL cp_print_key_finished_output(energies_unit, logger, root_section, "OPTIMIZE_INPUT%FORCE_MATCHING%COMPARE_FORCES")
425 :
426 : ! a restart file with the current values of the parameters
427 : restart_unit = cp_print_key_unit_nr(logger, root_section, "OPTIMIZE_INPUT%RESTART", extension=".restart", &
428 86 : file_position="REWIND", do_backup=.TRUE.)
429 86 : IF (restart_unit > 0) THEN
430 8 : oi_section => section_vals_get_subs_vals(root_section, "OPTIMIZE_INPUT")
431 8 : CALL section_vals_val_set(oi_section, "ITER_START_VAL", i_val=oi_env%iter_start_val + ostate%nf)
432 8 : variable_section => section_vals_get_subs_vals(oi_section, "VARIABLE")
433 24 : DO i_free_var = 1, n_free_var
434 : CALL section_vals_val_set(variable_section, "VALUE", i_rep_section=free_var_index(i_free_var), &
435 24 : r_val=free_vars(i_free_var))
436 : END DO
437 8 : CALL write_restart_header(restart_unit)
438 8 : CALL section_vals_write(root_section, unit_nr=restart_unit, hide_root=.TRUE.)
439 : END IF
440 86 : CALL cp_print_key_finished_output(restart_unit, logger, root_section, "OPTIMIZE_INPUT%RESTART")
441 :
442 : END IF
443 :
444 104 : IF (state == -1) EXIT
445 :
446 98 : CALL external_control(should_stop, "OPTIMIZE_INPUT", target_time=oi_env%target_time, start_time=oi_env%start_time)
447 :
448 98 : IF (should_stop) EXIT
449 :
450 : ! do a powell step if needed
451 98 : IF (para_env%is_source()) THEN
452 49 : CALL powell_optimize(ostate%nvar, free_vars, ostate)
453 : END IF
454 104 : CALL para_env%bcast(free_vars)
455 :
456 : END DO
457 :
458 : ! finally, get the best set of variables
459 6 : IF (para_env%is_source()) THEN
460 3 : ostate%state = 8
461 3 : CALL powell_optimize(ostate%nvar, free_vars, ostate)
462 : END IF
463 6 : CALL para_env%bcast(free_vars)
464 18 : DO i_free_var = 1, n_free_var
465 12 : WRITE (initial_variables(2, free_var_index(i_free_var)), *) free_vars(i_free_var)
466 18 : oi_env%variables(free_var_index(i_free_var))%value = free_vars(i_free_var)
467 : END DO
468 6 : IF (para_env%is_source()) THEN
469 3 : WRITE (output_unit, '(T2,A)') ''
470 3 : WRITE (output_unit, '(T2,A,T60,F20.12)') 'FORCE_MATCHING| optimal function value found so far:', ostate%fopt
471 3 : WRITE (output_unit, '(T2,A)') 'FORCE_MATCHING| optimal variables found so far:'
472 9 : DO i_var = 1, n_var
473 9 : WRITE (output_unit, '(T2,A,1X,E28.16)') TRIM(oi_env%variables(i_var)%label), oi_env%variables(i_var)%value
474 : END DO
475 : END IF
476 :
477 6 : CALL cp_rm_iter_level(logger%iter_info, "POWELL_OPT")
478 :
479 : ! deallocate for cleanup
480 6 : IF (ASSOCIATED(cell_traj)) DEALLOCATE (cell_traj)
481 6 : DEALLOCATE (pos, force, force_traj, pos_traj, force_var)
482 6 : DEALLOCATE (group_distribution, energy_traj, energy_var)
483 6 : CALL mpi_comm_minion%free()
484 6 : CALL mpi_comm_minion_primus%free()
485 :
486 6 : CALL timestop(handle)
487 :
488 12 : END SUBROUTINE force_matching
489 :
490 : ! **************************************************************************************************
491 : !> \brief reads the reference data for force matching results, the format of the files needs to be the CP2K xyz trajectory format
492 : !> \param oi_env ...
493 : !> \param para_env ...
494 : !> \param force_traj forces
495 : !> \param pos_traj position
496 : !> \param energy_traj energies, as extracted from the forces file
497 : !> \param cell_traj cell parameters, as extracted from a CP2K cell file
498 : !> \author Joost VandeVondele
499 : ! **************************************************************************************************
500 24 : SUBROUTINE read_reference_data(oi_env, para_env, force_traj, pos_traj, energy_traj, cell_traj)
501 : TYPE(oi_env_type) :: oi_env
502 : TYPE(mp_para_env_type), POINTER :: para_env
503 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: force_traj, pos_traj
504 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy_traj
505 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: cell_traj
506 :
507 : CHARACTER(len=*), PARAMETER :: routineN = 'read_reference_data'
508 :
509 : CHARACTER(len=default_path_length) :: filename
510 : CHARACTER(len=default_string_length) :: AA
511 : INTEGER :: cell_itimes, handle, i, iframe, &
512 : n_frames, n_frames_current, nread, &
513 : trj_itimes
514 : LOGICAL :: at_end, test_ok
515 : REAL(KIND=dp) :: cell_time, trj_epot, trj_time, vec(3), &
516 : vol
517 : TYPE(cp_parser_type) :: local_parser
518 :
519 6 : CALL timeset(routineN, handle)
520 :
521 : ! do IO of ref traj / frc / cell
522 :
523 : ! trajectory
524 6 : n_frames = 0
525 6 : n_frames_current = 0
526 6 : NULLIFY (pos_traj, energy_traj, force_traj)
527 6 : filename = oi_env%fm_env%ref_traj_file_name
528 6 : IF (filename .EQ. "") &
529 0 : CPABORT("The reference trajectory file name is empty")
530 6 : CALL parser_create(local_parser, filename, para_env=para_env)
531 : DO
532 312 : CALL parser_read_line(local_parser, 1, at_end=at_end)
533 312 : IF (at_end) EXIT
534 306 : READ (local_parser%input_line, FMT="(I8)") nread
535 306 : n_frames = n_frames + 1
536 :
537 306 : IF (n_frames > n_frames_current) THEN
538 18 : n_frames_current = 5*(n_frames_current + 10)/3
539 18 : CALL reallocate(pos_traj, 1, 3, 1, nread, 1, n_frames_current)
540 : END IF
541 :
542 : ! title line
543 306 : CALL parser_read_line(local_parser, 1)
544 :
545 : ! actual coordinates
546 2760 : DO i = 1, nread
547 2448 : CALL parser_read_line(local_parser, 1)
548 2448 : READ (local_parser%input_line(1:LEN_TRIM(local_parser%input_line)), *) AA, vec
549 10098 : pos_traj(:, i, n_frames) = vec*bohr
550 : END DO
551 :
552 : END DO
553 6 : CALL parser_release(local_parser)
554 :
555 6 : n_frames_current = n_frames
556 6 : CALL reallocate(energy_traj, 1, n_frames_current)
557 6 : CALL reallocate(force_traj, 1, 3, 1, nread, 1, n_frames_current)
558 6 : CALL reallocate(pos_traj, 1, 3, 1, nread, 1, n_frames_current)
559 :
560 : ! now force reference trajectory
561 6 : filename = oi_env%fm_env%ref_force_file_name
562 6 : IF (filename .EQ. "") &
563 0 : CPABORT("The reference force file name is empty")
564 6 : CALL parser_create(local_parser, filename, para_env=para_env)
565 312 : DO iframe = 1, n_frames
566 306 : CALL parser_read_line(local_parser, 1)
567 306 : READ (local_parser%input_line, FMT="(I8)") nread
568 :
569 : ! title line
570 306 : test_ok = .FALSE.
571 306 : CALL parser_read_line(local_parser, 1)
572 306 : READ (local_parser%input_line, FMT="(T6,I8,T23,F12.3,T41,F20.10)", ERR=999) trj_itimes, trj_time, trj_epot
573 : test_ok = .TRUE.
574 : 999 CONTINUE
575 : IF (.NOT. test_ok) THEN
576 0 : CPABORT("Could not parse the title line of the trajectory file")
577 : END IF
578 306 : energy_traj(iframe) = trj_epot
579 :
580 : ! actual forces, in a.u.
581 2760 : DO i = 1, nread
582 2448 : CALL parser_read_line(local_parser, 1)
583 2448 : READ (local_parser%input_line(1:LEN_TRIM(local_parser%input_line)), *) AA, vec
584 10098 : force_traj(:, i, iframe) = vec
585 : END DO
586 : END DO
587 6 : CALL parser_release(local_parser)
588 :
589 : ! and cell, which is optional
590 6 : NULLIFY (cell_traj)
591 6 : filename = oi_env%fm_env%ref_cell_file_name
592 6 : IF (filename .NE. "") THEN
593 6 : CALL parser_create(local_parser, filename, para_env=para_env)
594 18 : ALLOCATE (cell_traj(3, 3, n_frames))
595 312 : DO iframe = 1, n_frames
596 306 : CALL parser_read_line(local_parser, 1)
597 312 : CALL parse_cell_line(local_parser%input_line, cell_itimes, cell_time, cell_traj(:, :, iframe), vol)
598 : END DO
599 12 : CALL parser_release(local_parser)
600 : END IF
601 :
602 6 : CALL timestop(handle)
603 :
604 18 : END SUBROUTINE read_reference_data
605 :
606 : ! **************************************************************************************************
607 : !> \brief parses the input section, and stores in the optimize input environment
608 : !> \param oi_env optimize input environment
609 : !> \param root_section ...
610 : !> \author Joost VandeVondele
611 : ! **************************************************************************************************
612 12 : SUBROUTINE parse_input(oi_env, root_section)
613 : TYPE(oi_env_type) :: oi_env
614 : TYPE(section_vals_type), POINTER :: root_section
615 :
616 : CHARACTER(len=*), PARAMETER :: routineN = 'parse_input'
617 :
618 : INTEGER :: handle, ivar, n_var
619 : LOGICAL :: explicit
620 : TYPE(section_vals_type), POINTER :: fm_section, oi_section, variable_section
621 :
622 6 : CALL timeset(routineN, handle)
623 :
624 6 : CALL section_vals_val_get(root_section, "GLOBAL%PROJECT", c_val=oi_env%project_name)
625 6 : CALL section_vals_val_get(root_section, "GLOBAL%SEED", i_val=oi_env%seed)
626 : CALL cp2k_get_walltime(section=root_section, keyword_name="GLOBAL%WALLTIME", &
627 6 : walltime=oi_env%target_time)
628 :
629 6 : oi_section => section_vals_get_subs_vals(root_section, "OPTIMIZE_INPUT")
630 6 : variable_section => section_vals_get_subs_vals(oi_section, "VARIABLE")
631 :
632 6 : CALL section_vals_val_get(oi_section, "ACCURACY", r_val=oi_env%rhoend)
633 6 : CALL section_vals_val_get(oi_section, "STEP_SIZE", r_val=oi_env%rhobeg)
634 6 : CALL section_vals_val_get(oi_section, "MAX_FUN", i_val=oi_env%maxfun)
635 6 : CALL section_vals_val_get(oi_section, "ITER_START_VAL", i_val=oi_env%iter_start_val)
636 6 : CALL section_vals_val_get(oi_section, "RANDOMIZE_VARIABLES", r_val=oi_env%randomize_variables)
637 :
638 6 : CALL section_vals_get(variable_section, explicit=explicit, n_repetition=n_var)
639 6 : IF (explicit) THEN
640 30 : ALLOCATE (oi_env%variables(1:n_var))
641 18 : DO ivar = 1, n_var
642 : CALL section_vals_val_get(variable_section, "VALUE", i_rep_section=ivar, &
643 12 : r_val=oi_env%variables(ivar)%value)
644 : CALL section_vals_val_get(variable_section, "FIXED", i_rep_section=ivar, &
645 12 : l_val=oi_env%variables(ivar)%fixed)
646 : CALL section_vals_val_get(variable_section, "LABEL", i_rep_section=ivar, &
647 18 : c_val=oi_env%variables(ivar)%label)
648 : END DO
649 : END IF
650 :
651 6 : CALL section_vals_val_get(oi_section, "METHOD", i_val=oi_env%method)
652 12 : SELECT CASE (oi_env%method)
653 : CASE (opt_force_matching)
654 6 : fm_section => section_vals_get_subs_vals(oi_section, "FORCE_MATCHING")
655 6 : CALL section_vals_val_get(fm_section, "REF_TRAJ_FILE_NAME", c_val=oi_env%fm_env%ref_traj_file_name)
656 6 : CALL section_vals_val_get(fm_section, "REF_FORCE_FILE_NAME", c_val=oi_env%fm_env%ref_force_file_name)
657 6 : CALL section_vals_val_get(fm_section, "REF_CELL_FILE_NAME", c_val=oi_env%fm_env%ref_cell_file_name)
658 6 : CALL section_vals_val_get(fm_section, "OPTIMIZE_FILE_NAME", c_val=oi_env%fm_env%optimize_file_name)
659 6 : CALL section_vals_val_get(fm_section, "FRAME_START", i_val=oi_env%fm_env%frame_start)
660 6 : CALL section_vals_val_get(fm_section, "FRAME_STOP", i_val=oi_env%fm_env%frame_stop)
661 6 : CALL section_vals_val_get(fm_section, "FRAME_STRIDE", i_val=oi_env%fm_env%frame_stride)
662 6 : CALL section_vals_val_get(fm_section, "FRAME_COUNT", i_val=oi_env%fm_env%frame_count)
663 :
664 6 : CALL section_vals_val_get(fm_section, "GROUP_SIZE", i_val=oi_env%fm_env%group_size)
665 :
666 6 : CALL section_vals_val_get(fm_section, "ENERGY_WEIGHT", r_val=oi_env%fm_env%energy_weight)
667 6 : CALL section_vals_val_get(fm_section, "SHIFT_MM", r_val=oi_env%fm_env%shift_mm)
668 6 : CALL section_vals_val_get(fm_section, "SHIFT_QM", r_val=oi_env%fm_env%shift_qm)
669 6 : CALL section_vals_val_get(fm_section, "SHIFT_AVERAGE", l_val=oi_env%fm_env%shift_average)
670 : CASE DEFAULT
671 6 : CPABORT("")
672 : END SELECT
673 :
674 6 : CALL timestop(handle)
675 :
676 6 : END SUBROUTINE parse_input
677 :
678 0 : END MODULE optimize_input
|
