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