| 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,init_status |
| 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 |
| 37 |
|
|
| 38 |
|
implicit none |
| 39 |
|
PRIVATE |
| 40 |
+ |
SAVE |
| 41 |
|
|
| 42 |
|
|
| 42 |
– |
|
| 43 |
|
public :: nose_hoover_dynamics |
| 44 |
|
|
| 45 |
|
type (sim_status) :: nose_status |
| 46 |
|
|
| 47 |
< |
|
| 47 |
> |
|
| 48 |
|
contains |
| 49 |
|
|
| 50 |
|
|
| 58 |
|
logical :: not_done |
| 59 |
|
logical :: nmflag |
| 60 |
|
|
| 61 |
+ |
|
| 62 |
+ |
|
| 63 |
|
integer :: step, i |
| 64 |
|
|
| 65 |
|
! initialize status (ke,pe,time, etc.) |
| 139 |
|
|
| 140 |
|
|
| 141 |
|
! call nose_hoover chains |
| 142 |
< |
call apply_NHC_par() |
| 142 |
> |
call apply_NHC_par(dt2) |
| 143 |
> |
|
| 144 |
|
|
| 145 |
|
! advance velocities from t to t + dt/2 |
| 146 |
|
do i = 1, nlocal |
| 158 |
|
|
| 159 |
|
|
| 160 |
|
! advances velocities another half step from t+dt/2 to t + dt |
| 161 |
< |
subroutine nose_hooverb( ) |
| 161 |
> |
subroutine nose_hooverb() |
| 162 |
|
use velocity, ONLY : remove_cm_momentum,remove_angular_momentum |
| 163 |
|
integer i |
| 164 |
|
|
| 174 |
|
end do |
| 175 |
|
|
| 176 |
|
! apply nose-hoover thermostat |
| 177 |
< |
call apply_NHC_par() |
| 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 |
| 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 |
|
|
| 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() |
| 190 |
> |
subroutine apply_NHC_par(this_time) |
| 191 |
|
use velocity, ONLY : calc_temp |
| 192 |
< |
real(kind = DP) :: system_ke |
| 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 |
| 193 |
– |
real(kind = DP) :: zeta |
| 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 |
< |
|
| 199 |
< |
call calc_temp(system_temp,system_ke) |
| 201 |
> |
integer :: j |
| 202 |
|
|
| 203 |
< |
G_NkT = natoms*ndim*K_BOLTZ*target_temp |
| 204 |
< |
! and aren't the system_temp and system_ke similar to chris's kt and ke? |
| 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 / (omega * omega) |
| 212 |
> |
! Q_mass = G_NkT / (omeg * omeg) |
| 213 |
> |
Q_mass = 50000._DP |
| 214 |
|
|
| 215 |
|
! advance the zeta term to zeta(t + dt) |
| 216 |
< |
zeta = zeta + dt*((system_temp*2 - G_NkT)/Q_mass) |
| 209 |
< |
zetaScale = zeta * 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 |
< |
vi(1:ndim) = v(1:ndim,i)*zetaScale |
| 231 |
< |
! similar thing here for our angular momenta |
| 232 |
< |
|
| 233 |
< |
v(1:ndim,i) = v(1:ndim,i) - vi(1:ndim) |
| 217 |
< |
! similar thing here for our angular momenta |
| 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 |
