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             : !> \par History
      10             : !>      CJM 20-Feb-2001: Now npt_ifo is allocated to zero when not used
      11             : !>      CJM 11-apr-2001: adding routines to thermostat ao_type
      12             : !>      CJM 02-Aug-2003: renamed
      13             : !> \author CJM
      14             : ! **************************************************************************************************
      15             : MODULE extended_system_dynamics
      16             : 
      17             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      18             :                                               get_atomic_kind
      19             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      20             :    USE extended_system_types,           ONLY: lnhc_parameters_type,&
      21             :                                               map_info_type,&
      22             :                                               npt_info_type
      23             :    USE kinds,                           ONLY: dp
      24             :    USE message_passing,                 ONLY: mp_comm_type
      25             :    USE molecule_kind_types,             ONLY: molecule_kind_type
      26             :    USE molecule_types,                  ONLY: molecule_type
      27             :    USE particle_types,                  ONLY: particle_type
      28             :    USE shell_potential_types,           ONLY: shell_kind_type
      29             :    USE thermostat_utils,                ONLY: ke_region_baro,&
      30             :                                               ke_region_particles,&
      31             :                                               ke_region_shells,&
      32             :                                               vel_rescale_baro,&
      33             :                                               vel_rescale_particles,&
      34             :                                               vel_rescale_shells
      35             : #include "../../base/base_uses.f90"
      36             : 
      37             :    IMPLICIT NONE
      38             : 
      39             :    PRIVATE
      40             :    LOGICAL, PARAMETER :: debug_this_module = .FALSE.
      41             :    PUBLIC :: shell_scale_comv, &
      42             :              lnhc_particles, &
      43             :              lnhc_barostat, &
      44             :              lnhc_shells
      45             : 
      46             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'extended_system_dynamics'
      47             : 
      48             : CONTAINS
      49             : 
      50             : ! **************************************************************************************************
      51             : !> \brief ...
      52             : !> \param nhc ...
      53             : !> \param npt ...
      54             : !> \param group ...
      55             : !> \date 13-DEC-2000
      56             : !> \par History
      57             : !>      none
      58             : !> \author CJM
      59             : ! **************************************************************************************************
      60        3276 :    SUBROUTINE lnhc_barostat(nhc, npt, group)
      61             : 
      62             :       TYPE(lnhc_parameters_type), POINTER                :: nhc
      63             :       TYPE(npt_info_type), DIMENSION(:, :), &
      64             :          INTENT(INOUT)                                   :: npt
      65             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
      66             : 
      67             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'lnhc_barostat'
      68             : 
      69             :       INTEGER                                            :: handle
      70             :       TYPE(map_info_type), POINTER                       :: map_info
      71             : 
      72        3276 :       CALL timeset(routineN, handle)
      73        3276 :       map_info => nhc%map_info
      74             : 
      75             :       ! Compute the kinetic energy of the barostat
      76        3276 :       CALL ke_region_baro(map_info, npt, group)
      77             : 
      78             :       ! Calculate forces on the Nose-Hoover Thermostat and apply chains
      79        3276 :       CALL do_nhc(nhc, map_info)
      80             : 
      81             :       ! Now scale the particle velocities
      82        3276 :       CALL vel_rescale_baro(map_info, npt)
      83             : 
      84        3276 :       CALL timestop(handle)
      85        3276 :    END SUBROUTINE lnhc_barostat
      86             : 
      87             : ! **************************************************************************************************
      88             : !> \brief ...
      89             : !> \param nhc ...
      90             : !> \param molecule_kind_set ...
      91             : !> \param molecule_set ...
      92             : !> \param particle_set ...
      93             : !> \param local_molecules ...
      94             : !> \param group ...
      95             : !> \param shell_adiabatic ...
      96             : !> \param shell_particle_set ...
      97             : !> \param core_particle_set ...
      98             : !> \param vel ...
      99             : !> \param shell_vel ...
     100             : !> \param core_vel ...
     101             : !> \date 14-NOV-2000
     102             : !> \par History
     103             : !>      none
     104             : ! **************************************************************************************************
     105       26488 :    SUBROUTINE lnhc_particles(nhc, molecule_kind_set, molecule_set, &
     106             :                              particle_set, local_molecules, group, shell_adiabatic, &
     107       13244 :                              shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
     108             : 
     109             :       TYPE(lnhc_parameters_type), POINTER                :: nhc
     110             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     111             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     112             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     113             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     114             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
     115             :       LOGICAL, INTENT(IN), OPTIONAL                      :: shell_adiabatic
     116             :       TYPE(particle_type), OPTIONAL, POINTER             :: shell_particle_set(:), &
     117             :                                                             core_particle_set(:)
     118             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: vel(:, :), shell_vel(:, :), &
     119             :                                                             core_vel(:, :)
     120             : 
     121             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'lnhc_particles'
     122             : 
     123             :       INTEGER                                            :: handle
     124             :       LOGICAL                                            :: my_shell_adiabatic
     125             :       TYPE(map_info_type), POINTER                       :: map_info
     126             : 
     127       13244 :       CALL timeset(routineN, handle)
     128       13244 :       my_shell_adiabatic = .FALSE.
     129       13244 :       IF (PRESENT(shell_adiabatic)) my_shell_adiabatic = shell_adiabatic
     130       13244 :       map_info => nhc%map_info
     131             : 
     132             :       ! Compute the kinetic energy for the region to thermostat
     133             :       CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
     134       19866 :                                local_molecules, molecule_set, group, vel)
     135             : 
     136             :       ! Calculate forces on the Nose-Hoover Thermostat and apply chains
     137       13244 :       CALL do_nhc(nhc, map_info)
     138             : 
     139             :       ! Now scale the particle velocities
     140             :       CALL vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, &
     141             :                                  local_molecules, my_shell_adiabatic, shell_particle_set, core_particle_set, &
     142       67522 :                                  vel, shell_vel, core_vel)
     143             : 
     144       13244 :       CALL timestop(handle)
     145       13244 :    END SUBROUTINE lnhc_particles
     146             : 
     147             : ! **************************************************************************************************
     148             : !> \brief ...
     149             : !> \param nhc ...
     150             : !> \param atomic_kind_set ...
     151             : !> \param particle_set ...
     152             : !> \param local_particles ...
     153             : !> \param group ...
     154             : !> \param shell_particle_set ...
     155             : !> \param core_particle_set ...
     156             : !> \param vel ...
     157             : !> \param shell_vel ...
     158             : !> \param core_vel ...
     159             : !> \date 14-NOV-2000
     160             : !> \par History
     161             : !>      none
     162             : ! **************************************************************************************************
     163        2720 :    SUBROUTINE lnhc_shells(nhc, atomic_kind_set, particle_set, local_particles, &
     164        1360 :                           group, shell_particle_set, core_particle_set, vel, shell_vel, core_vel)
     165             : 
     166             :       TYPE(lnhc_parameters_type), POINTER                :: nhc
     167             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
     168             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     169             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     170             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
     171             :       TYPE(particle_type), OPTIONAL, POINTER             :: shell_particle_set(:), &
     172             :                                                             core_particle_set(:)
     173             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: vel(:, :), shell_vel(:, :), &
     174             :                                                             core_vel(:, :)
     175             : 
     176             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'lnhc_shells'
     177             : 
     178             :       INTEGER                                            :: handle
     179             :       TYPE(map_info_type), POINTER                       :: map_info
     180             : 
     181        1360 :       CALL timeset(routineN, handle)
     182        1360 :       map_info => nhc%map_info
     183             : 
     184             :       ! Compute the kinetic energy of the region to thermostat
     185             :       CALL ke_region_shells(map_info, particle_set, atomic_kind_set, local_particles, &
     186        4080 :                             group, core_particle_set, shell_particle_set, core_vel, shell_vel)
     187             : 
     188             :       ! Calculate forces on the Nose-Hoover Thermostat and apply chains
     189        1360 :       CALL do_nhc(nhc, map_info)
     190             : 
     191             :       ! Now scale the particle velocities
     192             :       CALL vel_rescale_shells(map_info, atomic_kind_set, particle_set, local_particles, &
     193        3400 :                               shell_particle_set, core_particle_set, shell_vel, core_vel, vel)
     194             : 
     195        1360 :       CALL timestop(handle)
     196        1360 :    END SUBROUTINE lnhc_shells
     197             : 
     198             : ! **************************************************************************************************
     199             : !> \brief ...
     200             : !> \param nhc ...
     201             : !> \param map_info ...
     202             : !> \author 10.2007 [tlaino] - Teodoro Laino - University of Zurich
     203             : ! **************************************************************************************************
     204       17880 :    SUBROUTINE do_nhc(nhc, map_info)
     205             :       TYPE(lnhc_parameters_type), POINTER                :: nhc
     206             :       TYPE(map_info_type), POINTER                       :: map_info
     207             : 
     208             :       INTEGER                                            :: imap, n
     209             : 
     210             : ! Force on the first bead in every thermostat chain
     211             : 
     212      619284 :       DO n = 1, nhc%loc_num_nhc
     213      601404 :          imap = nhc%map_info%map_index(n)
     214      601404 :          IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
     215      619284 :          nhc%nvt(1, n)%f = (map_info%s_kin(imap) - nhc%nvt(1, n)%nkt)/nhc%nvt(1, n)%mass
     216             :       END DO
     217             : 
     218             :       ! Perform multiple time stepping using Yoshida
     219       17880 :       CALL multiple_step_yoshida(nhc)
     220             : 
     221       17880 :    END SUBROUTINE do_nhc
     222             : 
     223             : ! **************************************************************************************************
     224             : !> \brief ...
     225             : !> \param atomic_kind_set ...
     226             : !> \param local_particles ...
     227             : !> \param particle_set ...
     228             : !> \param com_vel ...
     229             : !> \param shell_vel ...
     230             : !> \param core_vel ...
     231             : !> \date 14-NOV-2000
     232             : !> \par History
     233             : !>      none
     234             : ! **************************************************************************************************
     235          80 :    SUBROUTINE shell_scale_comv(atomic_kind_set, local_particles, particle_set, &
     236          80 :                                com_vel, shell_vel, core_vel)
     237             : 
     238             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind_set(:)
     239             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     240             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     241             :       REAL(KIND=dp), INTENT(IN)                          :: com_vel(:, :)
     242             :       REAL(KIND=dp), INTENT(INOUT)                       :: shell_vel(:, :), core_vel(:, :)
     243             : 
     244             :       INTEGER                                            :: iparticle, iparticle_kind, &
     245             :                                                             iparticle_local, nparticle_kind, &
     246             :                                                             nparticle_local, shell_index
     247             :       LOGICAL                                            :: is_shell
     248             :       REAL(KIND=dp)                                      :: fac_massc, fac_masss, mass, vc(3), vs(3)
     249             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     250             :       TYPE(shell_kind_type), POINTER                     :: shell
     251             : 
     252          80 :       nparticle_kind = SIZE(atomic_kind_set)
     253             : 
     254         240 :       DO iparticle_kind = 1, nparticle_kind
     255         160 :          atomic_kind => atomic_kind_set(iparticle_kind)
     256             :          CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass, &
     257         160 :                               shell_active=is_shell, shell=shell)
     258         240 :          IF (is_shell) THEN
     259         160 :             fac_masss = shell%mass_shell/mass
     260         160 :             fac_massc = shell%mass_core/mass
     261         160 :             nparticle_local = local_particles%n_el(iparticle_kind)
     262        4000 :             DO iparticle_local = 1, nparticle_local
     263        3840 :                iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
     264        3840 :                shell_index = particle_set(iparticle)%shell_index
     265       15360 :                vs(1:3) = shell_vel(1:3, shell_index)
     266       15360 :                vc(1:3) = core_vel(1:3, shell_index)
     267        3840 :                shell_vel(1, shell_index) = com_vel(1, iparticle) + fac_massc*(vs(1) - vc(1))
     268        3840 :                shell_vel(2, shell_index) = com_vel(2, iparticle) + fac_massc*(vs(2) - vc(2))
     269        3840 :                shell_vel(3, shell_index) = com_vel(3, iparticle) + fac_massc*(vs(3) - vc(3))
     270        3840 :                core_vel(1, shell_index) = com_vel(1, iparticle) + fac_masss*(vc(1) - vs(1))
     271        3840 :                core_vel(2, shell_index) = com_vel(2, iparticle) + fac_masss*(vc(2) - vs(2))
     272        4000 :                core_vel(3, shell_index) = com_vel(3, iparticle) + fac_masss*(vc(3) - vs(3))
     273             :             END DO
     274             :          END IF ! is_shell
     275             :       END DO ! iparticle_kind
     276          80 :    END SUBROUTINE shell_scale_comv
     277             : 
     278             : ! **************************************************************************************************
     279             : !> \brief ...
     280             : !> \param nhc ...
     281             : !> \date 14-NOV-2000
     282             : !> \par History
     283             : !>      none
     284             : ! **************************************************************************************************
     285       17880 :    SUBROUTINE multiple_step_yoshida(nhc)
     286             : 
     287             :       TYPE(lnhc_parameters_type), POINTER                :: nhc
     288             : 
     289             :       INTEGER                                            :: imap, inc, inhc, iyosh, n, nx1, nx2
     290             :       REAL(KIND=dp)                                      :: scale
     291             :       TYPE(map_info_type), POINTER                       :: map_info
     292             : 
     293       17880 :       nx1 = SIZE(nhc%nvt, 1)
     294       17880 :       nx2 = SIZE(nhc%nvt, 2)
     295       17880 :       map_info => nhc%map_info
     296             :       ! perform multiple time stepping using Yoshida
     297       53640 :       NCLOOP: DO inc = 1, nhc%nc
     298      160920 :          YOSH: DO iyosh = 1, nhc%nyosh
     299             : 
     300             :             ! update velocity on the last thermostat in the chain    ! O1
     301             :             nhc%nvt(nhc%nhc_len, :)%v = nhc%nvt(nhc%nhc_len, :)%v + &
     302     3715704 :                                         nhc%nvt(nhc%nhc_len, :)%f*0.25_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact
     303             : 
     304             :             ! update velocity of other thermostats on chain (from nhc_len-1 to 1)  ! O2
     305     3715704 :             DO n = 1, nhc%loc_num_nhc
     306     3608424 :                IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
     307    13484928 :                DO inhc = nhc%nhc_len - 1, 1, -1
     308     9813384 :                   scale = EXP(-0.125_dp*nhc%nvt(inhc + 1, n)%v*nhc%dt_yosh(iyosh)*nhc%dt_fact)
     309     9813384 :                   nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v*scale ! scale
     310             :                   nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v + &
     311     9813384 :                                        nhc%nvt(inhc, n)%f*0.25_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact ! shift
     312    13421808 :                   nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v*scale ! scale
     313             :                END DO
     314             :             END DO
     315             : 
     316             :             ! the core of the operator ----- START------
     317             :             ! update nhc positions
     318             :             nhc%nvt(:, :)%eta = nhc%nvt(:, :)%eta + &
     319    17225832 :                                 0.5_dp*nhc%nvt(:, :)%v*nhc%dt_yosh(iyosh)*nhc%dt_fact
     320             : 
     321             :             ! now accumulate the scale factor for particle velocities
     322     3715704 :             DO n = 1, nhc%loc_num_nhc
     323     3608424 :                imap = nhc%map_info%map_index(n)
     324     3608424 :                IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
     325     3715704 :                map_info%v_scale(imap) = map_info%v_scale(imap)*EXP(-0.5_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact*nhc%nvt(1, n)%v)
     326             :             END DO
     327             :             ! the core of the operator ------ END ------
     328             : 
     329             :             ! update the force on first thermostat again (since particle velocities changed)
     330     3715704 :             DO n = 1, nhc%loc_num_nhc
     331     3608424 :                imap = nhc%map_info%map_index(n)
     332     3608424 :                IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
     333             :                nhc%nvt(1, n)%f = (map_info%s_kin(imap)*map_info%v_scale(imap)* &
     334     3715704 :                                   map_info%v_scale(imap) - nhc%nvt(1, n)%nkt)/nhc%nvt(1, n)%mass
     335             :             END DO
     336             : 
     337             :             ! update velocity of other thermostats on chain (from 1 to nhc_len-1)  ! O2
     338      372696 :             DO inhc = 1, nhc%nhc_len - 1
     339    10167120 :                DO n = 1, nhc%loc_num_nhc
     340     9901704 :                   IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
     341     9813384 :                   scale = EXP(-0.125_dp*nhc%nvt(inhc + 1, n)%v*nhc%dt_yosh(iyosh)*nhc%dt_fact)
     342     9813384 :                   nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v*scale ! scale
     343             :                   nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v + &
     344     9813384 :                                        nhc%nvt(inhc, n)%f*0.25_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact ! shift
     345    10167120 :                   nhc%nvt(inhc, n)%v = nhc%nvt(inhc, n)%v*scale ! scale
     346             :                END DO
     347             : 
     348             :                ! updating the forces on all the thermostats
     349    10274400 :                DO n = 1, nhc%loc_num_nhc
     350     9901704 :                   IF (nhc%nvt(1, n)%nkt == 0.0_dp) CYCLE
     351             :                   nhc%nvt(inhc + 1, n)%f = (nhc%nvt(inhc, n)%mass*nhc%nvt(inhc, n)%v &
     352    10167120 :                                             *nhc%nvt(inhc, n)%v - nhc%nvt(inhc + 1, n)%nkt)/nhc%nvt(inhc + 1, n)%mass
     353             :                END DO
     354             :             END DO
     355             :             ! update velocity on last thermostat                             ! O1
     356             :             nhc%nvt(nhc%nhc_len, :)%v = nhc%nvt(nhc%nhc_len, :)%v + &
     357     3751464 :                                         nhc%nvt(nhc%nhc_len, :)%f*0.25_dp*nhc%dt_yosh(iyosh)*nhc%dt_fact
     358             :          END DO YOSH
     359             :       END DO NCLOOP
     360       17880 :    END SUBROUTINE multiple_step_yoshida
     361             : 
     362             : END MODULE extended_system_dynamics