Line data    Source code

       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \author Teodoro Laino [tlaino] 10.2007- University of Zurich
      10             : ! **************************************************************************************************
      11             : MODULE csvr_system_mapping
      12             : 
      13             :    USE csvr_system_types,               ONLY: csvr_system_type,&
      14             :                                               csvr_thermo_create
      15             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      16             :    USE extended_system_types,           ONLY: debug_isotropic_limit,&
      17             :                                               map_info_type
      18             :    USE input_constants,                 ONLY: &
      19             :         do_thermo_communication, do_thermo_no_communication, do_thermo_only_master, &
      20             :         isokin_ensemble, langevin_ensemble, npe_f_ensemble, npe_i_ensemble, &
      21             :         nph_uniaxial_damped_ensemble, nph_uniaxial_ensemble, npt_f_ensemble, npt_i_ensemble, &
      22             :         npt_ia_ensemble, nve_ensemble, nvt_ensemble, reftraj_ensemble
      23             :    USE kinds,                           ONLY: dp
      24             :    USE message_passing,                 ONLY: mp_para_env_type
      25             :    USE molecule_kind_types,             ONLY: molecule_kind_type
      26             :    USE molecule_types,                  ONLY: global_constraint_type,&
      27             :                                               molecule_type
      28             :    USE simpar_types,                    ONLY: simpar_type
      29             :    USE thermostat_mapping,              ONLY: init_baro_map_info,&
      30             :                                               thermostat_mapping_region
      31             :    USE thermostat_types,                ONLY: thermostat_info_type
      32             : #include "../../base/base_uses.f90"
      33             : 
      34             :    IMPLICIT NONE
      35             : 
      36             :    PRIVATE
      37             : 
      38             :    ! *** Global parameters ***
      39             : 
      40             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'csvr_system_mapping'
      41             : 
      42             :    PUBLIC :: csvr_to_particle_mapping, csvr_to_barostat_mapping, &
      43             :              csvr_to_shell_mapping
      44             : 
      45             : CONTAINS
      46             : 
      47             : ! **************************************************************************************************
      48             : !> \brief Creates the thermostatting for the barostat
      49             : !> \param simpar ...
      50             : !> \param csvr ...
      51             : !> \author Teodoro Laino [tlaino] 10.2007- University of Zurich
      52             : ! **************************************************************************************************
      53          32 :    SUBROUTINE csvr_to_barostat_mapping(simpar, csvr)
      54             :       TYPE(simpar_type), POINTER                         :: simpar
      55             :       TYPE(csvr_system_type), POINTER                    :: csvr
      56             : 
      57             :       INTEGER                                            :: i, ndeg
      58             :       TYPE(map_info_type), POINTER                       :: map_info
      59             : 
      60          32 :       SELECT CASE (simpar%ensemble)
      61             :       CASE DEFAULT
      62           0 :          CPABORT('Never reach this point!')
      63             :       CASE (npt_i_ensemble, npt_f_ensemble, npt_ia_ensemble)
      64          32 :          map_info => csvr%map_info
      65          32 :          map_info%dis_type = do_thermo_only_master
      66             : 
      67             :          ! Counting the total number of thermostats ( 1 for NPT_I, NPT_IA, and NPT_F )
      68          32 :          csvr%loc_num_csvr = 1
      69          32 :          csvr%glob_num_csvr = 1
      70          32 :          IF (simpar%ensemble == npt_f_ensemble) THEN
      71           4 :             ndeg = 9
      72             :          ELSE
      73          28 :             ndeg = 1
      74             :          END IF
      75             : 
      76          32 :          CALL init_baro_map_info(map_info, ndeg, csvr%loc_num_csvr)
      77          32 :          CALL csvr_thermo_create(csvr)
      78             : 
      79             :          ! Now that we know how many there are stick this into csvr%nkt
      80             :          ! (number of degrees of freedom times k_B T )
      81          96 :          DO i = 1, csvr%loc_num_csvr
      82          32 :             csvr%nvt(i)%nkt = simpar%temp_baro_ext*ndeg
      83          32 :             csvr%nvt(i)%degrees_of_freedom = ndeg
      84          32 :             IF (debug_isotropic_limit) THEN
      85             :                csvr%nvt(i)%nkt = simpar%temp_baro_ext
      86             :                csvr%nvt(i)%degrees_of_freedom = 1
      87             :             END IF
      88             :          END DO
      89             :       END SELECT
      90             : 
      91          32 :    END SUBROUTINE csvr_to_barostat_mapping
      92             : 
      93             : ! **************************************************************************************************
      94             : !> \brief Creates the thermostatting maps
      95             : !> \param thermostat_info ...
      96             : !> \param simpar ...
      97             : !> \param local_molecules ...
      98             : !> \param molecule_set ...
      99             : !> \param molecule_kind_set ...
     100             : !> \param csvr ...
     101             : !> \param para_env ...
     102             : !> \param gci ...
     103             : !> \author Teodoro Laino [tlaino] 10.2007- University of Zurich
     104             : ! **************************************************************************************************
     105         122 :    SUBROUTINE csvr_to_particle_mapping(thermostat_info, simpar, local_molecules, &
     106             :                                        molecule_set, molecule_kind_set, csvr, para_env, gci)
     107             : 
     108             :       TYPE(thermostat_info_type), POINTER                :: thermostat_info
     109             :       TYPE(simpar_type), POINTER                         :: simpar
     110             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     111             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     112             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     113             :       TYPE(csvr_system_type), POINTER                    :: csvr
     114             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     115             :       TYPE(global_constraint_type), POINTER              :: gci
     116             : 
     117             :       INTEGER                                            :: i, imap, j, natoms_local, &
     118             :                                                             sum_of_thermostats
     119         122 :       INTEGER, DIMENSION(:), POINTER                     :: deg_of_freedom, massive_atom_list
     120             :       REAL(KIND=dp)                                      :: fac
     121             :       TYPE(map_info_type), POINTER                       :: map_info
     122             : 
     123         122 :       NULLIFY (massive_atom_list, deg_of_freedom)
     124         122 :       SELECT CASE (simpar%ensemble)
     125             :       CASE DEFAULT
     126           0 :          CPABORT('Unknown ensemble!')
     127             :       CASE (nve_ensemble, isokin_ensemble, npe_f_ensemble, npe_i_ensemble, nph_uniaxial_ensemble, &
     128             :             nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble)
     129           0 :          CPABORT('Never reach this point!')
     130             :       CASE (nvt_ensemble, npt_i_ensemble, npt_f_ensemble, npt_ia_ensemble)
     131             : 
     132             :          CALL setup_csvr_thermostat(csvr, thermostat_info, deg_of_freedom, &
     133             :                                     massive_atom_list, molecule_kind_set, local_molecules, molecule_set, &
     134         122 :                                     para_env, natoms_local, simpar, sum_of_thermostats, gci)
     135             : 
     136             :          ! Sum up the number of degrees of freedom on each thermostat.
     137             :          ! first: initialize the target
     138         122 :          map_info => csvr%map_info
     139        1910 :          map_info%s_kin = 0.0_dp
     140         488 :          DO i = 1, 3
     141       89012 :             DO j = 1, natoms_local
     142       88890 :                map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
     143             :             END DO
     144             :          END DO
     145             : 
     146             :          ! If thermostats are replicated but molecules distributed, we have to
     147             :          ! sum s_kin over all processors
     148         354 :          IF (map_info%dis_type == do_thermo_communication) CALL para_env%sum(map_info%s_kin)
     149             : 
     150             :          ! We know the total number of system thermostats.
     151         122 :          IF ((sum_of_thermostats == 1) .AND. (map_info%dis_type /= do_thermo_no_communication)) THEN
     152          62 :             fac = map_info%s_kin(1) - deg_of_freedom(1) - simpar%nfree_rot_transl
     153          62 :             IF (fac == 0.0_dp) THEN
     154           0 :                CPABORT('Zero degrees of freedom. Nothing to thermalize!')
     155             :             END IF
     156          62 :             csvr%nvt(1)%nkt = simpar%temp_ext*fac
     157          62 :             csvr%nvt(1)%degrees_of_freedom = FLOOR(fac)
     158             :          ELSE
     159        1778 :             DO i = 1, csvr%loc_num_csvr
     160        1718 :                imap = map_info%map_index(i)
     161        1718 :                fac = (map_info%s_kin(imap) - deg_of_freedom(i))
     162        1718 :                csvr%nvt(i)%nkt = simpar%temp_ext*fac
     163        1778 :                csvr%nvt(i)%degrees_of_freedom = FLOOR(fac)
     164             :             END DO
     165             :          END IF
     166             : 
     167         122 :          DEALLOCATE (deg_of_freedom)
     168         244 :          DEALLOCATE (massive_atom_list)
     169             :       END SELECT
     170             : 
     171         122 :    END SUBROUTINE csvr_to_particle_mapping
     172             : 
     173             : ! **************************************************************************************************
     174             : !> \brief Main general setup for CSVR thermostats
     175             : !> \param csvr ...
     176             : !> \param thermostat_info ...
     177             : !> \param deg_of_freedom ...
     178             : !> \param massive_atom_list ...
     179             : !> \param molecule_kind_set ...
     180             : !> \param local_molecules ...
     181             : !> \param molecule_set ...
     182             : !> \param para_env ...
     183             : !> \param natoms_local ...
     184             : !> \param simpar ...
     185             : !> \param sum_of_thermostats ...
     186             : !> \param gci ...
     187             : !> \param shell ...
     188             : !> \author Teodoro Laino [tlaino] - University of Zurich - 10.2007
     189             : ! **************************************************************************************************
     190         256 :    SUBROUTINE setup_csvr_thermostat(csvr, thermostat_info, deg_of_freedom, &
     191             :                                     massive_atom_list, molecule_kind_set, local_molecules, molecule_set, &
     192             :                                     para_env, natoms_local, simpar, sum_of_thermostats, gci, shell)
     193             : 
     194             :       TYPE(csvr_system_type), POINTER                    :: csvr
     195             :       TYPE(thermostat_info_type), POINTER                :: thermostat_info
     196             :       INTEGER, DIMENSION(:), POINTER                     :: deg_of_freedom, massive_atom_list
     197             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     198             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     199             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     200             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     201             :       INTEGER, INTENT(OUT)                               :: natoms_local
     202             :       TYPE(simpar_type), POINTER                         :: simpar
     203             :       INTEGER, INTENT(OUT)                               :: sum_of_thermostats
     204             :       TYPE(global_constraint_type), POINTER              :: gci
     205             :       LOGICAL, INTENT(IN), OPTIONAL                      :: shell
     206             : 
     207             :       INTEGER                                            :: nkind, number, region
     208             :       LOGICAL                                            :: do_shell
     209             :       TYPE(map_info_type), POINTER                       :: map_info
     210             : 
     211         128 :       do_shell = .FALSE.
     212         128 :       IF (PRESENT(shell)) do_shell = shell
     213         128 :       map_info => csvr%map_info
     214             : 
     215         128 :       nkind = SIZE(molecule_kind_set)
     216         128 :       sum_of_thermostats = thermostat_info%sum_of_thermostats
     217         128 :       map_info%dis_type = thermostat_info%dis_type
     218         128 :       number = thermostat_info%number_of_thermostats
     219         128 :       region = csvr%region
     220             : 
     221             :       CALL thermostat_mapping_region(map_info, deg_of_freedom, massive_atom_list, &
     222             :                                      molecule_kind_set, local_molecules, molecule_set, para_env, natoms_local, &
     223             :                                      simpar, number, region, gci, do_shell, thermostat_info%map_loc_thermo_gen, &
     224         128 :                                      sum_of_thermostats)
     225             : 
     226             :       ! This is the local number of available thermostats
     227         128 :       csvr%loc_num_csvr = number
     228         128 :       csvr%glob_num_csvr = sum_of_thermostats
     229         128 :       CALL csvr_thermo_create(csvr)
     230             : 
     231         128 :    END SUBROUTINE setup_csvr_thermostat
     232             : 
     233             : ! **************************************************************************************************
     234             : !> \brief ...
     235             : !> \param thermostat_info ...
     236             : !> \param simpar ...
     237             : !> \param local_molecules ...
     238             : !> \param molecule_set ...
     239             : !> \param molecule_kind_set ...
     240             : !> \param csvr ...
     241             : !> \param para_env ...
     242             : !> \param gci ...
     243             : !> \author Teodoro Laino [tlaino] - University of Zurich - 10.2007
     244             : ! **************************************************************************************************
     245           6 :    SUBROUTINE csvr_to_shell_mapping(thermostat_info, simpar, local_molecules, &
     246             :                                     molecule_set, molecule_kind_set, csvr, para_env, gci)
     247             : 
     248             :       TYPE(thermostat_info_type), POINTER                :: thermostat_info
     249             :       TYPE(simpar_type), POINTER                         :: simpar
     250             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     251             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     252             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     253             :       TYPE(csvr_system_type), POINTER                    :: csvr
     254             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     255             :       TYPE(global_constraint_type), POINTER              :: gci
     256             : 
     257             :       INTEGER                                            :: i, imap, j, nshell_local, &
     258             :                                                             sum_of_thermostats
     259           6 :       INTEGER, DIMENSION(:), POINTER                     :: deg_of_freedom, massive_shell_list
     260             :       TYPE(map_info_type), POINTER                       :: map_info
     261             : 
     262           6 :       NULLIFY (massive_shell_list, deg_of_freedom)
     263             : 
     264           6 :       SELECT CASE (simpar%ensemble)
     265             :       CASE DEFAULT
     266           0 :          CPABORT('Unknown ensemble!')
     267             :       CASE (isokin_ensemble, nph_uniaxial_ensemble, &
     268             :             nph_uniaxial_damped_ensemble, reftraj_ensemble, langevin_ensemble)
     269           0 :          CPABORT('Never reach this point!')
     270             :       CASE (nve_ensemble, npe_f_ensemble, npe_i_ensemble, nvt_ensemble, npt_i_ensemble, npt_f_ensemble, &
     271             :             npt_ia_ensemble)
     272             : 
     273             :          CALL setup_csvr_thermostat(csvr, thermostat_info, deg_of_freedom, massive_shell_list, &
     274             :                                     molecule_kind_set, local_molecules, molecule_set, para_env, nshell_local, &
     275           6 :                                     simpar, sum_of_thermostats, gci, shell=.TRUE.)
     276             : 
     277           6 :          map_info => csvr%map_info
     278             :          ! Sum up the number of degrees of freedom on each thermostat.
     279             :          ! first: initialize the target
     280         200 :          map_info%s_kin = 0.0_dp
     281         294 :          DO j = 1, nshell_local
     282        1158 :             DO i = 1, 3
     283        1152 :                map_info%p_kin(i, j)%point = map_info%p_kin(i, j)%point + 1
     284             :             END DO
     285             :          END DO
     286             : 
     287             :          ! If thermostats are replicated but molecules distributed, we have to
     288             :          ! sum s_kin over all processors
     289          10 :          IF (map_info%dis_type == do_thermo_communication) CALL para_env%sum(map_info%s_kin)
     290             : 
     291             :          ! Now that we know how many there are stick this into csvr%nkt
     292             :          ! (number of degrees of freedom times k_B T )
     293         200 :          DO i = 1, csvr%loc_num_csvr
     294         194 :             imap = map_info%map_index(i)
     295         194 :             csvr%nvt(i)%nkt = simpar%temp_sh_ext*map_info%s_kin(imap)
     296         200 :             csvr%nvt(i)%degrees_of_freedom = FLOOR(map_info%s_kin(imap))
     297             :          END DO
     298             : 
     299           6 :          DEALLOCATE (deg_of_freedom)
     300           6 :          DEALLOCATE (massive_shell_list)
     301             :       END SELECT
     302             : 
     303           6 :    END SUBROUTINE csvr_to_shell_mapping
     304             : 
     305             : END MODULE csvr_system_mapping