1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \par History 8!> \author 9! ************************************************************************************************** 10MODULE gle_system_dynamics 11 12 USE cp_para_types, ONLY: cp_para_env_type 13 USE distribution_1d_types, ONLY: distribution_1d_type 14 USE extended_system_types, ONLY: map_info_type 15 USE gle_system_types, ONLY: gle_thermo_create,& 16 gle_type 17 USE input_constants, ONLY: & 18 do_thermo_communication, do_thermo_no_communication, isokin_ensemble, langevin_ensemble, & 19 npe_f_ensemble, npe_i_ensemble, nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, & 20 npt_f_ensemble, npt_i_ensemble, nve_ensemble, nvt_ensemble, reftraj_ensemble 21 USE input_section_types, ONLY: section_vals_get,& 22 section_vals_get_subs_vals,& 23 section_vals_remove_values,& 24 section_vals_type,& 25 section_vals_val_get 26 USE kinds, ONLY: dp 27 USE message_passing, ONLY: mp_sum 28 USE molecule_kind_types, ONLY: molecule_kind_type 29 USE molecule_types, ONLY: global_constraint_type,& 30 molecule_type 31 USE parallel_rng_types, ONLY: next_random_number,& 32 read_rng_stream,& 33 rng_record_length,& 34 rng_stream_type 35 USE particle_types, ONLY: particle_type 36 USE simpar_types, ONLY: simpar_type 37 USE thermostat_mapping, ONLY: thermostat_mapping_region 38 USE thermostat_types, ONLY: thermostat_info_type 39 USE thermostat_utils, ONLY: ke_region_particles,& 40 momentum_region_particles,& 41 vel_rescale_particles 42#include "../../base/base_uses.f90" 43 44 IMPLICIT NONE 45 PRIVATE 46 47 PUBLIC :: gle_particles, & 48 initialize_gle_part, & 49 gle_cholesky_stab, & 50 gle_matrix_exp, & 51 restart_gle 52 53 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gle_system_dynamics' 54 55CONTAINS 56 57! ************************************************************************************************** 58!> \brief ... 59!> \param gle ... 60!> \param molecule_kind_set ... 61!> \param molecule_set ... 62!> \param particle_set ... 63!> \param local_molecules ... 64!> \param group ... 65!> \param shell_adiabatic ... 66!> \param shell_particle_set ... 67!> \param core_particle_set ... 68!> \param vel ... 69!> \param shell_vel ... 70!> \param core_vel ... 71!> \date 72!> \par History 73! ************************************************************************************************** 74 SUBROUTINE gle_particles(gle, molecule_kind_set, molecule_set, particle_set, local_molecules, & 75 group, shell_adiabatic, shell_particle_set, core_particle_set, vel, shell_vel, core_vel) 76 77 TYPE(gle_type), POINTER :: gle 78 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 79 TYPE(molecule_type), POINTER :: molecule_set(:) 80 TYPE(particle_type), POINTER :: particle_set(:) 81 TYPE(distribution_1d_type), POINTER :: local_molecules 82 INTEGER, INTENT(IN) :: group 83 LOGICAL, INTENT(IN), OPTIONAL :: shell_adiabatic 84 TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), & 85 core_particle_set(:) 86 REAL(KIND=dp), INTENT(INOUT), OPTIONAL :: vel(:, :), shell_vel(:, :), & 87 core_vel(:, :) 88 89 CHARACTER(len=*), PARAMETER :: routineN = 'gle_particles', routineP = moduleN//':'//routineN 90 91 INTEGER :: handle, iadd, ideg, imap, ndim, num 92 LOGICAL :: my_shell_adiabatic, present_vel 93 REAL(dp) :: alpha, beta, rr 94 REAL(dp), DIMENSION(:, :), POINTER :: a_mat, e_tmp, h_tmp, s_tmp 95 TYPE(map_info_type), POINTER :: map_info 96 TYPE(rng_stream_type), POINTER :: rng_stream 97 98 CALL timeset(routineN, handle) 99 my_shell_adiabatic = .FALSE. 100 IF (PRESENT(shell_adiabatic)) my_shell_adiabatic = shell_adiabatic 101 NULLIFY (rng_stream) 102 present_vel = PRESENT(vel) 103 ndim = gle%ndim 104 ALLOCATE (s_tmp(ndim, gle%loc_num_gle)) 105 s_tmp = 0.0_dp 106 ALLOCATE (e_tmp(ndim, gle%loc_num_gle)) 107 ALLOCATE (h_tmp(ndim, gle%loc_num_gle)) 108 109 map_info => gle%map_info 110 CALL ke_region_particles(map_info, particle_set, molecule_kind_set, & 111 local_molecules, molecule_set, group, vel) 112 DO ideg = 1, gle%loc_num_gle 113 imap = gle%map_info%map_index(ideg) 114 gle%nvt(ideg)%kin_energy = map_info%s_kin(imap) 115 END DO 116 117 CALL momentum_region_particles(map_info, particle_set, molecule_kind_set, & 118 local_molecules, molecule_set, group, vel) 119 120 DO ideg = 1, gle%loc_num_gle 121 imap = gle%map_info%map_index(ideg) 122 IF (gle%nvt(ideg)%nkt == 0.0_dp) CYCLE 123 gle%nvt(ideg)%s(1) = map_info%s_kin(imap) 124 s_tmp(1, imap) = map_info%s_kin(imap) 125 rng_stream => gle%nvt(ideg)%gaussian_rng_stream 126 rr = next_random_number(rng_stream) 127 e_tmp(1, imap) = rr 128 DO iadd = 2, ndim 129 s_tmp(iadd, imap) = gle%nvt(ideg)%s(iadd) 130 rr = next_random_number(rng_stream) 131 e_tmp(iadd, imap) = rr 132 END DO 133 END DO 134 num = gle%loc_num_gle 135 a_mat => gle%gle_s 136 alpha = 1.0_dp 137 beta = 0.0_dp 138 CALL DGEMM('N', 'N', ndim, num, ndim, alpha, a_mat(1, 1), ndim, e_tmp(1, 1), ndim, beta, h_tmp(1, 1), ndim) 139! 140 a_mat => gle%gle_t 141 beta = 1.0_dp 142 CALL dgemm("N", "N", ndim, num, ndim, alpha, a_mat(1, 1), ndim, s_tmp(1, 1), ndim, beta, h_tmp(1, 1), ndim) 143 144 DO ideg = 1, gle%loc_num_gle 145 imap = gle%map_info%map_index(ideg) 146 IF (gle%nvt(ideg)%nkt == 0.0_dp) CYCLE 147 148 map_info%v_scale(imap) = h_tmp(1, imap)/s_tmp(1, imap) 149 DO iadd = 2, ndim 150 gle%nvt(ideg)%s(iadd) = h_tmp(iadd, ideg) 151 END DO 152 END DO 153 154 CALL vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, & 155 local_molecules, my_shell_adiabatic, shell_particle_set, core_particle_set, & 156 vel, shell_vel, core_vel) 157 158 CALL ke_region_particles(map_info, particle_set, molecule_kind_set, & 159 local_molecules, molecule_set, group, vel) 160 DO ideg = 1, gle%loc_num_gle 161 imap = gle%map_info%map_index(ideg) 162 gle%nvt(ideg)%thermostat_energy = gle%nvt(ideg)%thermostat_energy + & 163 0.5_dp*(gle%nvt(ideg)%kin_energy - map_info%s_kin(imap)) 164 END DO 165 166 DEALLOCATE (e_tmp, s_tmp, h_tmp) 167 168 CALL timestop(handle) 169 END SUBROUTINE gle_particles 170 171! ************************************************************************************************** 172!> \brief ... 173!> \param thermostat_info ... 174!> \param simpar ... 175!> \param local_molecules ... 176!> \param molecule ... 177!> \param molecule_kind_set ... 178!> \param para_env ... 179!> \param gle ... 180!> \param gle_section ... 181!> \param gci ... 182!> \param save_mem ... 183!> \author 184! ************************************************************************************************** 185 SUBROUTINE initialize_gle_part(thermostat_info, simpar, local_molecules, & 186 molecule, molecule_kind_set, para_env, gle, gle_section, & 187 gci, save_mem) 188 189 TYPE(thermostat_info_type), POINTER :: thermostat_info 190 TYPE(simpar_type), POINTER :: simpar 191 TYPE(distribution_1d_type), POINTER :: local_molecules 192 TYPE(molecule_type), POINTER :: molecule(:) 193 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 194 TYPE(cp_para_env_type), POINTER :: para_env 195 TYPE(gle_type), POINTER :: gle 196 TYPE(section_vals_type), POINTER :: gle_section 197 TYPE(global_constraint_type), POINTER :: gci 198 LOGICAL, INTENT(IN) :: save_mem 199 200 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_gle_part', & 201 routineP = moduleN//':'//routineN 202 203 LOGICAL :: restart 204 REAL(dp) :: Mtmp(gle%ndim, gle%ndim) 205 206 restart = .FALSE. 207 208 CALL gle_to_particle_mapping(thermostat_info, simpar, local_molecules, & 209 molecule, molecule_kind_set, gle, para_env, gci) 210 211 IF (gle%ndim /= 0) THEN 212 CALL init_gle_variables(gle) 213 END IF 214 CALL restart_gle(gle, gle_section, save_mem, restart) 215 216 ! here we should have read a_mat and c_mat; whe can therefore compute S and T 217 ! deterministic part of the propagator 218 CALL gle_matrix_exp((-simpar%dt*0.5_dp)*gle%a_mat, gle%ndim, 15, 15, gle%gle_t) 219 ! stochastic part 220 Mtmp = gle%c_mat - MATMUL(gle%gle_t, MATMUL(gle%c_mat, TRANSPOSE(gle%gle_t))) 221 CALL gle_cholesky_stab(Mtmp, gle%gle_s, gle%ndim) 222 223 END SUBROUTINE initialize_gle_part 224 225! ************************************************************************************************** 226!> \brief ... 227!> \param M ... 228!> \param n ... 229!> \param j ... 230!> \param k ... 231!> \param EM ... 232!> \author Michele 233!> routine to compute matrix exponential via scale & square 234! ************************************************************************************************** 235 SUBROUTINE gle_matrix_exp(M, n, j, k, EM) 236 237 INTEGER, INTENT(in) :: n 238 REAL(dp), INTENT(in) :: M(n, n) 239 INTEGER, INTENT(in) :: j, k 240 REAL(dp), INTENT(out) :: EM(n, n) 241 242 INTEGER :: i, p 243 REAL(dp) :: SM(n, n), tc(j + 1) 244 245 tc(1) = 1._dp 246 DO i = 1, j 247 tc(i + 1) = tc(i)/REAL(i, KIND=dp) 248 ENDDO 249 250 !scale 251 SM = M*(1._dp/2._dp**k) 252 EM = 0._dp 253 DO i = 1, n 254 EM(i, i) = tc(j + 1) 255 ENDDO 256 257 !taylor exp of scaled matrix 258 DO p = j, 1, -1 259 EM = MATMUL(SM, EM) 260 DO i = 1, n 261 EM(i, i) = EM(i, i) + tc(p) 262 ENDDO 263 ENDDO 264 265 !square 266 DO p = 1, k 267 EM = MATMUL(EM, EM) 268 ENDDO 269 END SUBROUTINE gle_matrix_exp 270 271! ************************************************************************************************** 272!> \brief ... 273!> \param SST ... 274!> \param S ... 275!> \param n ... 276!> \author Michele 277!> "stabilized" cholesky to deal with small & negative diagonal elements, 278!> in practice a LDL^T decomposition is computed, and negative els. are zeroed 279!> \par History 280!> 05.11.2015: Bug fix: In rare cases D(j) is negative due to numerics [Felix Uhl] 281! ************************************************************************************************** 282 SUBROUTINE gle_cholesky_stab(SST, S, n) 283 INTEGER, INTENT(in) :: n 284 REAL(dp), INTENT(out) :: S(n, n) 285 REAL(dp), INTENT(in) :: SST(n, n) 286 287 INTEGER :: i, j, k 288 REAL(dp) :: D(n), L(n, n) 289 290 D = 0._dp 291 L = 0._dp 292 S = 0._dp 293 DO i = 1, n 294 L(i, i) = 1.0_dp 295 D(i) = SST(i, i) 296 DO j = 1, i - 1 297 L(i, j) = SST(i, j); 298 DO k = 1, j - 1 299 L(i, j) = L(i, j) - L(i, k)*L(j, k)*D(k) 300 ENDDO 301 IF (ABS(D(j)) > EPSILON(1.0_dp)) L(i, j) = L(i, j)/D(j) 302 ENDDO 303 DO k = 1, i - 1 304 D(i) = D(i) - L(i, k)*L(i, k)*D(k) 305 END DO 306 ENDDO 307 DO i = 1, n 308 DO j = 1, i 309 IF ((ABS(D(j)) > EPSILON(1.0_dp)) .AND. (D(j) > 0.0_dp)) THEN 310 S(i, j) = S(i, j) + L(i, j)*SQRT(D(j)) 311 END IF 312 ENDDO 313 ENDDO 314 315 END SUBROUTINE gle_cholesky_stab 316 317! ************************************************************************************************** 318!> \brief ... 319!> \param thermostat_info ... 320!> \param simpar ... 321!> \param local_molecules ... 322!> \param molecule_set ... 323!> \param molecule_kind_set ... 324!> \param gle ... 325!> \param para_env ... 326!> \param gci ... 327!> \author 328! ************************************************************************************************** 329 SUBROUTINE gle_to_particle_mapping(thermostat_info, simpar, local_molecules, & 330 molecule_set, molecule_kind_set, gle, para_env, gci) 331 332 TYPE(thermostat_info_type), POINTER :: thermostat_info 333 TYPE(simpar_type), POINTER :: simpar 334 TYPE(distribution_1d_type), POINTER :: local_molecules 335 TYPE(molecule_type), POINTER :: molecule_set(:) 336 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 337 TYPE(gle_type), POINTER :: gle 338 TYPE(cp_para_env_type), POINTER :: para_env 339 TYPE(global_constraint_type), POINTER :: gci 340 341 CHARACTER(LEN=*), PARAMETER :: routineN = 'gle_to_particle_mapping', & 342 routineP = moduleN//':'//routineN 343 344 INTEGER :: i, imap, j, mal_size, natoms_local, & 345 nkind, number, region, & 346 sum_of_thermostats 347 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list 348 LOGICAL :: do_shell 349 REAL(KIND=dp) :: fac 350 TYPE(map_info_type), POINTER :: map_info 351 352 do_shell = .FALSE. 353 NULLIFY (massive_atom_list, deg_of_freedom) 354 SELECT CASE (simpar%ensemble) 355 CASE DEFAULT 356 CPABORT("Unknown ensemble!") 357 CASE (nve_ensemble, isokin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_ensemble, & 358 nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble) 359 CPABORT("Never reach this point!") 360 CASE (nvt_ensemble, npt_i_ensemble, npt_f_ensemble) 361 362 map_info => gle%map_info 363 nkind = SIZE(molecule_kind_set) 364 sum_of_thermostats = thermostat_info%sum_of_thermostats 365 map_info%dis_type = thermostat_info%dis_type 366 number = thermostat_info%number_of_thermostats 367 region = gle%region 368 369 CALL thermostat_mapping_region(map_info, deg_of_freedom, massive_atom_list, & 370 molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, & 371 simpar, number, region, gci, do_shell, thermostat_info%map_loc_thermo_gen, & 372 sum_of_thermostats) 373 374 ! This is the local number of available thermostats 375 gle%loc_num_gle = number 376 gle%glob_num_gle = sum_of_thermostats 377 mal_size = SIZE(massive_atom_list) 378 CPASSERT(mal_size /= 0) 379 CALL gle_thermo_create(gle, mal_size) 380 gle%mal(1:mal_size) = massive_atom_list(1:mal_size) 381 382 ! Sum up the number of degrees of freedom on each thermostat. 383 ! first: initialize the target 384 map_info%s_kin = 0.0_dp 385 DO i = 1, 3 386 DO j = 1, natoms_local 387 map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1 388 END DO 389 END DO 390 391 ! If thermostats are replicated but molecules distributed, we have to 392 ! sum s_kin over all processors 393 IF (map_info%dis_type == do_thermo_communication) CALL mp_sum(map_info%s_kin, para_env%group) 394 395 ! We know the total number of system thermostats. 396 IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN 397 fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl 398 IF (fac == 0.0_dp) THEN 399 CPABORT("Zero degrees of freedom. Nothing to thermalize!") 400 END IF 401 gle%nvt(1)%nkt = simpar%temp_ext*fac 402 gle%nvt(1)%degrees_of_freedom = FLOOR(fac) 403 ELSE 404 DO i = 1, gle%loc_num_gle 405 imap = map_info%map_index(i) 406 fac = (map_info%s_kin(imap) - deg_of_freedom(i)) 407 gle%nvt(i)%nkt = simpar%temp_ext*fac 408 gle%nvt(i)%degrees_of_freedom = FLOOR(fac) 409 END DO 410 END IF 411 DEALLOCATE (deg_of_freedom) 412 DEALLOCATE (massive_atom_list) 413 END SELECT 414 415 END SUBROUTINE gle_to_particle_mapping 416 417! ************************************************************************************************** 418!> \brief ... 419!> \param gle ... 420!> \param gle_section ... 421!> \param save_mem ... 422!> \param restart ... 423! ************************************************************************************************** 424 SUBROUTINE restart_gle(gle, gle_section, save_mem, restart) 425 426 TYPE(gle_type), POINTER :: gle 427 TYPE(section_vals_type), POINTER :: gle_section 428 LOGICAL, INTENT(IN) :: save_mem 429 LOGICAL, INTENT(OUT) :: restart 430 431 CHARACTER(len=*), PARAMETER :: routineN = 'restart_gle', routineP = moduleN//':'//routineN 432 433 CHARACTER(LEN=rng_record_length) :: rng_record 434 INTEGER :: glob_num, i, ind, j, loc_num, n_rep 435 LOGICAL :: explicit 436 REAL(KIND=dp), DIMENSION(:), POINTER :: buffer 437 TYPE(map_info_type), POINTER :: map_info 438 TYPE(section_vals_type), POINTER :: work_section 439 440 NULLIFY (buffer, work_section) 441 442 explicit = .FALSE. 443 restart = .FALSE. 444 445 IF (ASSOCIATED(gle_section)) THEN 446 work_section => section_vals_get_subs_vals(gle_section, "S") 447 CALL section_vals_get(work_section, explicit=explicit) 448 restart = explicit 449 END IF 450 451 IF (restart) THEN 452 map_info => gle%map_info 453 CALL section_vals_val_get(gle_section, "S%_DEFAULT_KEYWORD_", r_vals=buffer) 454 DO i = 1, SIZE(gle%nvt, 1) 455 ind = map_info%index(i) 456 ind = (ind - 1)*(gle%ndim) 457 DO j = 1, SIZE(gle%nvt(i)%s, 1) 458 ind = ind + 1 459 gle%nvt(i)%s(j) = buffer(ind) 460 END DO 461 END DO 462 463 IF (save_mem) THEN 464 NULLIFY (work_section) 465 work_section => section_vals_get_subs_vals(gle_section, "S") 466 CALL section_vals_remove_values(work_section) 467 END IF 468 469 ! Possibly restart the initial thermostat energy 470 work_section => section_vals_get_subs_vals(section_vals=gle_section, & 471 subsection_name="THERMOSTAT_ENERGY") 472 CALL section_vals_get(work_section, explicit=explicit) 473 IF (explicit) THEN 474 CALL section_vals_val_get(section_vals=work_section, keyword_name="_DEFAULT_KEYWORD_", & 475 n_rep_val=n_rep) 476 IF (n_rep == gle%glob_num_gle) THEN 477 DO i = 1, gle%loc_num_gle 478 ind = map_info%index(i) 479 CALL section_vals_val_get(section_vals=work_section, keyword_name="_DEFAULT_KEYWORD_", & 480 i_rep_val=ind, r_val=gle%nvt(i)%thermostat_energy) 481 END DO 482 ELSE 483 CALL cp_abort(__LOCATION__, & 484 "Number of thermostat energies not equal to the number of "// & 485 "total thermostats!") 486 END IF 487 END IF 488 489 ! Possibly restart the random number generators for the different thermostats 490 work_section => section_vals_get_subs_vals(section_vals=gle_section, & 491 subsection_name="RNG_INIT") 492 493 CALL section_vals_get(work_section, explicit=explicit) 494 IF (explicit) THEN 495 CALL section_vals_val_get(section_vals=work_section, keyword_name="_DEFAULT_KEYWORD_", & 496 n_rep_val=n_rep) 497 498 glob_num = gle%glob_num_gle 499 loc_num = gle%loc_num_gle 500 IF (n_rep == glob_num) THEN 501 DO i = 1, loc_num 502 ind = map_info%index(i) 503 CALL section_vals_val_get(section_vals=work_section, keyword_name="_DEFAULT_KEYWORD_", & 504 i_rep_val=ind, c_val=rng_record) 505 CALL read_rng_stream(rng_stream=gle%nvt(i)%gaussian_rng_stream, & 506 rng_record=rng_record) 507 END DO 508 ELSE 509 CALL cp_abort(__LOCATION__, & 510 "Number pf restartable stream not equal to the number of "// & 511 "total thermostats!") 512 END IF 513 END IF 514 END IF 515 516 END SUBROUTINE restart_gle 517 518! ************************************************************************************************** 519!> \brief ... 520!> \param gle ... 521! ************************************************************************************************** 522 SUBROUTINE init_gle_variables(gle) 523 524 TYPE(gle_type), POINTER :: gle 525 526 CHARACTER(len=*), PARAMETER :: routineN = 'init_gle_variables', & 527 routineP = moduleN//':'//routineN 528 529 INTEGER :: i, j 530 REAL(dp) :: rr(gle%ndim), cc(gle%ndim, gle%ndim) 531 TYPE(rng_stream_type), POINTER :: rng_stream 532 533 CALL gle_cholesky_stab(gle%c_mat, cc, gle%ndim) 534 DO i = 1, gle%loc_num_gle 535 rng_stream => gle%nvt(i)%gaussian_rng_stream 536 DO j = 1, gle%ndim 537 ! here s should be properly initialized, when it is not read from restart 538 rr(j) = next_random_number(rng_stream) 539 END DO 540 gle%nvt(i)%s = MATMUL(cc, rr) 541 END DO 542 543 END SUBROUTINE init_gle_variables 544 545END MODULE gle_system_dynamics 546