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 10.2007 [tlaino] - Teodoro Laino - University of Zurich
      10             : ! **************************************************************************************************
      11             : MODULE csvr_system_dynamics
      12             : 
      13             :    USE atomic_kind_types,               ONLY: atomic_kind_type
      14             :    USE csvr_system_types,               ONLY: csvr_system_type
      15             :    USE csvr_system_utils,               ONLY: rescaling_factor
      16             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      17             :    USE extended_system_types,           ONLY: map_info_type,&
      18             :                                               npt_info_type
      19             :    USE kinds,                           ONLY: dp
      20             :    USE message_passing,                 ONLY: mp_comm_type
      21             :    USE molecule_kind_types,             ONLY: molecule_kind_type
      22             :    USE molecule_types,                  ONLY: molecule_type
      23             :    USE parallel_rng_types,              ONLY: rng_stream_type
      24             :    USE particle_types,                  ONLY: particle_type
      25             :    USE thermostat_utils,                ONLY: ke_region_baro,&
      26             :                                               ke_region_particles,&
      27             :                                               ke_region_shells,&
      28             :                                               vel_rescale_baro,&
      29             :                                               vel_rescale_particles,&
      30             :                                               vel_rescale_shells
      31             : #include "../../base/base_uses.f90"
      32             : 
      33             :    IMPLICIT NONE
      34             : 
      35             :    PRIVATE
      36             :    LOGICAL, PARAMETER :: debug_this_module = .FALSE.
      37             :    PUBLIC :: csvr_particles, &
      38             :              csvr_barostat, &
      39             :              csvr_shells
      40             : 
      41             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'csvr_system_dynamics'
      42             : 
      43             : CONTAINS
      44             : 
      45             : ! **************************************************************************************************
      46             : !> \brief ...
      47             : !> \param csvr ...
      48             : !> \param npt ...
      49             : !> \param group ...
      50             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
      51             : ! **************************************************************************************************
      52        1084 :    SUBROUTINE csvr_barostat(csvr, npt, group)
      53             : 
      54             :       TYPE(csvr_system_type), POINTER                    :: csvr
      55             :       TYPE(npt_info_type), DIMENSION(:, :), &
      56             :          INTENT(INOUT)                                   :: npt
      57             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
      58             : 
      59             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'csvr_barostat'
      60             : 
      61             :       INTEGER                                            :: handle
      62             :       TYPE(map_info_type), POINTER                       :: map_info
      63             : 
      64        1084 :       CALL timeset(routineN, handle)
      65        1084 :       map_info => csvr%map_info
      66             : 
      67             :       ! Compute the kinetic energy of the barostat
      68        1084 :       CALL ke_region_baro(map_info, npt, group)
      69             : 
      70             :       ! Apply the Canonical Sampling through Velocity Rescaling
      71        1084 :       CALL do_csvr(csvr, map_info)
      72             : 
      73             :       ! Now scale the particle velocities
      74        1084 :       CALL vel_rescale_baro(map_info, npt)
      75             : 
      76             :       ! Re-Compute the kinetic energy of the barostat
      77        1084 :       CALL ke_region_baro(map_info, npt, group)
      78             : 
      79             :       ! Compute thermostat energy
      80        1084 :       CALL do_csvr_eval_energy(csvr, map_info)
      81             : 
      82        1084 :       CALL timestop(handle)
      83        1084 :    END SUBROUTINE csvr_barostat
      84             : 
      85             : ! **************************************************************************************************
      86             : !> \brief ...
      87             : !> \param csvr ...
      88             : !> \param molecule_kind_set ...
      89             : !> \param molecule_set ...
      90             : !> \param particle_set ...
      91             : !> \param local_molecules ...
      92             : !> \param group ...
      93             : !> \param shell_adiabatic ...
      94             : !> \param shell_particle_set ...
      95             : !> \param core_particle_set ...
      96             : !> \param vel ...
      97             : !> \param shell_vel ...
      98             : !> \param core_vel ...
      99             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
     100             : ! **************************************************************************************************
     101        9992 :    SUBROUTINE csvr_particles(csvr, molecule_kind_set, molecule_set, &
     102             :                              particle_set, local_molecules, group, shell_adiabatic, &
     103        4996 :                              shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
     104             : 
     105             :       TYPE(csvr_system_type), POINTER                    :: csvr
     106             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     107             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     108             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     109             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     110             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
     111             :       LOGICAL, INTENT(IN), OPTIONAL                      :: shell_adiabatic
     112             :       TYPE(particle_type), OPTIONAL, POINTER             :: shell_particle_set(:), &
     113             :                                                             core_particle_set(:)
     114             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: vel(:, :), shell_vel(:, :), &
     115             :                                                             core_vel(:, :)
     116             : 
     117             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'csvr_particles'
     118             : 
     119             :       INTEGER                                            :: handle
     120             :       LOGICAL                                            :: my_shell_adiabatic
     121             :       TYPE(map_info_type), POINTER                       :: map_info
     122             : 
     123        4996 :       CALL timeset(routineN, handle)
     124        4996 :       my_shell_adiabatic = .FALSE.
     125        4996 :       IF (PRESENT(shell_adiabatic)) my_shell_adiabatic = shell_adiabatic
     126        4996 :       map_info => csvr%map_info
     127             : 
     128             :       ! Compute the kinetic energy for the region to thermostat
     129             :       CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
     130        7494 :                                local_molecules, molecule_set, group, vel)
     131             : 
     132             :       ! Apply the Canonical Sampling through Velocity Rescaling
     133        4996 :       CALL do_csvr(csvr, map_info)
     134             : 
     135             :       ! Now scale the particle velocities
     136             :       CALL vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, &
     137             :                                  local_molecules, my_shell_adiabatic, shell_particle_set, core_particle_set, &
     138       27158 :                                  vel, shell_vel, core_vel)
     139             : 
     140             :       ! Re-Compute the kinetic energy for the region to thermostat
     141             :       CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
     142        7494 :                                local_molecules, molecule_set, group, vel)
     143             : 
     144             :       ! Compute thermostat energy
     145        4996 :       CALL do_csvr_eval_energy(csvr, map_info)
     146             : 
     147        4996 :       CALL timestop(handle)
     148        4996 :    END SUBROUTINE csvr_particles
     149             : 
     150             : ! **************************************************************************************************
     151             : !> \brief ...
     152             : !> \param csvr ...
     153             : !> \param atomic_kind_set ...
     154             : !> \param particle_set ...
     155             : !> \param local_particles ...
     156             : !> \param group ...
     157             : !> \param shell_particle_set ...
     158             : !> \param core_particle_set ...
     159             : !> \param vel ...
     160             : !> \param shell_vel ...
     161             : !> \param core_vel ...
     162             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
     163             : ! **************************************************************************************************
     164         480 :    SUBROUTINE csvr_shells(csvr, atomic_kind_set, particle_set, local_particles, &
     165         240 :                           group, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
     166             : 
     167             :       TYPE(csvr_system_type), POINTER                    :: csvr
     168             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
     169             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     170             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     171             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
     172             :       TYPE(particle_type), OPTIONAL, POINTER             :: shell_particle_set(:), &
     173             :                                                             core_particle_set(:)
     174             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: vel(:, :), shell_vel(:, :), &
     175             :                                                             core_vel(:, :)
     176             : 
     177             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'csvr_shells'
     178             : 
     179             :       INTEGER                                            :: handle
     180             :       TYPE(map_info_type), POINTER                       :: map_info
     181             : 
     182         240 :       CALL timeset(routineN, handle)
     183         240 :       map_info => csvr%map_info
     184             : 
     185             :       ! Compute the kinetic energy of the region to thermostat
     186             :       CALL ke_region_shells(map_info, particle_set, atomic_kind_set, local_particles, &
     187         720 :                             group, core_particle_set, shell_particle_set, core_vel, shell_vel)
     188             : 
     189             :       ! Apply the Canonical Sampling through Velocity Rescaling
     190         240 :       CALL do_csvr(csvr, map_info)
     191             : 
     192             :       ! Now scale the particle velocities
     193             :       CALL vel_rescale_shells(map_info, atomic_kind_set, particle_set, local_particles, &
     194         600 :                               shell_particle_set, core_particle_set, shell_vel, core_vel, vel)
     195             : 
     196             :       ! Re-Compute the kinetic energy of the region to thermostat
     197             :       CALL ke_region_shells(map_info, particle_set, atomic_kind_set, local_particles, &
     198         480 :                             group, core_particle_set, shell_particle_set, core_vel, shell_vel)
     199             : 
     200             :       ! Compute thermostat energy
     201         240 :       CALL do_csvr_eval_energy(csvr, map_info)
     202             : 
     203         240 :       CALL timestop(handle)
     204         240 :    END SUBROUTINE csvr_shells
     205             : 
     206             : ! **************************************************************************************************
     207             : !> \brief ...
     208             : !> \param csvr ...
     209             : !> \param map_info ...
     210             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
     211             : ! **************************************************************************************************
     212        6320 :    SUBROUTINE do_csvr(csvr, map_info)
     213             :       TYPE(csvr_system_type), POINTER                    :: csvr
     214             :       TYPE(map_info_type), POINTER                       :: map_info
     215             : 
     216             :       INTEGER                                            :: i, imap, ndeg
     217             :       REAL(KIND=dp)                                      :: kin_energy, kin_target, taut
     218             :       TYPE(rng_stream_type), POINTER                     :: rng_stream
     219             : 
     220      120604 :       DO i = 1, csvr%loc_num_csvr
     221      114284 :          imap = map_info%map_index(i)
     222      114284 :          kin_energy = map_info%s_kin(imap)
     223      114284 :          csvr%nvt(i)%region_kin_energy = kin_energy
     224      114284 :          kin_energy = kin_energy*0.5_dp
     225      114284 :          kin_target = csvr%nvt(i)%nkt*0.5_dp
     226      114284 :          ndeg = csvr%nvt(i)%degrees_of_freedom
     227      114284 :          taut = csvr%tau_csvr/csvr%dt_fact
     228      114284 :          rng_stream => csvr%nvt(i)%gaussian_rng_stream
     229             :          map_info%v_scale(imap) = rescaling_factor(kin_energy, kin_target, ndeg, taut, &
     230      120604 :                                                    rng_stream)
     231             :       END DO
     232             : 
     233        6320 :    END SUBROUTINE do_csvr
     234             : 
     235             : ! **************************************************************************************************
     236             : !> \brief ...
     237             : !> \param csvr ...
     238             : !> \param map_info ...
     239             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
     240             : ! **************************************************************************************************
     241        6320 :    SUBROUTINE do_csvr_eval_energy(csvr, map_info)
     242             :       TYPE(csvr_system_type), POINTER                    :: csvr
     243             :       TYPE(map_info_type), POINTER                       :: map_info
     244             : 
     245             :       INTEGER                                            :: i, imap
     246             :       REAL(KIND=dp)                                      :: kin_energy_ar, kin_energy_br
     247             : 
     248      120604 :       DO i = 1, csvr%loc_num_csvr
     249      114284 :          imap = map_info%map_index(i)
     250      114284 :          kin_energy_br = csvr%nvt(i)%region_kin_energy
     251      114284 :          kin_energy_ar = map_info%s_kin(imap)
     252             :          csvr%nvt(i)%thermostat_energy = csvr%nvt(i)%thermostat_energy + &
     253      120604 :                                          0.5_dp*(kin_energy_br - kin_energy_ar)
     254             :       END DO
     255             : 
     256        6320 :    END SUBROUTINE do_csvr_eval_energy
     257             : 
     258             : END MODULE csvr_system_dynamics