1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief methods to setup replicas of the same system differing only by atom 8!> positions and velocities (as used in path integral or nudged elastic 9!> band for example) 10!> \par History 11!> 09.2005 created [fawzi] 12!> \author fawzi 13! ************************************************************************************************** 14MODULE replica_methods 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_files, ONLY: close_file,& 17 open_file 18 USE cp_log_handling, ONLY: cp_get_default_logger,& 19 cp_logger_get_default_io_unit,& 20 cp_logger_type,& 21 cp_to_string 22 USE cp_output_handling, ONLY: cp_add_iter_level 23 USE cp_para_env, ONLY: cp_cart_create,& 24 cp_para_env_create 25 USE cp_para_types, ONLY: cp_para_cart_type,& 26 cp_para_env_type 27 USE cp_result_types, ONLY: cp_result_create,& 28 cp_result_retain 29 USE cp_subsys_types, ONLY: cp_subsys_get,& 30 cp_subsys_set,& 31 cp_subsys_type 32 USE f77_interface, ONLY: calc_force,& 33 create_force_env,& 34 f_env_add_defaults,& 35 f_env_rm_defaults,& 36 f_env_type,& 37 get_nparticle,& 38 get_pos,& 39 set_vel 40 USE force_env_types, ONLY: force_env_get,& 41 use_qs_force 42 USE input_section_types, ONLY: section_type,& 43 section_vals_type,& 44 section_vals_val_get,& 45 section_vals_val_set,& 46 section_vals_write 47 USE kinds, ONLY: default_path_length,& 48 dp 49 USE message_passing, ONLY: mp_cart_create,& 50 mp_cart_sub,& 51 mp_comm_null,& 52 mp_sum 53 USE qs_environment_types, ONLY: get_qs_env,& 54 qs_environment_type,& 55 set_qs_env 56 USE qs_wf_history_methods, ONLY: wfi_create,& 57 wfi_create_for_kp 58 USE qs_wf_history_types, ONLY: wfi_retain 59 USE replica_types, ONLY: rep_env_sync,& 60 rep_env_sync_results,& 61 rep_envs_add_rep_env,& 62 rep_envs_get_rep_env,& 63 replica_env_type 64#include "./base/base_uses.f90" 65 66 IMPLICIT NONE 67 PRIVATE 68 69 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 70 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'replica_methods' 71 INTEGER, SAVE, PRIVATE :: last_rep_env_id = 0 72 73 PUBLIC :: rep_env_create, rep_env_calc_e_f 74 75CONTAINS 76 77! ************************************************************************************************** 78!> \brief creates a replica environment together with its force environment 79!> \param rep_env the replica environment that will be created 80!> \param para_env the parallel environment that will contain the replicas 81!> \param input the input used to initialize the force environment 82!> \param input_declaration ... 83!> \param nrep the number of replicas to calculate 84!> \param prep the number of processors for each replica 85!> \param sync_v if the velocity should be synchronized (defaults to false) 86!> \param keep_wf_history if wf history should be kept on a per replica 87!> basis (defaults to true for QS jobs) 88!> \param row_force to use the new mapping to the cart with rows 89!> working on force instead of columns. 90!> \author fawzi 91! ************************************************************************************************** 92 SUBROUTINE rep_env_create(rep_env, para_env, input, input_declaration, nrep, prep, & 93 sync_v, keep_wf_history, row_force) 94 TYPE(replica_env_type), POINTER :: rep_env 95 TYPE(cp_para_env_type), POINTER :: para_env 96 TYPE(section_vals_type), POINTER :: input 97 TYPE(section_type), POINTER :: input_declaration 98 INTEGER :: nrep, prep 99 LOGICAL, INTENT(in), OPTIONAL :: sync_v, keep_wf_history, row_force 100 101 CHARACTER(len=*), PARAMETER :: routineN = 'rep_env_create', routineP = moduleN//':'//routineN 102 103 CHARACTER(len=default_path_length) :: input_file_path, output_file_path 104 INTEGER :: comm_cart, comm_f, comm_inter_rep, forcedim, i, i0, ierr, ip, ir, irep, lp, & 105 my_prep, new_env_id, nparticle, nrep_local, unit_nr 106 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: gridinfo 107 INTEGER, DIMENSION(2) :: dims, pos 108 LOGICAL, DIMENSION(2) :: rdim 109 TYPE(cp_logger_type), POINTER :: logger 110 TYPE(cp_para_cart_type), POINTER :: cart 111 TYPE(cp_para_env_type), POINTER :: para_env_f, para_env_full, & 112 para_env_inter_rep 113 114 CPASSERT(.NOT. ASSOCIATED(rep_env)) 115 CPASSERT(ASSOCIATED(input_declaration)) 116 117 NULLIFY (cart, para_env_f, para_env_inter_rep) 118 logger => cp_get_default_logger() 119 unit_nr = cp_logger_get_default_io_unit(logger) 120 new_env_id = -1 121 forcedim = 1 122 IF (PRESENT(row_force)) THEN 123 IF (row_force) forcedim = 2 124 END IF 125 my_prep = MIN(prep, para_env%num_pe) 126 dims(3 - forcedim) = MIN(para_env%num_pe/my_prep, nrep) 127 dims(forcedim) = my_prep 128 IF ((dims(1)*dims(2) /= para_env%num_pe) .AND. (unit_nr > 0)) THEN 129 WRITE (unit_nr, FMT="(T2,A)") "REPLICA| WARNING: number of processors is not divisible by the number of replicas" 130 WRITE (unit_nr, FMT="(T2,A,I0,A)") "REPLICA| ", para_env%num_pe - dims(1)*dims(2), " MPI process(es) will be idle" 131 END IF 132 CALL mp_cart_create(comm_old=para_env%group, ndims=2, dims=dims, pos=pos, comm_cart=comm_cart) 133 IF (comm_cart /= mp_comm_null) THEN 134 CALL cp_cart_create(cart, comm_cart, ndims=2, owns_group=.TRUE.) 135 NULLIFY (para_env_full) 136 CALL cp_para_env_create(para_env_full, comm_cart, owns_group=.FALSE.) 137 rdim(3 - forcedim) = .FALSE. 138 rdim(forcedim) = .TRUE. 139 CALL mp_cart_sub(comm=comm_cart, rdim=rdim, sub_comm=comm_f) 140 CALL cp_para_env_create(para_env_f, comm_f, owns_group=.TRUE.) 141 rdim(3 - forcedim) = .TRUE. 142 rdim(forcedim) = .FALSE. 143 CALL mp_cart_sub(comm=comm_cart, rdim=rdim, sub_comm=comm_inter_rep) 144 CALL cp_para_env_create(para_env_inter_rep, comm_inter_rep, & 145 owns_group=.TRUE.) 146 ALLOCATE (rep_env) 147 END IF 148 ALLOCATE (gridinfo(2, 0:para_env%num_pe - 1)) 149 gridinfo = 0 150 gridinfo(:, para_env%mepos) = pos 151 CALL mp_sum(gridinfo, para_env%group) 152 IF (unit_nr > 0) THEN 153 WRITE (unit_nr, FMT="(T2,A,T71,I10)") "REPLICA| layout of the replica grid, number of groups ", para_env_inter_rep%num_pe 154 WRITE (unit_nr, FMT="(T2,A,T71,I10)") "REPLICA| layout of the replica grid, size of each group", para_env_f%num_pe 155 WRITE (unit_nr, FMT="(T2,A)", ADVANCE="NO") "REPLICA| MPI process to grid (group,rank) correspondence:" 156 DO i = 0, para_env%num_pe - 1 157 IF (MODULO(i, 4) == 0) WRITE (unit_nr, *) 158 WRITE (unit_nr, FMT='(A3,I4,A3,I4,A1,I4,A1)', ADVANCE="NO") & 159 " (", i, " : ", gridinfo(3 - forcedim, i), ",", & 160 gridinfo(forcedim, i), ")" 161 END DO 162 WRITE (unit_nr, *) 163 ENDIF 164 DEALLOCATE (gridinfo) 165 IF (ASSOCIATED(rep_env)) THEN 166 last_rep_env_id = last_rep_env_id + 1 167 rep_env%id_nr = last_rep_env_id 168 rep_env%ref_count = 1 169 rep_env%nrep = nrep 170 rep_env%sync_v = .FALSE. 171 IF (PRESENT(sync_v)) rep_env%sync_v = sync_v 172 rep_env%keep_wf_history = .TRUE. 173 IF (PRESENT(keep_wf_history)) rep_env%keep_wf_history = keep_wf_history 174 NULLIFY (rep_env%wf_history) 175 NULLIFY (rep_env%results) 176 177 rep_env%force_dim = forcedim 178 rep_env%my_rep_group = cart%mepos(3 - forcedim) 179 ALLOCATE (rep_env%inter_rep_rank(0:para_env_inter_rep%num_pe - 1), & 180 rep_env%force_rank(0:para_env_f%num_pe - 1)) 181 rep_env%inter_rep_rank = 0 182 rep_env%inter_rep_rank(rep_env%my_rep_group) = para_env_inter_rep%mepos 183 CALL mp_sum(rep_env%inter_rep_rank, para_env_inter_rep%group) 184 rep_env%force_rank = 0 185 rep_env%force_rank(cart%mepos(forcedim)) = para_env_f%mepos 186 CALL mp_sum(rep_env%force_rank, para_env_f%group) 187 188 CALL section_vals_val_get(input, "GLOBAL%PROJECT_NAME", & 189 c_val=input_file_path) 190 rep_env%original_project_name = input_file_path 191 ! By default replica_env handles files for each replica 192 ! with the structure PROJECT_NAME-r-N where N is the 193 ! number of the local replica.. 194 lp = LEN_TRIM(input_file_path) 195 input_file_path(lp + 1:LEN(input_file_path)) = "-r-"// & 196 ADJUSTL(cp_to_string(rep_env%my_rep_group)) 197 lp = LEN_TRIM(input_file_path) 198 ! Setup new project name 199 CALL section_vals_val_set(input, "GLOBAL%PROJECT_NAME", & 200 c_val=input_file_path) 201 ! Redirect the output of each replica on a same local file 202 output_file_path = input_file_path(1:lp)//".out" 203 CALL section_vals_val_set(input, "GLOBAL%OUTPUT_FILE_NAME", & 204 c_val=TRIM(output_file_path)) 205 206 ! Dump an input file to warm-up new force_eval structures and 207 ! delete them immediately afterwards.. 208 input_file_path(lp + 1:LEN(input_file_path)) = ".inp" 209 IF (para_env_f%ionode) THEN 210 CALL open_file(file_name=TRIM(input_file_path), file_status="UNKNOWN", & 211 file_form="FORMATTED", file_action="WRITE", & 212 unit_number=unit_nr) 213 CALL section_vals_write(input, unit_nr, hide_root=.TRUE.) 214 CALL close_file(unit_nr) 215 END IF 216 CALL create_force_env(new_env_id, input_declaration, input_file_path, & 217 output_file_path, para_env_f%group, ierr=ierr) 218 CPASSERT(ierr == 0) 219 220 ! Delete input files.. 221 IF (para_env_f%ionode) THEN 222 CALL open_file(file_name=TRIM(input_file_path), file_status="OLD", & 223 file_form="FORMATTED", file_action="READ", unit_number=unit_nr) 224 CALL close_file(unit_number=unit_nr, file_status="DELETE") 225 END IF 226 227 rep_env%f_env_id = new_env_id 228 CALL get_nparticle(new_env_id, nparticle, ierr) 229 CPASSERT(ierr == 0) 230 rep_env%nparticle = nparticle 231 rep_env%ndim = 3*nparticle 232 ALLOCATE (rep_env%replica_owner(nrep)) 233 234 i0 = nrep/para_env_inter_rep%num_pe 235 ir = MODULO(nrep, para_env_inter_rep%num_pe) 236 DO ip = 0, para_env_inter_rep%num_pe - 1 237 DO i = i0*ip + MIN(ip, ir) + 1, i0*(ip + 1) + MIN(ip + 1, ir) 238 rep_env%replica_owner(i) = ip 239 END DO 240 END DO 241 242 nrep_local = i0 243 IF (rep_env%my_rep_group < ir) nrep_local = nrep_local + 1 244 ALLOCATE (rep_env%local_rep_indices(nrep_local), & 245 rep_env%rep_is_local(nrep)) 246 nrep_local = 0 247 rep_env%rep_is_local = .FALSE. 248 DO irep = 1, nrep 249 IF (rep_env%replica_owner(irep) == rep_env%my_rep_group) THEN 250 nrep_local = nrep_local + 1 251 rep_env%local_rep_indices(nrep_local) = irep 252 rep_env%rep_is_local(irep) = .TRUE. 253 END IF 254 END DO 255 CPASSERT(nrep_local == SIZE(rep_env%local_rep_indices)) 256 257 rep_env%cart => cart 258 rep_env%para_env => para_env_full 259 rep_env%para_env_f => para_env_f 260 rep_env%para_env_inter_rep => para_env_inter_rep 261 262 ALLOCATE (rep_env%r(rep_env%ndim, nrep), rep_env%v(rep_env%ndim, nrep), & 263 rep_env%f(rep_env%ndim + 1, nrep)) 264 265 rep_env%r = 0._dp 266 rep_env%f = 0._dp 267 rep_env%v = 0._dp 268 CALL set_vel(rep_env%f_env_id, rep_env%v(:, 1), rep_env%ndim, ierr) 269 CPASSERT(ierr == 0) 270 DO i = 1, nrep 271 IF (rep_env%rep_is_local(i)) THEN 272 CALL get_pos(rep_env%f_env_id, rep_env%r(:, i), rep_env%ndim, ierr) 273 CPASSERT(ierr == 0) 274 END IF 275 END DO 276 END IF 277 IF (ASSOCIATED(rep_env)) THEN 278 CALL rep_envs_add_rep_env(rep_env) 279 CALL rep_env_init_low(rep_env%id_nr, ierr) 280 CPASSERT(ierr == 0) 281 END IF 282 END SUBROUTINE rep_env_create 283 284! ************************************************************************************************** 285!> \brief finishes the low level initialization of the replica env 286!> \param rep_env_id id_nr of the replica environment that should be initialized 287!> \param ierr will be non zero if there is an initialization error 288!> \author fawzi 289! ************************************************************************************************** 290 SUBROUTINE rep_env_init_low(rep_env_id, ierr) 291 INTEGER, INTENT(in) :: rep_env_id 292 INTEGER, INTENT(out) :: ierr 293 294 CHARACTER(len=*), PARAMETER :: routineN = 'rep_env_init_low', & 295 routineP = moduleN//':'//routineN 296 297 INTEGER :: i, in_use, stat 298 LOGICAL :: do_kpoints, has_unit_metric 299 TYPE(cp_logger_type), POINTER :: logger 300 TYPE(cp_subsys_type), POINTER :: subsys 301 TYPE(dft_control_type), POINTER :: dft_control 302 TYPE(f_env_type), POINTER :: f_env 303 TYPE(qs_environment_type), POINTER :: qs_env 304 TYPE(replica_env_type), POINTER :: rep_env 305 306 rep_env => rep_envs_get_rep_env(rep_env_id, ierr=stat) 307 IF (.NOT. ASSOCIATED(rep_env)) & 308 CPABORT("could not find rep_env with id_nr"//cp_to_string(rep_env_id)) 309 NULLIFY (qs_env, dft_control, subsys) 310 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env) 311 logger => cp_get_default_logger() 312 logger%iter_info%iteration(1) = rep_env%my_rep_group 313 CALL cp_add_iter_level(iteration_info=logger%iter_info, & 314 level_name="REPLICA_EVAL") 315 !wf interp 316 IF (rep_env%keep_wf_history) THEN 317 CALL force_env_get(f_env%force_env, in_use=in_use) 318 IF (in_use == use_qs_force) THEN 319 CALL force_env_get(f_env%force_env, qs_env=qs_env) 320 CALL get_qs_env(qs_env, dft_control=dft_control) 321 ALLOCATE (rep_env%wf_history(SIZE(rep_env%local_rep_indices))) 322 DO i = 1, SIZE(rep_env%wf_history) 323 NULLIFY (rep_env%wf_history(i)%wf_history) 324 IF (i == 1) THEN 325 CALL get_qs_env(qs_env, & 326 wf_history=rep_env%wf_history(i)%wf_history) 327 CALL wfi_retain(rep_env%wf_history(i)%wf_history) 328 ELSE 329 CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric, & 330 do_kpoints=do_kpoints) 331 CALL wfi_create(rep_env%wf_history(i)%wf_history, & 332 interpolation_method_nr= & 333 dft_control%qs_control%wf_interpolation_method_nr, & 334 extrapolation_order=dft_control%qs_control%wf_extrapolation_order, & 335 has_unit_metric=has_unit_metric) 336 IF (do_kpoints) THEN 337 CALL wfi_create_for_kp(rep_env%wf_history(i)%wf_history) 338 END IF 339 END IF 340 END DO 341 ELSE 342 rep_env%keep_wf_history = .FALSE. 343 END IF 344 END IF 345 ALLOCATE (rep_env%results(rep_env%nrep)) 346 DO i = 1, rep_env%nrep 347 NULLIFY (rep_env%results(i)%results) 348 IF (i == 1) THEN 349 CALL force_env_get(f_env%force_env, subsys=subsys) 350 CALL cp_subsys_get(subsys, results=rep_env%results(i)%results) 351 CALL cp_result_retain(rep_env%results(i)%results) 352 ELSE 353 CALL cp_result_create(rep_env%results(i)%results) 354 END IF 355 END DO 356 CALL rep_env_sync(rep_env, rep_env%r) 357 CALL rep_env_sync(rep_env, rep_env%v) 358 CALL rep_env_sync(rep_env, rep_env%f) 359 360 CALL f_env_rm_defaults(f_env, ierr) 361 CPASSERT(ierr == 0) 362 END SUBROUTINE rep_env_init_low 363 364! ************************************************************************************************** 365!> \brief evaluates the forces 366!> \param rep_env the replica environment on which you want to evaluate the 367!> forces 368!> \param calc_f if true calculates also the forces, if false only the 369!> energy 370!> \author fawzi 371!> \note 372!> indirect through f77_int_low to work around fortran madness 373! ************************************************************************************************** 374 SUBROUTINE rep_env_calc_e_f(rep_env, calc_f) 375 TYPE(replica_env_type), POINTER :: rep_env 376 LOGICAL, OPTIONAL :: calc_f 377 378 CHARACTER(len=*), PARAMETER :: routineN = 'rep_env_calc_e_f', & 379 routineP = moduleN//':'//routineN 380 381 INTEGER :: handle, ierr, my_calc_f 382 383 CALL timeset(routineN, handle) 384 CPASSERT(ASSOCIATED(rep_env)) 385 CPASSERT(rep_env%ref_count > 0) 386 my_calc_f = 0 387 IF (PRESENT(calc_f)) THEN 388 IF (calc_f) my_calc_f = 1 389 END IF 390 CALL rep_env_calc_e_f_low(rep_env%id_nr, my_calc_f, ierr) 391 CPASSERT(ierr == 0) 392 CALL timestop(handle) 393 END SUBROUTINE rep_env_calc_e_f 394 395! ************************************************************************************************** 396!> \brief calculates energy and force, internal private method 397!> \param rep_env_id the id if the replica environment in which energy and 398!> forces have to be evaluated 399!> \param calc_f if nonzero calculates also the forces along with the 400!> energy 401!> \param ierr if an error happens this will be nonzero 402!> \author fawzi 403!> \note 404!> low level wrapper to export this function in f77_int_low and work 405!> around the handling of circular dependencies in fortran 406! ************************************************************************************************** 407 RECURSIVE SUBROUTINE rep_env_calc_e_f_low(rep_env_id, calc_f, ierr) 408 INTEGER, INTENT(in) :: rep_env_id, calc_f 409 INTEGER, INTENT(out) :: ierr 410 411 CHARACTER(len=*), PARAMETER :: routineN = 'rep_env_calc_e_f_low', & 412 routineP = moduleN//':'//routineN 413 414 TYPE(f_env_type), POINTER :: f_env 415 TYPE(replica_env_type), POINTER :: rep_env 416 417 rep_env => rep_envs_get_rep_env(rep_env_id, ierr) 418 IF (ASSOCIATED(rep_env)) THEN 419 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env) 420 CALL rep_env_calc_e_f_int(rep_env, calc_f /= 0) 421 CALL f_env_rm_defaults(f_env, ierr) 422 ELSE 423 ierr = 111 424 END IF 425 END SUBROUTINE rep_env_calc_e_f_low 426 427! ************************************************************************************************** 428!> \brief calculates energy and force, internal private method 429!> \param rep_env the replica env to update 430!> \param calc_f if the force should be calculated as well (defaults to true) 431!> \author fawzi 432!> \note 433!> this is the where the real work is done 434! ************************************************************************************************** 435 SUBROUTINE rep_env_calc_e_f_int(rep_env, calc_f) 436 TYPE(replica_env_type), POINTER :: rep_env 437 LOGICAL, OPTIONAL :: calc_f 438 439 CHARACTER(len=*), PARAMETER :: routineN = 'rep_env_calc_e_f_int', & 440 routineP = moduleN//':'//routineN 441 442 INTEGER :: i, ierr, irep, md_iter, my_calc_f, ndim 443 TYPE(cp_logger_type), POINTER :: logger 444 TYPE(cp_subsys_type), POINTER :: subsys 445 TYPE(f_env_type), POINTER :: f_env 446 TYPE(qs_environment_type), POINTER :: qs_env 447 448 NULLIFY (f_env, qs_env, subsys) 449 CPASSERT(ASSOCIATED(rep_env)) 450 CPASSERT(rep_env%ref_count > 0) 451 my_calc_f = 3*rep_env%nparticle 452 IF (PRESENT(calc_f)) THEN 453 IF (.NOT. calc_f) my_calc_f = 0 454 END IF 455 456 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env) 457 logger => cp_get_default_logger() 458 ! md_iter=logger%iter_info%iteration(2)+1 459 md_iter = logger%iter_info%iteration(2) 460 CALL f_env_rm_defaults(f_env, ierr) 461 CPASSERT(ierr == 0) 462 DO i = 1, SIZE(rep_env%local_rep_indices) 463 irep = rep_env%local_rep_indices(i) 464 ndim = 3*rep_env%nparticle 465 IF (rep_env%sync_v) THEN 466 CALL set_vel(rep_env%f_env_id, rep_env%v(:, irep), ndim, ierr) 467 CPASSERT(ierr == 0) 468 END IF 469 470 logger%iter_info%iteration(1) = irep 471 logger%iter_info%iteration(2) = md_iter 472 473 IF (rep_env%keep_wf_history) THEN 474 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env) 475 CALL force_env_get(f_env%force_env, qs_env=qs_env) 476 CALL set_qs_env(qs_env, & 477 wf_history=rep_env%wf_history(i)%wf_history) 478 CALL f_env_rm_defaults(f_env, ierr) 479 CPASSERT(ierr == 0) 480 END IF 481 482 CALL f_env_add_defaults(f_env_id=rep_env%f_env_id, f_env=f_env) 483 CALL force_env_get(f_env%force_env, subsys=subsys) 484 CALL cp_subsys_set(subsys, results=rep_env%results(irep)%results) 485 CALL f_env_rm_defaults(f_env, ierr) 486 CPASSERT(ierr == 0) 487 CALL calc_force(rep_env%f_env_id, rep_env%r(:, irep), ndim, & 488 rep_env%f(ndim + 1, irep), rep_env%f(:ndim, irep), & 489 my_calc_f, ierr) 490 CPASSERT(ierr == 0) 491 END DO 492 CALL rep_env_sync(rep_env, rep_env%f) 493 CALL rep_env_sync_results(rep_env, rep_env%results) 494 495 END SUBROUTINE rep_env_calc_e_f_int 496 497END MODULE replica_methods 498