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