chemshell_qchem / blob
1 SUBROUTINE nvt_nhc(lshmov, idnode, imcon, mxnode, NHLen, dof, ntot, natms, &
2 & nscons, ntcons, etadot, eta, consv, engke, taut, sigma, tolnce, tstep, &
3 & vircon, lashap, lishap, listcon, listme, listot, lstfrz, buffer, cell, &
4 & dxt, dyt, dzt, fxx, fyy, fzz, prmcon, txx, tyy, tzz, vxx, vyy, vzz, &
5 & weight, xxt, xxx, yyt, yyy, zzt, zzz, stress, &
6 & listrcrd, rcrd, cfor, ntrcrd)
7 !===============================================================================
8 ! Subroutine for integrating equations of motion in molecular dynamics to be used
9 ! with ChemShell/DL_POLY.
10 !
11 ! Implements the Nose-Hoover chain thermostat
12 ! (Martyna, Klein, Tuckerman, J. Chem. Phys 1992, 97, 2635-2643).
13 !
14 ! Leap-frog-type integration scheme according to the "LF-2" integrator from
15 ! Jang, Voth, J. Chem. Phys. 1997, 107, 9514-9526.
16 !
17 ! Code tries to follow the nvt_h1 routine from dl_poly.
18 !
19 ! Hans Martin Senn (hms), MPI fuer Kohlenforschung
20 ! - Initial version (June 2003)
21 ! - Remove explicit declaration of length of INTEGER and LOGICAL to allow for compiler switches
22 ! for 64-bit Linux system (Jan 2004)
23 !
24 !===============================================================================
25 IMPLICIT NONE
26 INCLUDE 'dl_params.f90.inc'
28 INTEGER, INTENT(in) :: ntot, ntcons, ntrcrd
29 INTEGER, INTENT(in) :: NHLen ! Length of NH chain
31 REAL(8), INTENT(in) :: taut, sigma, dof
33 ! These are INOUT only because they are passed on to rdshake_1
34 ! where we don't (want to) know what is going to happen to them
35 INTEGER, INTENT(inout) :: idnode, imcon, mxnode, natms, nscons
36 INTEGER, INTENT(inout) :: lashap(mxproc), lishap(mxlshp), listcon(nscons,3)
37 INTEGER, INTENT(inout) :: listme(mxatms), listot(mxatms), lstfrz(natms), listrcrd(6)
39 REAL(8), INTENT(inout) :: rcrd, cfor, tolnce, tstep
40 REAL(8), INTENT(inout) :: cell(9), buffer(10000)
41 REAL(8), INTENT(inout) :: dxt(mxcons), dyt(mxcons), dzt(mxcons)
42 REAL(8), INTENT(inout) :: prmcon(mxtcon), weight(natms)
43 REAL(8), INTENT(inout) :: txx(mxatms), tyy(mxatms), tzz(mxatms)
44 REAL(8), INTENT(inout) :: xxt(mxatms)
45 REAL(8), INTENT(inout) :: yyt(mxatms)
46 REAL(8), INTENT(inout) :: zzt(mxatms)
48 LOGICAL, INTENT(inout) :: lshmov
50 ! "True" INOUTs
51 REAL(8), INTENT(inout) :: stress(9), etadot(NHLen), eta(NHLen)
52 REAL(8), INTENT(inout) :: xxx(ntot), yyy(ntot), zzz(ntot)
53 REAL(8), INTENT(inout) :: fxx(natms), fyy(natms), fzz(natms)
54 REAL(8), INTENT(inout) :: vxx(natms), vyy(natms), vzz(natms)
56 REAL(8), INTENT(out) :: vircon, consv, engke
58 ! Local variables
59 INTEGER :: nbuff, iatm0, iatm1, Nat1
61 REAL(8) :: stres1(9) = 0.d0, viracc, Ekin, Q(NHLen), G(NHLen)
62 REAL(8) :: d(nscons,3), R(ntot,3), etaIn(NHLen), etadot1(NHLen+1)
63 REAL(8) :: etadotIn(NHLen), eta0(NHLen), etadot0(NHLen)
64 REAL(8), ALLOCATABLE :: F0(:,:), vmh(:,:), vph(:,:), M(:,:), Rp(:,:), v0(:,:)
66 LOGICAL :: safe
67 !-------------------------------------------------------------------------------
68 ! Preliminaries
69 !-------------------------------------------------------------------------------
71 #ifdef VAMPIR
72 CALL VTBEGIN(53, ierr)
73 #endif
75 safe = .FALSE.
76 nbuff = mxbuff
78 ! Constraint virial
79 vircon = 0.d0
81 ! Set up block indices: First and last atom on local node
82 iatm0 = (idnode*natms)/mxnode + 1
83 iatm1 = ((idnode + 1)*natms)/mxnode
85 ! Local number of atoms
86 Nat1 = iatm1 - iatm0 + 1
88 ! Allocate arrays for local data
89 ALLOCATE(vmh(Nat1,3))
90 ALLOCATE(vph(Nat1,3))
91 ALLOCATE(F0(Nat1,3))
92 ALLOCATE(M(Nat1,3))
93 ALLOCATE(v0(Nat1,3))
95 ! Store full initial positions R(0)
96 R(:,1) = xxx
97 R(:,2) = yyy
98 R(:,3) = zzz
100 ! Store local initial velocities v(-Dt/2)
101 vmh(:,1) = vxx(iatm0:iatm1)
102 vmh(:,2) = vyy(iatm0:iatm1)
103 vmh(:,3) = vzz(iatm0:iatm1)
105 ! Store local initial forces F(0)
106 F0(:,1) = fxx(iatm0:iatm1)
107 F0(:,2) = fyy(iatm0:iatm1)
108 F0(:,3) = fzz(iatm0:iatm1)
110 ! Store initial thermostat variables
111 etaIn = eta
112 etadotIn = etadot
114 ! Build local mass matrix
115 M = SPREAD(weight(iatm0:iatm1), DIM=2, NCOPIES=3)
116 !-------------------------------------------------------------------------------
117 ! Bond vectors for SHAKE
118 !-------------------------------------------------------------------------------
119 IF ((ntcons > 0).OR.(ntrcrd > 0)) THEN
120 ! Construct current bond vectors for all atoms
121 d(1:nscons,:) = R(listcon(1:nscons,2),:) - R(listcon(1:nscons,3),:)
123 ! Periodic boundary condition for bond vectors
124 CALL images(imcon, 0, 1, nscons, cell, d(:,1), d(:,2), d(:,3))
125 end IF
127 !-------------------------------------------------------------------------------
128 ! Propagate system variables
129 !-------------------------------------------------------------------------------
130 ! Propagate local system velocities v(-Dt/2) -> v(+Dt/2)
131 vph = vmh*EXP(-tstep*etadot(1)) + tstep*EXP(-tstep*etadot(1)/2.d0)*F0/M
133 ! Propagate local system positions R(0) -> R(+Dt)
134 R(iatm0:iatm1,:) = R(iatm0:iatm1,:) + tstep*vph
136 !-------------------------------------------------------------------------------
137 ! Handle constraints
138 !-------------------------------------------------------------------------------
139 IF ((ntcons > 0).OR.(ntrcrd > 0)) THEN
140 ! Apply constraints procedures
142 ! Store propagated local positions
143 ALLOCATE(Rp(Nat1,3))
144 Rp = R(iatm0:iatm1,:)
146 ! Exchange of configuration data
147 IF (mxnode > 1) THEN
148 CALL merge(idnode, mxnode, natms, nbuff, R(:,1), R(:,2), R(:,3), buffer)
149 ENDIF
151 ! Apply constraints corrections (global)
152 CALL rdshake_1(safe, lshmov, idnode, imcon, mxnode, natms, nscons,&
153 & tolnce, tstep, viracc, lashap, lishap, listcon, listme,&
154 & listot, lstfrz, buffer, cell, dxt, d(:,1), dyt, d(:,2), dzt, d(:,3),&
155 & prmcon, txx, tyy, tzz, weight, xxt, R(:,1), yyt, R(:,2), zzt, R(:,3),&
156 & stres1, listrcrd, rcrd, cfor)
158 ! Contribution to virial
159 vircon = vircon + viracc
161 ! Force correction (local)
162 ! N.B. Rp: propagated positions
163 ! R: propagated AND SHAKE-corrected positions
164 F0 = F0 + M/tstep**2 * (R(iatm0:iatm1,:) - Rp(:,:))
166 DEALLOCATE(Rp)
168 ! Recalculate propagated velocities
169 vph = vmh*EXP(-tstep*etadot(1)) + tstep*EXP(-tstep*etadot(1)/2.d0)*F0/M
170 ENDIF
172 !-------------------------------------------------------------------------------
173 ! Propagate thermostat variables
174 !-------------------------------------------------------------------------------
175 ! Propagate ODD thermostat positions eta(-Dt/2) -> eta(+Dt/2)
176 eta(1:NHLen:2) = eta(1:NHLen:2) + tstep*etadot(1:NHLen:2)
178 ! Thermostat masses (target kinetic energy sigma = dof*k*T/2)
179 Q(1) = 2.d0*sigma*taut**2
180 Q(2:NHLen) = Q(1)/dof
182 ! EVEN thermostat forces G(0)
183 G(2:NHLen:2) = Q(1:NHLen:2)*etadot(1:NHLen:2)**2 - 2.d0*sigma/dof
185 ! The last thermostat in the chain is "thermostatted" with zero friction
186 etadot1 = (/etadot, 0.d0/)
188 ! Propagate EVEN thermostat velocities etadot(-Dt/2) -> etadot(+Dt/2)
189 etadot(2:NHLen:2) = etadot(2:NHLen:2)*EXP(-tstep*etadot1(3:NHLen+1:2)) &
190 & + tstep*EXP(-tstep*etadot1(3:NHLen+1:2)/2.d0)*G(2:NHLen:2)/Q(2:NHLen:2)
192 ! Propagate EVEN thermostat positions eta(0) -> eta(+Dt)
193 eta(2:NHLen:2) = eta(2:NHLen:2) + tstep*etadot(2:NHLen:2)
195 ! System kinetic energy at +Dt/2
196 Ekin = 0.5d0*SUM(M*vph**2)
197 IF (mxnode > 1) CALL gdsum(Ekin, 1, buffer)
199 ! ODD thermostat forces G(+Dt/2)
200 G(1) = 2.d0*(Ekin - sigma)
201 G(3:NHLen:2) = Q(2:NHLen:2)*etadot(2:NHLen:2)**2 - 2.d0*sigma/dof
203 ! The last thermostat in the chain is "thermostatted" with zero friction
204 etadot1 = (/etadot, 0.d0/)
206 ! Propagate ODD thermostat velocities etadot(0) -> etadot(+Dt)
207 etadot(1:NHLen:2) = etadot(1:NHLen:2)*EXP(-tstep*etadot1(2:NHLen+1:2)) &
208 & + tstep*EXP(-tstep*etadot1(2:NHLen+1:2)/2.d0)*G(1:NHLen:2)/Q(1:NHLen:2)
209 !-------------------------------------------------------------------------------
210 ! Finish up
211 !-------------------------------------------------------------------------------
212 ! Calculate all quantities (including those defined at half-steps) at 0
213 v0 = 0.5d0*(vmh + vph)
215 eta0(1:NHLen:2) = 0.5d0*(etaIn(1:NHLen:2) + eta(1:NHLen:2))
216 eta0(2:NHLen:2) = etaIn(2:NHLen:2)
218 etadot0(1:NHLen:2) = etadotIn(1:NHLen:2)
219 etadot0(2:NHLen:2) = 0.5d0*(etadotIn(2:NHLen:2) + etadot(2:NHLen:2))
222 #ifdef STRESS
223 ! Kinetic contribution to stress tensor (local)
224 stres1(1) = stres1(1) + SUM(M(:,1)*v0(:,1)**2)
225 stres1(2) = stres1(2) + SUM(M(:,1)*v0(:,1)*v0(:,2))
226 stres1(3) = stres1(3) + SUM(M(:,1)*v0(:,1)*v0(:,3))
227 stres1(5) = stres1(5) + SUM(M(:,1)*v0(:,2)**2)
228 stres1(6) = stres1(6) + SUM(M(:,1)*v0(:,2)*v0(:,3))
229 stres1(9) = stres1(9) + SUM(M(:,1)*v0(:,3)**2)
230 #endif
232 ! System kinetic energy at 0
233 engke = 0.5d0*SUM(M*v0**2)
234 IF (mxnode > 1) CALL gdsum(engke, 1, buffer)
236 ! Thermostat contributions to conserved energy
237 consv = 0.5d0*SUM(Q*etadot0**2) + 2.d0*sigma*eta0(1) + 2.d0*sigma/dof*SUM(eta0(2:NHLen))
239 ! Periodic boundary conditions
240 CALL images(imcon, idnode, mxnode, natms, cell, R(:,1), R(:,2), R(:,3))
242 ! Put data back into dl_poly structures
243 xxx = R(:,1)
244 yyy = R(:,2)
245 zzz = R(:,3)
247 vxx(iatm0:iatm1) = vph(:,1)
248 vyy(iatm0:iatm1) = vph(:,2)
249 vzz(iatm0:iatm1) = vph(:,3)
251 fxx(iatm0:iatm1) = F0(:,1)
252 fyy(iatm0:iatm1) = F0(:,2)
253 fzz(iatm0:iatm1) = F0(:,3)
255 ! Global exchange of configuration data
256 IF(mxnode > 1)THEN
257 nbuff = mxbuff
258 CALL merge(idnode, mxnode, natms, nbuff, xxx, yyy, zzz, buffer)
259 CALL merge(idnode, mxnode, natms, nbuff, vxx, vyy, vzz, buffer)
261 IF((ntcons > 0).OR.(ntrcrd > 0)) THEN
262 CALL merge(idnode, mxnode, natms, nbuff, fxx, fyy, fzz, buffer)
263 ENDIF
264 ENDIF
266 #ifdef STRESS
267 ! Complete stress tensor
268 IF(mxnode > 1) CALL gdsum(stres1, 9, buffer)
270 stress(1) = stress(1) + stres1(1)
271 stress(2) = stress(2) + stres1(2)
272 stress(3) = stress(3) + stres1(3)
273 stress(4) = stress(4) + stres1(2)
274 stress(5) = stress(5) + stres1(5)
275 stress(6) = stress(6) + stres1(6)
276 stress(7) = stress(7) + stres1(3)
277 stress(8) = stress(8) + stres1(6)
278 stress(9) = stress(9) + stres1(9)
279 #endif
281 #ifdef VAMPIR
282 CALL VTEND(53, ierr)
283 #endif
285 DEALLOCATE(vmh)
286 DEALLOCATE(vph)
287 DEALLOCATE(F0)
288 DEALLOCATE(M)
289 DEALLOCATE(v0)
291 RETURN
292 END SUBROUTINE nvt_nhc
