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 Noam Bernstein [noamb] 02.2012
      10             : ! **************************************************************************************************
      11             : MODULE al_system_dynamics
      12             : 
      13             :    USE al_system_types,                 ONLY: al_system_type
      14             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      15             :                                               get_atomic_kind
      16             :    USE constraint_fxd,                  ONLY: fix_atom_control
      17             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      18             :    USE extended_system_types,           ONLY: map_info_type
      19             :    USE force_env_types,                 ONLY: force_env_type
      20             :    USE kinds,                           ONLY: dp
      21             :    USE message_passing,                 ONLY: mp_comm_type
      22             :    USE molecule_kind_types,             ONLY: molecule_kind_type
      23             :    USE molecule_types,                  ONLY: get_molecule,&
      24             :                                               molecule_type
      25             :    USE particle_types,                  ONLY: particle_type
      26             :    USE thermostat_utils,                ONLY: ke_region_particles,&
      27             :                                               vel_rescale_particles
      28             : #include "../../base/base_uses.f90"
      29             : 
      30             :    IMPLICIT NONE
      31             : 
      32             :    PRIVATE
      33             :    LOGICAL, PARAMETER :: debug_this_module = .FALSE.
      34             :    PUBLIC :: al_particles
      35             : 
      36             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'al_system_dynamics'
      37             : 
      38             : CONTAINS
      39             : 
      40             : ! **************************************************************************************************
      41             : !> \brief ...
      42             : !> \param al ...
      43             : !> \param force_env ...
      44             : !> \param molecule_kind_set ...
      45             : !> \param molecule_set ...
      46             : !> \param particle_set ...
      47             : !> \param local_molecules ...
      48             : !> \param local_particles ...
      49             : !> \param group ...
      50             : !> \param vel ...
      51             : !> \author Noam Bernstein [noamb] 02.2012
      52             : ! **************************************************************************************************
      53          32 :    SUBROUTINE al_particles(al, force_env, molecule_kind_set, molecule_set, &
      54          16 :                            particle_set, local_molecules, local_particles, group, vel)
      55             : 
      56             :       TYPE(al_system_type), POINTER                      :: al
      57             :       TYPE(force_env_type), POINTER                      :: force_env
      58             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
      59             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
      60             :       TYPE(particle_type), POINTER                       :: particle_set(:)
      61             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
      62             :       TYPE(mp_comm_type), INTENT(IN)                     :: group
      63             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: vel(:, :)
      64             : 
      65             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'al_particles'
      66             : 
      67             :       INTEGER                                            :: handle
      68             :       LOGICAL                                            :: my_shell_adiabatic
      69             :       TYPE(map_info_type), POINTER                       :: map_info
      70             : 
      71          16 :       CALL timeset(routineN, handle)
      72          16 :       my_shell_adiabatic = .FALSE.
      73          16 :       map_info => al%map_info
      74             : 
      75             :       IF (debug_this_module) &
      76             :          CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "INIT")
      77             : 
      78          16 :       IF (al%tau_nh <= 0.0_dp) THEN
      79             :          CALL al_OU_step(0.5_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
      80           0 :                          particle_set, local_molecules, local_particles, vel)
      81             :          IF (debug_this_module) &
      82             :             CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post OU")
      83             :       ELSE
      84             :          ! quarter step of Langevin using Ornstein-Uhlenbeck
      85             :          CALL al_OU_step(0.25_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
      86          24 :                          particle_set, local_molecules, local_particles, vel)
      87             :          IF (debug_this_module) &
      88             :             CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post 1st OU")
      89             : 
      90             :          ! Compute the kinetic energy for the region to thermostat for the (T dependent chi step)
      91             :          CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
      92          24 :                                   local_molecules, molecule_set, group, vel=vel)
      93             :          ! quarter step of chi, and set vel drag factors for a half step
      94          16 :          CALL al_NH_quarter_step(al, map_info, set_half_step_vel_factors=.TRUE.)
      95             : 
      96             :          ! Now scale the particle velocities for a NH half step
      97             :          CALL vel_rescale_particles(map_info, molecule_kind_set, molecule_set, particle_set, &
      98          24 :                                     local_molecules, my_shell_adiabatic, vel=vel)
      99             :          ! Recompute the kinetic energy for the region to thermostat (for the T dependent chi step)
     100             :          CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
     101          24 :                                   local_molecules, molecule_set, group, vel=vel)
     102             :          IF (debug_this_module) &
     103             :             CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post rescale_vel")
     104             : 
     105             :          ! quarter step of chi
     106          16 :          CALL al_NH_quarter_step(al, map_info, set_half_step_vel_factors=.FALSE.)
     107             : 
     108             :          ! quarter step of Langevin using Ornstein-Uhlenbeck
     109             :          CALL al_OU_step(0.25_dp, al, force_env, map_info, molecule_kind_set, molecule_set, &
     110          24 :                          particle_set, local_molecules, local_particles, vel)
     111             :          IF (debug_this_module) &
     112             :             CALL dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, "post 2nd OU")
     113             :       END IF
     114             : 
     115             :       ! Recompute the final kinetic energy for the region to thermostat
     116             :       CALL ke_region_particles(map_info, particle_set, molecule_kind_set, &
     117          24 :                                local_molecules, molecule_set, group, vel=vel)
     118             : 
     119          16 :       CALL timestop(handle)
     120          16 :    END SUBROUTINE al_particles
     121             : 
     122             : ! **************************************************************************************************
     123             : !> \brief ...
     124             : !> \param molecule_kind_set ...
     125             : !> \param molecule_set ...
     126             : !> \param local_molecules ...
     127             : !> \param particle_set ...
     128             : !> \param vel ...
     129             : !> \param label ...
     130             : ! **************************************************************************************************
     131           0 :    SUBROUTINE dump_vel(molecule_kind_set, molecule_set, local_molecules, particle_set, vel, label)
     132             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     133             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     134             :       TYPE(distribution_1d_type), POINTER                :: local_molecules
     135             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     136             :       REAL(dp), OPTIONAL                                 :: vel(:, :)
     137             :       CHARACTER(len=*)                                   :: label
     138             : 
     139             :       INTEGER                                            :: first_atom, ikind, imol, imol_local, &
     140             :                                                             ipart, last_atom, nmol_local
     141             :       TYPE(molecule_type), POINTER                       :: molecule
     142             : 
     143           0 :       DO ikind = 1, SIZE(molecule_kind_set)
     144           0 :          nmol_local = local_molecules%n_el(ikind)
     145           0 :          DO imol_local = 1, nmol_local
     146           0 :             imol = local_molecules%list(ikind)%array(imol_local)
     147           0 :             molecule => molecule_set(imol)
     148           0 :             CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
     149           0 :             DO ipart = first_atom, last_atom
     150           0 :                IF (PRESENT(vel)) THEN
     151           0 :                   WRITE (unit=*, fmt='("VEL ",A20," IPART ",I6," V ",3F20.10)') TRIM(label), ipart, vel(:, ipart)
     152             :                ELSE
     153           0 :                   WRITE (unit=*, fmt='("PARTICLE_SET%VEL ",A20," IPART ",I6," V ",3F20.10)') TRIM(label), &
     154           0 :                      ipart, particle_set(ipart)%v(:)
     155             :                END IF
     156             :             END DO
     157             :          END DO
     158             :       END DO
     159           0 :    END SUBROUTINE dump_vel
     160             : 
     161             : ! **************************************************************************************************
     162             : !> \brief ...
     163             : !> \param step ...
     164             : !> \param al ...
     165             : !> \param force_env ...
     166             : !> \param map_info ...
     167             : !> \param molecule_kind_set ...
     168             : !> \param molecule_set ...
     169             : !> \param particle_set ...
     170             : !> \param local_molecules ...
     171             : !> \param local_particles ...
     172             : !> \param vel ...
     173             : ! **************************************************************************************************
     174          32 :    SUBROUTINE al_OU_step(step, al, force_env, map_info, molecule_kind_set, molecule_set, &
     175          32 :                          particle_set, local_molecules, local_particles, vel)
     176             :       REAL(dp), INTENT(in)                               :: step
     177             :       TYPE(al_system_type), POINTER                      :: al
     178             :       TYPE(force_env_type), POINTER                      :: force_env
     179             :       TYPE(map_info_type), POINTER                       :: map_info
     180             :       TYPE(molecule_kind_type), POINTER                  :: molecule_kind_set(:)
     181             :       TYPE(molecule_type), POINTER                       :: molecule_set(:)
     182             :       TYPE(particle_type), POINTER                       :: particle_set(:)
     183             :       TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
     184             :       REAL(KIND=dp), INTENT(INOUT), OPTIONAL             :: vel(:, :)
     185             : 
     186             :       INTEGER :: first_atom, i, ii, ikind, imap, imol, imol_local, ipart, iparticle_kind, &
     187             :          iparticle_local, jj, last_atom, nmol_local, nparticle, nparticle_kind, nparticle_local
     188             :       LOGICAL                                            :: check, present_vel
     189             :       REAL(KIND=dp)                                      :: mass
     190          32 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: w(:, :)
     191             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     192             :       TYPE(molecule_type), POINTER                       :: molecule
     193             : 
     194          32 :       present_vel = PRESENT(vel)
     195             : 
     196             :       ![NB] not a big deal, but could this be done once at init time?
     197      396776 :       DO i = 1, al%loc_num_al
     198      396744 :          imap = map_info%map_index(i)
     199             :          ! drag on velocities
     200      396776 :          IF (al%tau_langevin > 0.0_dp) THEN
     201      396744 :             map_info%v_scale(imap) = EXP(-step*al%dt/al%tau_langevin)
     202      396744 :             map_info%s_kin(imap) = SQRT((al%nvt(i)%nkt/al%nvt(i)%degrees_of_freedom)*(1.0_dp - map_info%v_scale(imap)**2))
     203             :          ELSE
     204           0 :             map_info%v_scale(imap) = 1.0_dp
     205           0 :             map_info%s_kin(imap) = 0.0_dp
     206             :          END IF
     207             :          ! magnitude of random force, not including 1/sqrt(mass) part
     208             :       END DO
     209             : 
     210          32 :       nparticle = SIZE(particle_set)
     211          32 :       nparticle_kind = SIZE(local_particles%n_el)
     212          96 :       ALLOCATE (w(3, nparticle))
     213     1058016 :       w(:, :) = 0.0_dp
     214          32 :       check = (nparticle_kind <= SIZE(local_particles%n_el) .AND. nparticle_kind <= SIZE(local_particles%list))
     215           0 :       CPASSERT(check)
     216          32 :       check = ASSOCIATED(local_particles%local_particle_set)
     217          32 :       CPASSERT(check)
     218        5152 :       DO iparticle_kind = 1, nparticle_kind
     219        5120 :          nparticle_local = local_particles%n_el(iparticle_kind)
     220        5120 :          check = (nparticle_local <= SIZE(local_particles%list(iparticle_kind)%array))
     221        5120 :          CPASSERT(check)
     222      137400 :          DO iparticle_local = 1, nparticle_local
     223      132248 :             ipart = local_particles%list(iparticle_kind)%array(iparticle_local)
     224      132248 :             w(1, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
     225      132248 :             w(2, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
     226      137368 :             w(3, ipart) = local_particles%local_particle_set(iparticle_kind)%rng(iparticle_local)%stream%next(variance=1.0_dp)
     227             :          END DO
     228             :       END DO
     229             : 
     230          32 :       CALL fix_atom_control(force_env, w)
     231             : 
     232          32 :       ii = 0
     233       61320 :       DO ikind = 1, SIZE(molecule_kind_set)
     234       61288 :          nmol_local = local_molecules%n_el(ikind)
     235      101200 :          DO imol_local = 1, nmol_local
     236       39880 :             imol = local_molecules%list(ikind)%array(imol_local)
     237       39880 :             molecule => molecule_set(imol)
     238       39880 :             CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
     239      233416 :             DO ipart = first_atom, last_atom
     240      132248 :                ii = ii + 1
     241      132248 :                atomic_kind => particle_set(ipart)%atomic_kind
     242      132248 :                CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     243      172128 :                IF (present_vel) THEN
     244      264496 :                   DO jj = 1, 3
     245             :                      vel(jj, ipart) = vel(jj, ipart)*map_info%p_scale(jj, ii)%point + &
     246      264496 :                                       map_info%p_kin(jj, ii)%point/SQRT(mass)*w(jj, ipart)
     247             :                   END DO
     248             :                ELSE
     249      264496 :                   DO jj = 1, 3
     250             :                      particle_set(ipart)%v(jj) = particle_set(ipart)%v(jj)*map_info%p_scale(jj, ii)%point + &
     251      264496 :                                                  map_info%p_kin(jj, ii)%point/SQRT(mass)*w(jj, ipart)
     252             :                   END DO
     253             :                END IF
     254             :             END DO
     255             :          END DO
     256             :       END DO
     257             : 
     258          32 :       DEALLOCATE (w)
     259             : 
     260          32 :    END SUBROUTINE al_OU_step
     261             : 
     262             : ! **************************************************************************************************
     263             : !> \brief ...
     264             : !> \param al ...
     265             : !> \param map_info ...
     266             : !> \param set_half_step_vel_factors ...
     267             : !> \author Noam Bernstein [noamb] 02.2012
     268             : ! **************************************************************************************************
     269          32 :    SUBROUTINE al_NH_quarter_step(al, map_info, set_half_step_vel_factors)
     270             :       TYPE(al_system_type), POINTER                      :: al
     271             :       TYPE(map_info_type), POINTER                       :: map_info
     272             :       LOGICAL, INTENT(in)                                :: set_half_step_vel_factors
     273             : 
     274             :       INTEGER                                            :: i, imap
     275             :       REAL(KIND=dp)                                      :: decay, delta_K
     276             : 
     277             : ![NB] how to deal with dt_fact?
     278             : 
     279      396776 :       DO i = 1, al%loc_num_al
     280      396776 :          IF (al%nvt(i)%mass > 0.0_dp) THEN
     281      396744 :             imap = map_info%map_index(i)
     282      396744 :             delta_K = 0.5_dp*(map_info%s_kin(imap) - al%nvt(i)%nkt)
     283      396744 :             al%nvt(i)%chi = al%nvt(i)%chi + 0.5_dp*al%dt*delta_K/al%nvt(i)%mass
     284      396744 :             IF (set_half_step_vel_factors) THEN
     285      198372 :                decay = EXP(-0.5_dp*al%dt*al%nvt(i)%chi)
     286      198372 :                map_info%v_scale(imap) = decay
     287             :             END IF
     288             :          ELSE
     289           0 :             al%nvt(i)%chi = 0.0_dp
     290           0 :             IF (set_half_step_vel_factors) THEN
     291           0 :                map_info%v_scale(imap) = 1.0_dp
     292             :             END IF
     293             :          END IF
     294             :       END DO
     295             : 
     296          32 :    END SUBROUTINE al_NH_quarter_step
     297             : 
     298             : END MODULE al_system_dynamics