mirror of
https://github.com/nwchemgit/nwchem.git
synced 2026-07-28 14:15:30 -04:00
Nose-Hoover chain with arbitrary number of oscillators
This commit is contained in:
Daniel Mejia-Rodriguez
parent
76060f7e42
commit
d98fff175f
4 changed files with 166 additions and 58 deletions
|
|
@ -7,13 +7,14 @@ c
|
|||
common/qmd_int/lgeom,igeom,jdebug
|
||||
c
|
||||
c nose-hoover thermostat parameters
|
||||
integer maxchain
|
||||
parameter (maxchain = 2)
|
||||
integer maxchain,maxnsy,nnhos,nnhsy
|
||||
parameter (maxchain = 10, maxnsy = 7)
|
||||
logical do_nosehoover
|
||||
double precision r_nh(maxchain) ! chain positions
|
||||
double precision v_nh(maxchain) ! chain velocities
|
||||
double precision v_nh(maxchain) ! chain velocities
|
||||
double precision m_nh(maxchain) ! chain masses
|
||||
double precision g_nh(maxchain) ! bath couplings
|
||||
double precision syw_nh(maxnsy) ! Suzuki-Yoshida weights
|
||||
c
|
||||
c berendsen thermostat
|
||||
logical do_berendsen
|
||||
|
|
@ -31,7 +32,8 @@ c stochastic velocity rescaling
|
|||
logical do_svr
|
||||
c
|
||||
common/qmd_thermostats/do_berendsen,tau,
|
||||
& do_nosehoover,r_nh,v_nh,m_nh,g_nh,
|
||||
& do_nosehoover,r_nh,v_nh,m_nh,g_nh,syw_nh,
|
||||
& nnhos, nnhsy,
|
||||
& do_langevin,friction,
|
||||
& do_rescale,idum,
|
||||
& do_svr
|
||||
|
|
|
|||
|
|
@ -19,6 +19,7 @@ c
|
|||
#include "msgids.fh"
|
||||
#include "stdio.fh"
|
||||
#include "msgtypesf.h"
|
||||
#include "util_params.fh"
|
||||
c
|
||||
#include "qmd_common.fh"
|
||||
c
|
||||
|
|
@ -41,7 +42,7 @@ c
|
|||
character*32 integrator_nucl ! nuclear md integrator
|
||||
character*32 integrator_elec ! electronic dynamics integrator
|
||||
integer ichain
|
||||
double precision nh_mass(2) ! nose-hoover mass
|
||||
double precision nhfreq ! nose-hoover mass
|
||||
logical do_none
|
||||
logical restart
|
||||
logical do_namd
|
||||
|
|
@ -131,19 +132,38 @@ c thermostat
|
|||
end if
|
||||
c
|
||||
c nose-hoover mass parameter
|
||||
if (.not.rtdb_get(rtdb,'qmd:nh_mass',mt_dbl,2,nh_mass))
|
||||
if (.not.rtdb_get(rtdb,'qmd:nhfreq',mt_dbl,1,nhfreq))
|
||||
& call errquit(pname//'failed to read nh_mass',0,RTDB_ERR)
|
||||
if (.not.rtdb_get(rtdb,'qmd:nnhos',mt_int,1,nnhos))
|
||||
& call errquit(pname//'failed to read nh_nchain',0,RTDB_ERR)
|
||||
if (.not.rtdb_get(rtdb,'qmd:nnhsy',mt_int,1,nnhsy))
|
||||
& call errquit(pname//'failed to read nh_nsy',0,RTDB_ERR)
|
||||
if (do_nosehoover) then
|
||||
do ichain = 1,maxchain
|
||||
r_nh(ichain) = 0.d0 ! position of chains
|
||||
v_nh(ichain) = 0.d0 ! velocity of chains
|
||||
m_nh(ichain) = nh_mass(ichain) ! mass of chains
|
||||
g_nh(ichain) = 0.d0 ! bath coupling
|
||||
end do
|
||||
do ichain = 1,nnhos
|
||||
r_nh(ichain) = 0.d0 ! position of chains
|
||||
m_nh(ichain) = kb_au*targ_temp/nhfreq**2 ! mass of chains
|
||||
v_nh(ichain) = 0.d0 ! velocity of chains
|
||||
g_nh(ichain) = 0.d0 ! bath coupling
|
||||
end do
|
||||
|
||||
if (nnhsy.eq.1) then
|
||||
syw_nh(1) = 1d0
|
||||
elseif (nnhsy.eq.3) then
|
||||
syw_nh(1) = 1d0/(2d0 - 2d0**(1d0/3d0))
|
||||
syw_nh(2) = 1 - 2*syw_nh(1)
|
||||
syw_nh(3) = syw_nh(1)
|
||||
elseif (nnhsy.eq.7) then
|
||||
syw_nh(1) = 0.784513610477560d0
|
||||
syw_nh(2) = 0.235573213359357d0
|
||||
syw_nh(3) = -1.17767998417887d0
|
||||
syw_nh(4) = 1d0 - 2*sum(syw_nh(1:3))
|
||||
syw_nh(5) = syw_nh(3)
|
||||
syw_nh(6) = syw_nh(2)
|
||||
syw_nh(7) = syw_nh(1)
|
||||
endif
|
||||
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,130) "NH mass(1):",nh_mass(1)
|
||||
if (ga_nodeid().eq.0)
|
||||
& write(luout,130) "NH mass(2):",nh_mass(2)
|
||||
& write(luout,130) "NH frequency:",nhfreq
|
||||
end if
|
||||
c
|
||||
c langevin friction parameter
|
||||
|
|
|
|||
|
|
@ -31,9 +31,9 @@ c
|
|||
integer tvals(8)
|
||||
double precision targ_temp
|
||||
double precision friction
|
||||
integer idum
|
||||
integer idum, nnhos, nnhsy
|
||||
double precision tau
|
||||
double precision nh_mass(2)
|
||||
double precision nh_freq
|
||||
character*32 test
|
||||
character*32 thermostat
|
||||
character*32 integrator_nucl
|
||||
|
|
@ -83,10 +83,15 @@ c Berendsen tau parameter, also used for svr
|
|||
& call errquit(pname//'failed to write tau',0,RTDB_ERR)
|
||||
c
|
||||
c Nose-Hoover mass parameter
|
||||
nh_mass(1)=0.1d0
|
||||
nh_mass(2)=0.1d0
|
||||
if (.not.rtdb_put(rtdb,'qmd:nh_mass',mt_dbl,2,nh_mass(1)))
|
||||
& call errquit(pname//'failed to write nh_mass',0,RTDB_ERR)
|
||||
nh_freq = 0.013669d0
|
||||
nnhos = 3
|
||||
nnhsy = 3
|
||||
if (.not.rtdb_put(rtdb,'qmd:nhfreq',mt_dbl,1,nh_freq))
|
||||
& call errquit(pname//'failed to write tau',0,RTDB_ERR)
|
||||
if (.not.rtdb_put(rtdb,'qmd:nnhos',mt_dbl,1,nnhos))
|
||||
& call errquit(pname//'failed to write tau',0,RTDB_ERR)
|
||||
if (.not.rtdb_put(rtdb,'qmd:nnhsy',mt_dbl,1,nnhsy))
|
||||
& call errquit(pname//'failed to write tau',0,RTDB_ERR)
|
||||
c
|
||||
c Nuclear MD integrator
|
||||
integrator_nucl='velocity-verlet'
|
||||
|
|
@ -199,10 +204,48 @@ c Nuclear thermostat
|
|||
1 call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
|
||||
c Nose-Hoover: get masses
|
||||
if (thermostat.eq.'nose-hoover') then
|
||||
if (.not.inp_f(nh_mass(1))) nh_mass(1) = 0.1d0
|
||||
if (.not.inp_f(nh_mass(2))) nh_mass(2) = 0.1d0
|
||||
if (.not.rtdb_put(rtdb,'qmd:nh_mass',mt_dbl,2,nh_mass(1)))
|
||||
& call errquit(pname//'failed to write nh_mass',0,RTDB_ERR)
|
||||
200 if (.not.inp_a(test)) goto 10
|
||||
201 if (inp_compare(.false.,'nnhos',test)) then
|
||||
if (.not.inp_i(nnhos))
|
||||
& call errquit(pname//'failed to read number of oscillators',
|
||||
& 0,INPUT_ERR)
|
||||
if (.not.rtdb_put(rtdb,'qmd:nnhos',mt_int,1,nnhos))
|
||||
& call errquit(pname//'failed to write nnhos',0,RTDB_ERR)
|
||||
goto 200
|
||||
elseif(inp_compare(.false.,'frequency',test)) then
|
||||
if (.not.inp_f(nh_freq))
|
||||
& call errquit(pname//'failed to read frequency',0,INPUT_ERR)
|
||||
if (.not.rtdb_put(rtdb,'qmd:nhfreq',mt_dbl,1,nh_freq))
|
||||
& call errquit(pname//'failed to write nhfreq',0,RTDB_ERR)
|
||||
if (.not.inp_a(test)) goto 10
|
||||
if (inp_compare(.false.,'1/cm',test)) then
|
||||
nh_freq = nh_freq*4.5563352529120d-6
|
||||
elseif (inp_compare(.false.,'cm^-1',test)) then
|
||||
nh_freq = nh_freq*4.5563352529120d-6
|
||||
elseif (inp_compare(.false.,'ev',test)) then
|
||||
nh_freq = nh_freq*3.6749322175655d-2
|
||||
elseif (inp_compare(.false.,'au',test)) then
|
||||
nh_freq = nh_freq
|
||||
elseif (inp_compare(.false.,'ha',test)) then
|
||||
nh_freq = nh_freq
|
||||
else
|
||||
goto 201
|
||||
endif
|
||||
if (.not.rtdb_put(rtdb,'qmd:nhfreq',mt_dbl,1,nh_freq))
|
||||
& call errquit(pname//'failed to write nnhos',0,RTDB_ERR)
|
||||
goto 200
|
||||
elseif(inp_compare(.false.,'nsy',test)) then
|
||||
if (.not.inp_i(nnhsy))
|
||||
& call errquit(pname//'failed to number of Suzuki-Yoshida'//
|
||||
& ' points',0,INPUT_ERR)
|
||||
if (.not.rtdb_put(rtdb,'qmd:nnhsy',mt_int,1,nnhsy))
|
||||
& call errquit(pname//'failed to write nnhsy',0,RTDB_ERR)
|
||||
goto 200
|
||||
else
|
||||
call errquit(pname//'unknown nose-hoover option '//
|
||||
$ trim(test),0,INPUT_ERR)
|
||||
endif
|
||||
|
||||
c Berendsen or svr (Bussi, Donadio, Parrinello 2007): get tau
|
||||
else if (thermostat.eq.'berendsen'.or.thermostat.eq.'svr') then
|
||||
if (.not.inp_f(tau)) tau = 1.d3
|
||||
|
|
@ -342,10 +385,10 @@ c $Id$
|
|||
#include "mafdecls.fh"
|
||||
#include "errquit.fh"
|
||||
c
|
||||
integer nstep_nucl
|
||||
integer nstep_nucl, nnhos, nnhsy
|
||||
double precision dt_nucl
|
||||
double precision targ_temp
|
||||
double precision nh_mass(2)
|
||||
double precision nhfreq
|
||||
double precision tau
|
||||
character*32 thermostat
|
||||
character*32 integrator_nucl
|
||||
|
|
@ -378,11 +421,22 @@ c Nuclear thermostat
|
|||
& call errquit(pname//'failed to write thermostat',0,RTDB_ERR)
|
||||
endif
|
||||
c Nose-Hoover mass parameter
|
||||
if (.not.rtdb_get(rtdb,'qmd:nh_mass',mt_dbl,2,nh_mass(1)))then
|
||||
nh_mass(1)=0.1d0
|
||||
nh_mass(2)=0.1d0
|
||||
if (.not.rtdb_put(rtdb,'qmd:nh_mass',mt_dbl,2,nh_mass(1)))
|
||||
& call errquit(pname//'failed to write nh_mass',0,RTDB_ERR)
|
||||
if (.not.rtdb_get(rtdb,'qmd:nhfreq',mt_dbl,1,nhfreq))then
|
||||
nhfreq = 0.013669d0
|
||||
if (.not.rtdb_put(rtdb,'qmd:nhfreq',mt_dbl,1,nhfreq))
|
||||
& call errquit(pname//'failed to write nhfreq',0,RTDB_ERR)
|
||||
endif
|
||||
c Nose-Hoover number of oscillators
|
||||
if (.not.rtdb_get(rtdb,'qmd:nnhos',mt_int,1,nnhos))then
|
||||
nnhos = 3
|
||||
if (.not.rtdb_put(rtdb,'qmd:nnhos',mt_int,1,nnhos))
|
||||
& call errquit(pname//'failed to write nnhos',0,RTDB_ERR)
|
||||
endif
|
||||
c Nose-Hoover number of Suzuki-Yoshida points
|
||||
if (.not.rtdb_get(rtdb,'qmd:nnhsy',mt_int,1,nnhsy))then
|
||||
nnhsy = 3
|
||||
if (.not.rtdb_put(rtdb,'qmd:nnhsy',mt_int,1,nnhsy))
|
||||
& call errquit(pname//'failed to write nnhsy',0,RTDB_ERR)
|
||||
endif
|
||||
c Berendsen tau parameter, also used for svr
|
||||
if (.not.rtdb_get(rtdb,'qmd:tau',mt_dbl,1,tau)) then
|
||||
|
|
|
|||
|
|
@ -69,9 +69,9 @@ c
|
|||
c
|
||||
double precision dt2,dt4,dt8
|
||||
c
|
||||
integer i,j
|
||||
double precision scalefac
|
||||
double precision kbt
|
||||
integer i,j,ios,isy,ims,mdeg
|
||||
double precision scalefac,aa,kbt
|
||||
|
||||
c
|
||||
character*30 pname
|
||||
c
|
||||
|
|
@ -81,43 +81,75 @@ c preliminaries
|
|||
dt4 = 0.25d0*dt
|
||||
dt8 = 0.125d0*dt
|
||||
kbt = kb_au*targ_temp
|
||||
|
||||
mdeg = ndeg + 6
|
||||
if (do_linear) mdeg = mdeg - 1
|
||||
|
||||
scalefac = 1d0
|
||||
m_nh(1) = m_nh(1)*mdeg
|
||||
|
||||
g_nh(1) = (2d0*KE_in - mdeg*kbt)/m_nh(1)
|
||||
|
||||
do ims=1,3
|
||||
do isy=1,nnhsy
|
||||
v_nh(nnhos) = v_nh(nnhos) + g_nh(nnhos)*dt4*syw_nh(isy)/3
|
||||
do ios=1,nnhos-1
|
||||
aa = dexp(-dt8*syw_nh(isy)*v_nh(nnhos+1-ios)/3)
|
||||
v_nh(nnhos-ios) = v_nh(nnhos-ios)*aa**2 +
|
||||
$ dt4*syw_nh(isy)*g_nh(nnhos-ios)*aa/3
|
||||
enddo
|
||||
aa = dexp(-dt2*syw_nh(isy)*v_nh(1)/3)
|
||||
scalefac = scalefac*aa
|
||||
g_nh(1) = (2d0*KE_in*scalefac**2 - mdeg*kbt)/m_nh(1)
|
||||
do ios=1,nnhos
|
||||
r_nh(ios) = r_nh(ios) + v_nh(ios)*dt2*syw_nh(isy)/3
|
||||
enddo
|
||||
do ios=1,nnhos-1
|
||||
aa = dexp(-dt8*syw_nh(isy)*v_nh(ios+1)/3)
|
||||
v_nh(ios) = v_nh(ios)*aa**2 + dt4*syw_nh(isy)*g_nh(ios)*aa/3
|
||||
g_nh(ios+1) = (m_nh(ios)*v_nh(ios)**2 - kbt)/m_nh(ios+1)
|
||||
enddo
|
||||
v_nh(nnhos) = v_nh(nnhos) + g_nh(nnhos)*dt4*syw_nh(isy)/3
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do i=1,n
|
||||
v(1:3,i) = scalefac*v(1:3,i)
|
||||
enddo
|
||||
|
||||
m_nh(1) = m_nh(1)/mdeg
|
||||
|
||||
c
|
||||
c nose-hoover 2 chain formula: Frenkel & Smit
|
||||
g_nh(2) = (m_nh(1)*v_nh(1)*v_nh(1)-kbt)/m_nh(2)
|
||||
|
||||
c g_nh(2) = (m_nh(1)*v_nh(1)*v_nh(1)-kbt)/m_nh(2)
|
||||
c
|
||||
v_nh(2) = v_nh(2) + g_nh(2)*dt4
|
||||
v_nh(1) = v_nh(1) * dexp(-v_nh(2)*dt8)
|
||||
c v_nh(2) = v_nh(2) + g_nh(2)*dt4
|
||||
c v_nh(1) = v_nh(1) * dexp(-v_nh(2)*dt8)
|
||||
c
|
||||
g_nh(1) = (2.0d0*KE_in - ndeg*kbt)/m_nh(1)
|
||||
c g_nh(1) = (2.0d0*KE_in - ndeg*kbt)/m_nh(1)
|
||||
c
|
||||
v_nh(1) = v_nh(1) + g_nh(1)*dt4
|
||||
v_nh(1) = v_nh(1) * dexp(-v_nh(2)*dt8)
|
||||
c v_nh(1) = v_nh(1) + g_nh(1)*dt4
|
||||
c v_nh(1) = v_nh(1) * dexp(-v_nh(2)*dt8)
|
||||
c
|
||||
r_nh(1) = r_nh(1) + v_nh(1)*dt2
|
||||
r_nh(2) = r_nh(2) + v_nh(2)*dt2
|
||||
c r_nh(1) = r_nh(1) + v_nh(1)*dt2
|
||||
c r_nh(2) = r_nh(2) + v_nh(2)*dt2
|
||||
c
|
||||
scalefac = dexp(-v_nh(1)*dt2)
|
||||
c scalefac = dexp(-v_nh(1)*dt2)
|
||||
c
|
||||
c scale kinetic energy
|
||||
KE_in = KE_in*scalefac*scalefac
|
||||
c KE_in = KE_in*scalefac*scalefac
|
||||
c
|
||||
c nose-hoover 2 chain formula: Frenkel & Smit
|
||||
g_nh(1) = (2.0d0*KE_in - ndeg*kbt)/m_nh(1)
|
||||
c g_nh(1) = (2.0d0*KE_in - ndeg*kbt)/m_nh(1)
|
||||
c
|
||||
v_nh(1) = v_nh(1) * dexp(-v_nh(2)*dt8)
|
||||
v_nh(1) = v_nh(1) + g_nh(1)*dt4
|
||||
v_nh(1) = v_nh(1) * dexp(-v_nh(2)*dt8)
|
||||
c v_nh(1) = v_nh(1) * dexp(-v_nh(2)*dt8)
|
||||
c v_nh(1) = v_nh(1) + g_nh(1)*dt4
|
||||
c v_nh(1) = v_nh(1) * dexp(-v_nh(2)*dt8)
|
||||
c
|
||||
g_nh(2) = (m_nh(1)*v_nh(1)*v_nh(1)-kbt)/m_nh(2)
|
||||
c g_nh(2) = (m_nh(1)*v_nh(1)*v_nh(1)-kbt)/m_nh(2)
|
||||
c
|
||||
v_nh(2) = v_nh(2) + g_nh(2)*dt4
|
||||
c
|
||||
c scale velocity
|
||||
do i = 1,n
|
||||
v(1,i) = scalefac*v(1,i)
|
||||
v(2,i) = scalefac*v(2,i)
|
||||
v(3,i) = scalefac*v(3,i)
|
||||
end do
|
||||
c v_nh(2) = v_nh(2) + g_nh(2)*dt4
|
||||
c
|
||||
return
|
||||
end
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue