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!> CJM, 20-Feb-01 9!> JGH (10-Mar-2001) 10!> CJM (10-Apr-2001) 11!> \author CJM 12! ************************************************************************************************** 13MODULE extended_system_mapping 14 15 USE cp_para_types, ONLY: cp_para_env_type 16 USE distribution_1d_types, ONLY: distribution_1d_type 17 USE extended_system_types, ONLY: debug_isotropic_limit,& 18 lnhc_parameters_type,& 19 map_info_type 20 USE input_constants, ONLY: & 21 do_thermo_communication, do_thermo_no_communication, do_thermo_only_master, & 22 isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, & 23 nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, & 24 nve_ensemble, nvt_adiabatic_ensemble, nvt_ensemble, reftraj_ensemble 25 USE kinds, ONLY: dp 26 USE message_passing, ONLY: mp_sum 27 USE molecule_kind_types, ONLY: molecule_kind_type 28 USE molecule_types, ONLY: global_constraint_type,& 29 molecule_type 30 USE simpar_types, ONLY: simpar_type 31 USE thermostat_mapping, ONLY: adiabatic_mapping_region,& 32 init_baro_map_info,& 33 thermostat_mapping_region 34 USE thermostat_types, ONLY: thermostat_info_type 35#include "../../base/base_uses.f90" 36 37 IMPLICIT NONE 38 39 PRIVATE 40 41 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'extended_system_mapping' 42 43 PUBLIC :: nhc_to_particle_mapping, nhc_to_barostat_mapping, & 44 nhc_to_shell_mapping, nhc_to_particle_mapping_fast, & 45 nhc_to_particle_mapping_slow 46 47CONTAINS 48 49! ************************************************************************************************** 50!> \brief Creates the thermostatting for the barostat 51!> \param simpar ... 52!> \param nhc ... 53!> \par History 54!> CJM, 20-Feb-01 : nhc structure allocated to zero when not in use 55!> JGH (10-Mar-2001) : set nhc variables to zero when not in use 56!> \author CJM 57! ************************************************************************************************** 58 SUBROUTINE nhc_to_barostat_mapping(simpar, nhc) 59 60 TYPE(simpar_type), POINTER :: simpar 61 TYPE(lnhc_parameters_type), POINTER :: nhc 62 63 CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_barostat_mapping', & 64 routineP = moduleN//':'//routineN 65 66 INTEGER :: handle, i, number 67 TYPE(map_info_type), POINTER :: map_info 68 69 CALL timeset(routineN, handle) 70 71 SELECT CASE (simpar%ensemble) 72 CASE DEFAULT 73 CPABORT('Never reach this point!') 74 CASE (npt_i_ensemble, npt_f_ensemble) 75 map_info => nhc%map_info 76 map_info%dis_type = do_thermo_only_master 77 78 ! Counting the total number of thermostats ( 1 for both NPT_I and NPT_F ) 79 nhc%loc_num_nhc = 1 80 nhc%glob_num_nhc = 1 81 IF (simpar%ensemble == npt_f_ensemble) THEN 82 number = 9 83 ELSE 84 number = 1 85 ENDIF 86 87 CALL init_baro_map_info(map_info, number, nhc%loc_num_nhc) 88 89 ALLOCATE (nhc%nvt(nhc%nhc_len, nhc%loc_num_nhc)) 90 ! Now that we know how many there are stick this into nhc % nkt 91 ! (number of degrees of freedom times k_B T ) 92 DO i = 1, nhc%loc_num_nhc 93 nhc%nvt(1, i)%nkt = simpar%temp_ext*number 94 nhc%nvt(1, i)%degrees_of_freedom = number 95 IF (debug_isotropic_limit) THEN 96 nhc%nvt(1, i)%nkt = simpar%temp_ext 97 END IF 98 END DO 99 100 ! getting the number of degrees of freedom times k_B T for the rest of the chain 101 DO i = 2, nhc%nhc_len 102 nhc%nvt(i, :)%nkt = simpar%temp_ext 103 END DO 104 105 ! Let's clean the arrays 106 map_info%s_kin = 0.0_dp 107 map_info%v_scale = 0.0_dp 108 END SELECT 109 110 CALL timestop(handle) 111 112 END SUBROUTINE nhc_to_barostat_mapping 113 114! ************************************************************************************************** 115!> \brief Creates the thermostatting maps 116!> \param thermostat_info ... 117!> \param simpar ... 118!> \param local_molecules ... 119!> \param molecule_set ... 120!> \param molecule_kind_set ... 121!> \param nhc ... 122!> \param para_env ... 123!> \param gci ... 124!> \par History 125!> 29-Nov-00 (JGH) correct counting of DOF if constraints are off 126!> CJM, 20-Feb-01 : nhc structure allocated to zero when not in use 127!> JGH (10-Mar-2001) : set nhc variables to zero when not in use 128!> CJM(10-NOV-2001) : New parallelization with new molecule structures 129!> Teodoro Laino 09.2007 [tlaino] - University of Zurich - cleaning and updating 130!> \author CJM 131! ************************************************************************************************** 132 SUBROUTINE nhc_to_particle_mapping(thermostat_info, simpar, local_molecules, & 133 molecule_set, molecule_kind_set, nhc, para_env, gci) 134 135 TYPE(thermostat_info_type), POINTER :: thermostat_info 136 TYPE(simpar_type), POINTER :: simpar 137 TYPE(distribution_1d_type), POINTER :: local_molecules 138 TYPE(molecule_type), POINTER :: molecule_set(:) 139 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 140 TYPE(lnhc_parameters_type), POINTER :: nhc 141 TYPE(cp_para_env_type), POINTER :: para_env 142 TYPE(global_constraint_type), POINTER :: gci 143 144 CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_particle_mapping', & 145 routineP = moduleN//':'//routineN 146 147 INTEGER :: handle, i, imap, j, natoms_local, & 148 sum_of_thermostats 149 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list 150 REAL(KIND=dp) :: fac 151 TYPE(map_info_type), POINTER :: map_info 152 153 CALL timeset(routineN, handle) 154 155 NULLIFY (massive_atom_list, deg_of_freedom) 156 157 SELECT CASE (simpar%ensemble) 158 CASE DEFAULT 159 CPABORT('Unknown ensemble!') 160 CASE (nve_ensemble, isokin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_ensemble, & 161 nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble) 162 CPABORT('Never reach this point!') 163 CASE (nvt_ensemble, npt_i_ensemble, npt_f_ensemble) 164 165 CALL setup_nhc_thermostat(nhc, thermostat_info, deg_of_freedom, massive_atom_list, & 166 molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, & 167 simpar, sum_of_thermostats, gci) 168 169 ! Sum up the number of degrees of freedom on each thermostat. 170 ! first: initialize the target 171 map_info => nhc%map_info 172 map_info%s_kin = 0.0_dp 173 DO i = 1, 3 174 DO j = 1, natoms_local 175 map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1 176 END DO 177 END DO 178 179 ! if thermostats are replicated but molecules distributed, we have to 180 ! sum s_kin over all processors 181 IF (map_info%dis_type == do_thermo_communication) CALL mp_sum(map_info%s_kin, para_env%group) 182 183 ! We know the total number of system thermostats. 184 IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN 185 fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl 186 IF (fac == 0.0_dp) THEN 187 CPABORT('Zero degrees of freedom. Nothing to thermalize!') 188 END IF 189 nhc%nvt(1, 1)%nkt = simpar%temp_ext*fac 190 nhc%nvt(1, 1)%degrees_of_freedom = FLOOR(fac) 191 ELSE 192 DO i = 1, nhc%loc_num_nhc 193 imap = map_info%map_index(i) 194 fac = (map_info%s_kin(imap) - deg_of_freedom(i)) 195 nhc%nvt(1, i)%nkt = simpar%temp_ext*fac 196 nhc%nvt(1, i)%degrees_of_freedom = FLOOR(fac) 197 END DO 198 END IF 199 200 ! Getting the number of degrees of freedom times k_B T for the rest 201 ! of the chain 202 DO i = 2, nhc%nhc_len 203 nhc%nvt(i, :)%nkt = simpar%temp_ext 204 nhc%nvt(i, :)%degrees_of_freedom = 1 205 END DO 206 DEALLOCATE (deg_of_freedom) 207 DEALLOCATE (massive_atom_list) 208 209 ! Let's clean the arrays 210 map_info%s_kin = 0.0_dp 211 map_info%v_scale = 0.0_dp 212 END SELECT 213 214 CALL timestop(handle) 215 216 END SUBROUTINE nhc_to_particle_mapping 217 218! ************************************************************************************************** 219!> \brief Main general setup for Adiabatic Nose-Hoover thermostats 220!> \param nhc ... 221!> \param thermostat_info ... 222!> \param deg_of_freedom ... 223!> \param massive_atom_list ... 224!> \param molecule_kind_set ... 225!> \param local_molecules ... 226!> \param molecule_set ... 227!> \param para_env ... 228!> \param natoms_local ... 229!> \param simpar ... 230!> \param sum_of_thermostats ... 231!> \param gci ... 232!> \param shell ... 233!> \author CJM -PNNL -2011 234! ************************************************************************************************** 235 SUBROUTINE setup_adiabatic_thermostat(nhc, thermostat_info, deg_of_freedom, & 236 massive_atom_list, molecule_kind_set, local_molecules, molecule_set, & 237 para_env, natoms_local, simpar, sum_of_thermostats, gci, shell) 238 239 TYPE(lnhc_parameters_type), POINTER :: nhc 240 TYPE(thermostat_info_type), POINTER :: thermostat_info 241 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list 242 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 243 TYPE(distribution_1d_type), POINTER :: local_molecules 244 TYPE(molecule_type), POINTER :: molecule_set(:) 245 TYPE(cp_para_env_type), POINTER :: para_env 246 INTEGER, INTENT(OUT) :: natoms_local 247 TYPE(simpar_type), POINTER :: simpar 248 INTEGER, INTENT(OUT) :: sum_of_thermostats 249 TYPE(global_constraint_type), POINTER :: gci 250 LOGICAL, INTENT(IN), OPTIONAL :: shell 251 252 CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_adiabatic_thermostat', & 253 routineP = moduleN//':'//routineN 254 255 INTEGER :: handle, nkind, number, region 256 LOGICAL :: do_shell 257 TYPE(map_info_type), POINTER :: map_info 258 259 CALL timeset(routineN, handle) 260 261 do_shell = .FALSE. 262 IF (PRESENT(shell)) do_shell = shell 263 map_info => nhc%map_info 264 265 nkind = SIZE(molecule_kind_set) 266 sum_of_thermostats = thermostat_info%sum_of_thermostats 267 map_info%dis_type = thermostat_info%dis_type 268 number = thermostat_info%number_of_thermostats 269 region = nhc%region 270 271 CALL adiabatic_mapping_region(map_info, deg_of_freedom, massive_atom_list, & 272 molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, & 273 simpar, number, region, gci, do_shell, thermostat_info%map_loc_thermo_gen, & 274 sum_of_thermostats) 275 ALLOCATE (nhc%nvt(nhc%nhc_len, number)) 276 277 ! Now that we know how many there are stick this into nhc%nkt 278 ! (number of degrees of freedom times k_B T for the first thermostat 279 ! on the chain) 280 nhc%loc_num_nhc = number 281 nhc%glob_num_nhc = sum_of_thermostats 282 283 CALL timestop(handle) 284 285 END SUBROUTINE setup_adiabatic_thermostat 286 287! ************************************************************************************************** 288!> \brief Creates the thermostatting maps 289!> \param thermostat_info ... 290!> \param simpar ... 291!> \param local_molecules ... 292!> \param molecule_set ... 293!> \param molecule_kind_set ... 294!> \param nhc ... 295!> \param para_env ... 296!> \param gci ... 297!> \par History 298!> \author CJM 299! ************************************************************************************************** 300 SUBROUTINE nhc_to_particle_mapping_slow(thermostat_info, simpar, local_molecules, & 301 molecule_set, molecule_kind_set, nhc, para_env, gci) 302 303 TYPE(thermostat_info_type), POINTER :: thermostat_info 304 TYPE(simpar_type), POINTER :: simpar 305 TYPE(distribution_1d_type), POINTER :: local_molecules 306 TYPE(molecule_type), POINTER :: molecule_set(:) 307 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 308 TYPE(lnhc_parameters_type), POINTER :: nhc 309 TYPE(cp_para_env_type), POINTER :: para_env 310 TYPE(global_constraint_type), POINTER :: gci 311 312 CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_particle_mapping_slow', & 313 routineP = moduleN//':'//routineN 314 315 INTEGER :: handle, i, imap, j, natoms_local, & 316 sum_of_thermostats 317 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list 318 REAL(KIND=dp) :: fac 319 TYPE(map_info_type), POINTER :: map_info 320 321 CALL timeset(routineN, handle) 322 323 NULLIFY (massive_atom_list, deg_of_freedom) 324 325 SELECT CASE (simpar%ensemble) 326 CASE DEFAULT 327 CPABORT('Unknown ensemble!') 328 CASE (nvt_adiabatic_ensemble) 329 CALL setup_adiabatic_thermostat(nhc, thermostat_info, deg_of_freedom, massive_atom_list, & 330 molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, & 331 simpar, sum_of_thermostats, gci) 332 333 ! Sum up the number of degrees of freedom on each thermostat. 334 ! first: initialize the target 335 map_info => nhc%map_info 336 map_info%s_kin = 0.0_dp 337 DO i = 1, 3 338 DO j = 1, natoms_local 339 IF (ASSOCIATED(map_info%p_kin(i, j)%point)) & 340 map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1 341 END DO 342 END DO 343 344 ! if thermostats are replicated but molecules distributed, we have to 345 ! sum s_kin over all processors 346 IF (map_info%dis_type == do_thermo_communication) CALL mp_sum(map_info%s_kin, para_env%group) 347 348 ! We know the total number of system thermostats. 349 IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN 350 fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl 351 IF (fac == 0.0_dp) THEN 352 CPABORT('Zero degrees of freedom. Nothing to thermalize!') 353 END IF 354 nhc%nvt(1, 1)%nkt = simpar%temp_slow*fac 355 nhc%nvt(1, 1)%degrees_of_freedom = FLOOR(fac) 356 ELSE 357 DO i = 1, nhc%loc_num_nhc 358 imap = map_info%map_index(i) 359 fac = (map_info%s_kin(imap) - deg_of_freedom(i)) 360 nhc%nvt(1, i)%nkt = simpar%temp_slow*fac 361 nhc%nvt(1, i)%degrees_of_freedom = FLOOR(fac) 362 END DO 363 END IF 364 365 ! Getting the number of degrees of freedom times k_B T for the rest 366 ! of the chain 367 DO i = 2, nhc%nhc_len 368 nhc%nvt(i, :)%nkt = simpar%temp_slow 369 nhc%nvt(i, :)%degrees_of_freedom = 1 370 END DO 371 DEALLOCATE (deg_of_freedom) 372 DEALLOCATE (massive_atom_list) 373 374 ! Let's clean the arrays 375 map_info%s_kin = 0.0_dp 376 map_info%v_scale = 0.0_dp 377 END SELECT 378 379 CALL timestop(handle) 380 381 END SUBROUTINE nhc_to_particle_mapping_slow 382 383! ************************************************************************************************** 384!> \brief Creates the thermostatting maps 385!> \param thermostat_info ... 386!> \param simpar ... 387!> \param local_molecules ... 388!> \param molecule_set ... 389!> \param molecule_kind_set ... 390!> \param nhc ... 391!> \param para_env ... 392!> \param gci ... 393!> \par History 394!> \author CJM 395! ************************************************************************************************** 396 SUBROUTINE nhc_to_particle_mapping_fast(thermostat_info, simpar, local_molecules, & 397 molecule_set, molecule_kind_set, nhc, para_env, gci) 398 399 TYPE(thermostat_info_type), POINTER :: thermostat_info 400 TYPE(simpar_type), POINTER :: simpar 401 TYPE(distribution_1d_type), POINTER :: local_molecules 402 TYPE(molecule_type), POINTER :: molecule_set(:) 403 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 404 TYPE(lnhc_parameters_type), POINTER :: nhc 405 TYPE(cp_para_env_type), POINTER :: para_env 406 TYPE(global_constraint_type), POINTER :: gci 407 408 CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_particle_mapping_fast', & 409 routineP = moduleN//':'//routineN 410 411 INTEGER :: handle, i, imap, j, natoms_local, & 412 sum_of_thermostats 413 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list 414 REAL(KIND=dp) :: fac 415 TYPE(map_info_type), POINTER :: map_info 416 417 CALL timeset(routineN, handle) 418 419 NULLIFY (massive_atom_list, deg_of_freedom) 420 421 SELECT CASE (simpar%ensemble) 422 CASE DEFAULT 423 CPABORT('Unknown ensemble!') 424 CASE (nvt_adiabatic_ensemble) 425 CALL setup_adiabatic_thermostat(nhc, thermostat_info, deg_of_freedom, massive_atom_list, & 426 molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, & 427 simpar, sum_of_thermostats, gci) 428 429 ! Sum up the number of degrees of freedom on each thermostat. 430 ! first: initialize the target 431 map_info => nhc%map_info 432 map_info%s_kin = 0.0_dp 433 DO i = 1, 3 434 DO j = 1, natoms_local 435 IF (ASSOCIATED(map_info%p_kin(i, j)%point)) & 436 map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1 437 END DO 438 END DO 439 440 ! if thermostats are replicated but molecules distributed, we have to 441 ! sum s_kin over all processors 442 IF (map_info%dis_type == do_thermo_communication) CALL mp_sum(map_info%s_kin, para_env%group) 443 444 ! We know the total number of system thermostats. 445 IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN 446 fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl 447 IF (fac == 0.0_dp) THEN 448 CPABORT('Zero degrees of freedom. Nothing to thermalize!') 449 END IF 450 nhc%nvt(1, 1)%nkt = simpar%temp_fast*fac 451 nhc%nvt(1, 1)%degrees_of_freedom = FLOOR(fac) 452 ELSE 453 DO i = 1, nhc%loc_num_nhc 454 imap = map_info%map_index(i) 455 fac = (map_info%s_kin(imap) - deg_of_freedom(i)) 456 nhc%nvt(1, i)%nkt = simpar%temp_fast*fac 457 nhc%nvt(1, i)%degrees_of_freedom = FLOOR(fac) 458 END DO 459 END IF 460 461 ! Getting the number of degrees of freedom times k_B T for the rest 462 ! of the chain 463 DO i = 2, nhc%nhc_len 464 nhc%nvt(i, :)%nkt = simpar%temp_fast 465 nhc%nvt(i, :)%degrees_of_freedom = 1 466 END DO 467 DEALLOCATE (deg_of_freedom) 468 DEALLOCATE (massive_atom_list) 469 470 ! Let's clean the arrays 471 map_info%s_kin = 0.0_dp 472 map_info%v_scale = 0.0_dp 473 END SELECT 474 475 CALL timestop(handle) 476 477 END SUBROUTINE nhc_to_particle_mapping_fast 478 479! ************************************************************************************************** 480!> \brief Main general setup for Nose-Hoover thermostats 481!> \param nhc ... 482!> \param thermostat_info ... 483!> \param deg_of_freedom ... 484!> \param massive_atom_list ... 485!> \param molecule_kind_set ... 486!> \param local_molecules ... 487!> \param molecule_set ... 488!> \param para_env ... 489!> \param natoms_local ... 490!> \param simpar ... 491!> \param sum_of_thermostats ... 492!> \param gci ... 493!> \param shell ... 494!> \author Teodoro Laino [tlaino] - University of Zurich - 10.2007 495! ************************************************************************************************** 496 SUBROUTINE setup_nhc_thermostat(nhc, thermostat_info, deg_of_freedom, & 497 massive_atom_list, molecule_kind_set, local_molecules, molecule_set, & 498 para_env, natoms_local, simpar, sum_of_thermostats, gci, shell) 499 500 TYPE(lnhc_parameters_type), POINTER :: nhc 501 TYPE(thermostat_info_type), POINTER :: thermostat_info 502 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_atom_list 503 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 504 TYPE(distribution_1d_type), POINTER :: local_molecules 505 TYPE(molecule_type), POINTER :: molecule_set(:) 506 TYPE(cp_para_env_type), POINTER :: para_env 507 INTEGER, INTENT(OUT) :: natoms_local 508 TYPE(simpar_type), POINTER :: simpar 509 INTEGER, INTENT(OUT) :: sum_of_thermostats 510 TYPE(global_constraint_type), POINTER :: gci 511 LOGICAL, INTENT(IN), OPTIONAL :: shell 512 513 CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_nhc_thermostat', & 514 routineP = moduleN//':'//routineN 515 516 INTEGER :: handle, nkind, number, region 517 LOGICAL :: do_shell 518 TYPE(map_info_type), POINTER :: map_info 519 520 CALL timeset(routineN, handle) 521 522 do_shell = .FALSE. 523 IF (PRESENT(shell)) do_shell = shell 524 map_info => nhc%map_info 525 526 nkind = SIZE(molecule_kind_set) 527 sum_of_thermostats = thermostat_info%sum_of_thermostats 528 map_info%dis_type = thermostat_info%dis_type 529 number = thermostat_info%number_of_thermostats 530 region = nhc%region 531 532 CALL thermostat_mapping_region(map_info, deg_of_freedom, massive_atom_list, & 533 molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, & 534 simpar, number, region, gci, do_shell, thermostat_info%map_loc_thermo_gen, & 535 sum_of_thermostats) 536 537 ALLOCATE (nhc%nvt(nhc%nhc_len, number)) 538 539 ! Now that we know how many there are stick this into nhc%nkt 540 ! (number of degrees of freedom times k_B T for the first thermostat 541 ! on the chain) 542 nhc%loc_num_nhc = number 543 nhc%glob_num_nhc = sum_of_thermostats 544 545 CALL timestop(handle) 546 547 END SUBROUTINE setup_nhc_thermostat 548 549! ************************************************************************************************** 550!> \brief ... 551!> \param thermostat_info ... 552!> \param simpar ... 553!> \param local_molecules ... 554!> \param molecule_set ... 555!> \param molecule_kind_set ... 556!> \param nhc ... 557!> \param para_env ... 558!> \param gci ... 559! ************************************************************************************************** 560 SUBROUTINE nhc_to_shell_mapping(thermostat_info, simpar, local_molecules, & 561 molecule_set, molecule_kind_set, nhc, para_env, gci) 562 563 TYPE(thermostat_info_type), POINTER :: thermostat_info 564 TYPE(simpar_type), POINTER :: simpar 565 TYPE(distribution_1d_type), POINTER :: local_molecules 566 TYPE(molecule_type), POINTER :: molecule_set(:) 567 TYPE(molecule_kind_type), POINTER :: molecule_kind_set(:) 568 TYPE(lnhc_parameters_type), POINTER :: nhc 569 TYPE(cp_para_env_type), POINTER :: para_env 570 TYPE(global_constraint_type), POINTER :: gci 571 572 CHARACTER(LEN=*), PARAMETER :: routineN = 'nhc_to_shell_mapping', & 573 routineP = moduleN//':'//routineN 574 575 INTEGER :: handle, i, imap, j, nshell_local, & 576 sum_of_thermostats 577 INTEGER, DIMENSION(:), POINTER :: deg_of_freedom, massive_shell_list 578 TYPE(map_info_type), POINTER :: map_info 579 580 CALL timeset(routineN, handle) 581 582 NULLIFY (massive_shell_list, deg_of_freedom) 583 584 SELECT CASE (simpar%ensemble) 585 CASE DEFAULT 586 CPABORT('Unknown ensemble!') 587 CASE (isokin_ensemble, nph_uniaxial_ensemble, & 588 nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble) 589 CPABORT('Never reach this point!') 590 CASE (nve_ensemble, nvt_ensemble, npe_f_ensemble, npe_i_ensemble, npt_i_ensemble, npt_f_ensemble) 591 592 CALL setup_nhc_thermostat(nhc, thermostat_info, deg_of_freedom, massive_shell_list, & 593 molecule_kind_set, local_molecules, molecule_set, para_env, nshell_local, & 594 simpar, sum_of_thermostats, gci, shell=.TRUE.) 595 596 map_info => nhc%map_info 597 ! Sum up the number of degrees of freedom on each thermostat. 598 ! first: initialize the target, via p_kin init s_kin 599 map_info%s_kin = 0.0_dp 600 DO j = 1, nshell_local 601 DO i = 1, 3 602 map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1 603 END DO 604 END DO 605 606 ! If thermostats are replicated but molecules distributed, we have to 607 ! sum s_kin over all processors 608 IF (map_info%dis_type == do_thermo_communication) CALL mp_sum(map_info%s_kin, para_env%group) 609 610 ! Now that we know how many there are stick this into nhc%nkt 611 ! (number of degrees of freedom times k_B T ) 612 DO i = 1, nhc%loc_num_nhc 613 imap = map_info%map_index(i) 614 nhc%nvt(1, i)%nkt = simpar%temp_sh_ext*map_info%s_kin(imap) 615 nhc%nvt(1, i)%degrees_of_freedom = INT(map_info%s_kin(imap)) 616 END DO 617 618 ! Getting the number of degrees of freedom times k_B T for the rest of the chain 619 DO i = 2, nhc%nhc_len 620 nhc%nvt(i, :)%nkt = simpar%temp_sh_ext 621 nhc%nvt(i, :)%degrees_of_freedom = 1 622 END DO 623 DEALLOCATE (deg_of_freedom) 624 DEALLOCATE (massive_shell_list) 625 626 ! Let's clean the arrays 627 map_info%s_kin = 0.0_dp 628 map_info%v_scale = 0.0_dp 629 END SELECT 630 631 CALL timestop(handle) 632 633 END SUBROUTINE nhc_to_shell_mapping 634 635END MODULE extended_system_mapping 636