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