| 1 |
! MTK (Martina Tuckerman and Klein) Nosie' Hoover chains |
| 2 |
! See Molecular Physics, 1996, Vol. 87, No. 5, 1117-1157 for details of |
| 3 |
! the algorithm. |
| 4 |
! |
| 5 |
! author: Charles F. Vardeman II |
| 6 |
! date: 6-12-02 |
| 7 |
! |
| 8 |
! Interfaces: |
| 9 |
! Called from dynamics module: |
| 10 |
! public: nose_hoover_dynamics(): controls nose-hoover chains dynamics |
| 11 |
! for the length of the simulation time: |
| 12 |
! No Arguments: |
| 13 |
|
| 14 |
! private: |
| 15 |
! nose_hoovera: updates nose_hoover chains and then advances particle positions and velocities |
| 16 |
! by velocity verlet. |
| 17 |
! nose_hooverb: updates velocities, second half of velocity verlet. and thermostat velocities and |
| 18 |
! positions. |
| 19 |
! init_hoover: initalize hoover structures |
| 20 |
! apply_NHC_par: performs details of nose-hoover chains. |
| 21 |
|
| 22 |
|
| 23 |
module dynamics_nose_hoover |
| 24 |
use definitions, ONLY : DP,NDIM |
| 25 |
use constants |
| 26 |
use simulation |
| 27 |
use parameter |
| 28 |
use force_module, ONLY : calc_forces |
| 29 |
use force_utilities, ONLY : check |
| 30 |
use velocity, ONLY : scale_velocities |
| 31 |
use dynamics_utilities, ONLY : evolve_time, sim_status, init_dynamics_loop, & |
| 32 |
DT2,DTSQRT,DTSQ2,init_status |
| 33 |
#ifdef MPI |
| 34 |
use mpi_module, ONLY : mpi_allreduce,mpi_comm_world,mpi_err,& |
| 35 |
mpi_double_precision, mpi_sum |
| 36 |
#endif |
| 37 |
|
| 38 |
implicit none |
| 39 |
PRIVATE |
| 40 |
SAVE |
| 41 |
|
| 42 |
|
| 43 |
public :: nose_hoover_dynamics |
| 44 |
|
| 45 |
type (sim_status) :: nose_status |
| 46 |
|
| 47 |
|
| 48 |
contains |
| 49 |
|
| 50 |
|
| 51 |
subroutine nose_hoover_dynamics() |
| 52 |
use velocity, ONLY : calc_temp |
| 53 |
|
| 54 |
|
| 55 |
|
| 56 |
logical :: update_nlist |
| 57 |
logical :: do_pot |
| 58 |
logical :: not_done |
| 59 |
logical :: nmflag |
| 60 |
|
| 61 |
|
| 62 |
|
| 63 |
integer :: step, i |
| 64 |
|
| 65 |
! initialize status (ke,pe,time, etc.) |
| 66 |
call init_status(nose_status) |
| 67 |
call init_dynamics_loop( nose_status ) |
| 68 |
|
| 69 |
! Start the main simulation loop. |
| 70 |
|
| 71 |
do_pot = .true. |
| 72 |
not_done = .true. |
| 73 |
nmflag = .false. |
| 74 |
step = 0 |
| 75 |
|
| 76 |
dynamics: do |
| 77 |
if (.not.not_done) exit dynamics |
| 78 |
|
| 79 |
step = step + 1 |
| 80 |
|
| 81 |
|
| 82 |
if (checktemp.and.(mod(step,check_temp_steps).eq.0)) then |
| 83 |
call calc_temp(nose_status%temp) |
| 84 |
if (dabs(nose_status%temp-target_temp).gt.therm_variance) then |
| 85 |
call scale_velocities(target_temp) |
| 86 |
end if |
| 87 |
end if |
| 88 |
|
| 89 |
call nose_hoovera() |
| 90 |
|
| 91 |
call check(update_nlist) |
| 92 |
|
| 93 |
|
| 94 |
if (do_pot) then |
| 95 |
call calc_forces(update_nlist,nmflag,pe = nose_status%pot_e) |
| 96 |
else |
| 97 |
call calc_forces(update_nlist,nmflag) |
| 98 |
endif |
| 99 |
|
| 100 |
|
| 101 |
call nose_hooverb() |
| 102 |
|
| 103 |
call evolve_time(step,nose_status,do_pot,not_done) |
| 104 |
|
| 105 |
|
| 106 |
#ifdef MPI |
| 107 |
call mpi_barrier(mpi_comm_world,mpi_err) |
| 108 |
#endif |
| 109 |
end do dynamics |
| 110 |
|
| 111 |
end subroutine nose_hoover_dynamics |
| 112 |
|
| 113 |
|
| 114 |
! nose_hoovera(): advances velocities from T to T + DT/2 and positions from T to T + DT |
| 115 |
subroutine nose_hoovera() |
| 116 |
use velocity, ONLY : remove_cm_momentum,remove_angular_momentum |
| 117 |
|
| 118 |
|
| 119 |
|
| 120 |
|
| 121 |
integer :: i, dim |
| 122 |
real( kind = DP) :: ke, v2, instatemp, chi |
| 123 |
real( kind = DP) :: kebar |
| 124 |
|
| 125 |
real( kind = DP ), dimension(ndim) :: v_dim_i |
| 126 |
|
| 127 |
real(kind=DP) :: ke_local |
| 128 |
|
| 129 |
! units for time are femtosec (10^-15 sec) |
| 130 |
! units for distance are angstroms (10^-10 m) |
| 131 |
! units for velocity are angstroms femtosec^-1 |
| 132 |
! units for mass are a.m.u. (1.661 * 10^-27 kg) |
| 133 |
! units for force are kcal mol^-1 angstrom^-1 |
| 134 |
! |
| 135 |
! converter will put the final terms into angstroms. |
| 136 |
! or angstrom/femtosecond. |
| 137 |
|
| 138 |
! converter = 1.0d0/2.390664d3 |
| 139 |
|
| 140 |
|
| 141 |
! call nose_hoover chains |
| 142 |
call apply_NHC_par(dt2) |
| 143 |
|
| 144 |
|
| 145 |
! advance velocities from t to t + dt/2 |
| 146 |
do i = 1, nlocal |
| 147 |
v(1:ndim,i) = v(1:ndim,i) + & |
| 148 |
(dt2*f(1:ndim,i)/mass(i))*KCALMOL_TO_AMUA2FS2 |
| 149 |
end do |
| 150 |
|
| 151 |
! advance positions from t to t+dt |
| 152 |
do i = 1, nlocal |
| 153 |
q(1:ndim,i) = q(1:ndim,i) + dt*v(1:ndim,i) |
| 154 |
end do |
| 155 |
|
| 156 |
|
| 157 |
end subroutine nose_hoovera |
| 158 |
|
| 159 |
|
| 160 |
! advances velocities another half step from t+dt/2 to t + dt |
| 161 |
subroutine nose_hooverb() |
| 162 |
use velocity, ONLY : remove_cm_momentum,remove_angular_momentum |
| 163 |
integer i |
| 164 |
|
| 165 |
! units for time are femtosec (10^-15 sec) |
| 166 |
! units for distance are angstroms (10^-10 m) |
| 167 |
! units for velocity are angstroms femtosec^-1 |
| 168 |
! units for mass are a.m.u. (1.661 * 10^-27 kg) |
| 169 |
! units for force are kcal mol^-1 angstrom^-1 |
| 170 |
|
| 171 |
do i = 1, nlocal |
| 172 |
v(1:ndim,i) = v(1:ndim,i) + & |
| 173 |
(dt2*f(1:ndim,i)/mass(i))*KCALMOL_TO_AMUA2FS2 |
| 174 |
end do |
| 175 |
|
| 176 |
! apply nose-hoover thermostat |
| 177 |
call apply_NHC_par(dt2) |
| 178 |
|
| 179 |
! if (sim_type == "cluster") then |
| 180 |
! call remove_cm_momentum() |
| 181 |
! call remove_angular_momentum() |
| 182 |
! endif |
| 183 |
end subroutine nose_hooverb |
| 184 |
|
| 185 |
|
| 186 |
! does the Nose-Hoover part of the integration from t = 0 to t=DT/2 |
| 187 |
! borrowed from Chris Fennel's code of thermostating |
| 188 |
! by way of Hoover, Phys.Rev.A 1985, Vol.31 (3) 1695-1697 |
| 189 |
|
| 190 |
subroutine apply_NHC_par(this_time) |
| 191 |
use velocity, ONLY : calc_temp |
| 192 |
real(kind = DP) :: system_ke, ke_temp |
| 193 |
real(kind = DP) :: system_temp |
| 194 |
real(kind = DP) :: G_NkT |
| 195 |
real(kind = DP) :: Q_mass |
| 196 |
real(kind = DP) :: zetaScale |
| 197 |
real(kind = DP),dimension(ndim) :: vi |
| 198 |
real(kind = DP),intent(out) :: this_time |
| 199 |
|
| 200 |
integer :: i |
| 201 |
integer :: j |
| 202 |
|
| 203 |
call calc_temp(system_temp, system_ke) |
| 204 |
|
| 205 |
G_NkT = natoms*ndim * 8.31451d-7 * target_temp |
| 206 |
ke_temp = system_ke*KCALMOL_TO_AMUA2FS2 |
| 207 |
|
| 208 |
! write(*,*) 'System_temp: ', system_temp |
| 209 |
! write(*,*) 'ke_temp: ',ke_temp |
| 210 |
|
| 211 |
! Q_mass is a function of the omega parameter |
| 212 |
! Q_mass = G_NkT / (omeg * omeg) |
| 213 |
Q_mass = 50000._DP |
| 214 |
|
| 215 |
! advance the zeta term to zeta(t + dt) |
| 216 |
nzeta = nzeta + this_time*((ke_temp*2 - G_NkT)/Q_mass) |
| 217 |
|
| 218 |
! write(*,*) 'zeta: ',zeta |
| 219 |
! write(*,*) 'dt: ',dt |
| 220 |
! write(*,*) 'ke_temp: ',ke_temp |
| 221 |
! write(*,*) 'G_NkT: ',G_NkT |
| 222 |
! write(*,*) 'Q_mass: ',Q_mass |
| 223 |
|
| 224 |
zetaScale = nzeta * this_time |
| 225 |
|
| 226 |
! write(*,*) 'zetaScale: ',zetaScale |
| 227 |
|
| 228 |
! perform thermostating on velocities and angular momenta |
| 229 |
do i = 1, nlocal |
| 230 |
do j = 1, ndim |
| 231 |
vi(j) = v(j,i)*zetaScale |
| 232 |
v(j,i) = v(j,i) - vi(j) |
| 233 |
enddo |
| 234 |
enddo |
| 235 |
|
| 236 |
end subroutine apply_NHC_par |
| 237 |
|
| 238 |
|
| 239 |
|
| 240 |
!!$ ! other code borrowed from Chris Fennell's code of |
| 241 |
!!$ ! G.J. Martyna et al. Mol. Phys., 1996, Vol. 87, No. 5, 1117-1157 |
| 242 |
!!$ |
| 243 |
!!$ subroutine apply_NHC_par() |
| 244 |
!!$ use velocity, ONLY : calc_temp |
| 245 |
!!$ real(kind = DP) :: scale |
| 246 |
!!$ real(kind = DP) :: system_ke |
| 247 |
!!$ real(kind = DP) :: system_temp |
| 248 |
!!$ real(kind = DP) :: G_NkT |
| 249 |
!!$ real(kind = DP) :: omega2 |
| 250 |
!!$ real(kind = DP) :: value1, value2, value3, value4 |
| 251 |
!!$ real(kind = DP) :: wdti2(n_ysval), wdti4(n_ysval), wdti8(n_ysval) |
| 252 |
!!$ real(kind = DP) :: G_val(N_nos), vval(N_nos), xival(N_nos) |
| 253 |
!!$ real(kind = DP) :: Q_mass |
| 254 |
!!$ |
| 255 |
!!$ integer :: i, j |
| 256 |
!!$ |
| 257 |
!!$ scale = 1._DP |
| 258 |
!!$ |
| 259 |
!!$ call calc_temp(system_temp,system_ke) |
| 260 |
!!$ |
| 261 |
!!$ G_NkT = natoms*ndim*K_BOLTZ*my_temp |
| 262 |
!!$ |
| 263 |
!!$ ! we need to get the w numbers for these values |
| 264 |
!!$ ! but i think that chris has calculated the values |
| 265 |
!!$ ! because they are values of constants |
| 266 |
!!$ value1 = 0.67560359598d0 * dt |
| 267 |
!!$ value2 = -0.85120719196d0 * dt |
| 268 |
!!$ value3 = 0.207245385897d0 * dt |
| 269 |
!!$ value4 = -0.328981543589d0 * dt |
| 270 |
!!$ if (n_ysval.eq.3) then |
| 271 |
!!$ wdti2(1) = value1 |
| 272 |
!!$ wdti2(2) = value2 |
| 273 |
!!$ wdti2(3) = value1 |
| 274 |
!!$ wdti4(1) = 0.5_DP * wdti2(1) |
| 275 |
!!$ wdti4(2) = 0.5_DP * wdti2(2) |
| 276 |
!!$ wdti4(3) = wdti4(1) |
| 277 |
!!$ wdti8(1) = 0.5_DP * wdti4(1) |
| 278 |
!!$ wdti8(2) = 0.5_DP * wdti4(2) |
| 279 |
!!$ wdti8(3) = wdti8(1) |
| 280 |
!!$ else |
| 281 |
!!$ write(*,*) 'Error: n_syval must be 3' |
| 282 |
!!$ stop |
| 283 |
!!$ endif |
| 284 |
!!$ |
| 285 |
!!$ omega2 = omega * omega |
| 286 |
!!$ |
| 287 |
!!$ ! set the Q_mass values |
| 288 |
!!$ Q_mass(1) = G_NkT/omega2 |
| 289 |
!!$ do i = 2, N_nos-1 |
| 290 |
!!$ Q_mass(i) = system_temp/omega2 |
| 291 |
!!$ enddo |
| 292 |
!!$ |
| 293 |
!!$ qmass(N_nos) = 2 * omega |
| 294 |
!!$ |
| 295 |
!!$ ! zero the intial xivals and vvals |
| 296 |
!!$ do i = 1, N_nos |
| 297 |
!!$ xival(i) = 0._DP |
| 298 |
!!$ vval(i) = 0._DP |
| 299 |
!!$ enddo |
| 300 |
!!$ |
| 301 |
!!$ ! Start the Nose_Hoover chain stepping |
| 302 |
!!$ ! update the forces |
| 303 |
!!$ G_val(1) = (system_ke - G_NkT)/Q_mass(1) |
| 304 |
!!$ |
| 305 |
!!$ ! start the multiple time stepping |
| 306 |
!!$ do i = 1, n_cval |
| 307 |
!!$ do j = 1, n_ysval |
| 308 |
!!$ |
| 309 |
!!$ ! update the thermostat velocities |
| 310 |
!!$ vval(N_nos) = vval(N_nos) + G_val(N_nos)*wdti4(j)/i |
| 311 |
!!$ do mcount = 1, N_nos-1 |
| 312 |
!!$ aa = exp((-wdti8(j)/i) * vval((N_nos+1)-mcount)) |
| 313 |
!!$ vval(N_nos-mcount) = vval(N_nos-mcount)*aa*aa & |
| 314 |
!!$ + wdti4(j)*G_val(N_nos)*wdti4(j)/(i*i) |
| 315 |
!!$ enddo |
| 316 |
!!$ |
| 317 |
!!$ ! update particle velocities |
| 318 |
!!$ aa = exp((-wdti2(j)/i)*vval(1)) |
| 319 |
!!$ scale = scale*aa |
| 320 |
!!$ |
| 321 |
!!$ ! update the forces |
| 322 |
!!$ G_val(1) = (scale*scale*ke -G_NkT)/Q_mass(1) |
| 323 |
!!$ |
| 324 |
!!$ ! update the thermostat positions |
| 325 |
!!$ do mcount = 1, N_nos-1 |
| 326 |
!!$ xival(mcount) = xival(mcount) + vval(mcount)*(wdti2(j)/i) |
| 327 |
!!$ enddo |
| 328 |
!!$ |
| 329 |
!!$ ! update the thermostat velocities |
| 330 |
!!$ do mcount = 1, N_nos-1 |
| 331 |
!!$ aa = exp((-wdti8(j)/i)*vval(mcount+1)) |
| 332 |
!!$ vval(mcount) = vval(mcount)*aa*aa & |
| 333 |
!!$ + (wdti4(j)/i)*G_val(mcount)*aa |
| 334 |
!!$ G_val(mcount+1) = (Q_mass(mcount)*vval(mcount)*vval(mcount) & |
| 335 |
!!$ - system_temp)/Q_mass(mcount+1) |
| 336 |
!!$ enddo |
| 337 |
!!$ |
| 338 |
!!$ vval(N_nos) = vval(N_nos) + G_val(N_nos)*(wdti4(j)/i) |
| 339 |
!!$ enddo |
| 340 |
!!$ enddo |
| 341 |
!!$ |
| 342 |
!!$ ! scale both linear velocities and angular momenta |
| 343 |
!!$ do i = 1, nlocal |
| 344 |
!!$ v(1:ndim,i) = v(1:ndim,i)*scale |
| 345 |
!!$ !something here for what we use for angular momenta |
| 346 |
!!$ enddo |
| 347 |
!!$ |
| 348 |
!!$ end subroutine apply_NHC_par |
| 349 |
|
| 350 |
|
| 351 |
|
| 352 |
end module dynamics_nose_hoover |
