-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathNhc_Integrate_Cent.F90
136 lines (102 loc) · 3.89 KB
/
Nhc_Integrate_Cent.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
Subroutine Nhc_Integrate_Cent
use Parameters
use utility, only: norm_seq
implicit none
Double Precision :: skin, scale, dt_ys, vfact, pvfact
integer :: i, iys, iatom, Inhc
skin = 0.d0
do iatom = 1, Natom
skin = skin + fictmass(iatom,1) * norm_seq( vur(:,iatom,1) )
enddo
scale = 1.d0
! /* update the force */
fbc11(1) = (skin - gnkt)/qmcent11(1)
do Inhc = 2, nnhc
fbc11(Inhc) = (qmcent11(Inhc-1)*vbc11(Inhc-1)*vbc11(Inhc-1) &
- gkt)/qmcent11(Inhc)
enddo
! /* start multiple time step integration */
do iys = 1, nys
! /* set time increment at this loop */
dt_ys = dt*ysweight(iys)
! /* update the thermostat velocities */
vbc11(nnhc) = vbc11(nnhc) + 0.25d0*fbc11(nnhc)*dt_ys
do Inhc = 1, nnhc-1
vfact=dexp(-0.125d0* vbc11(nnhc-Inhc+1)*dt_ys)
vbc11(nnhc-Inhc) = vbc11(nnhc-Inhc)*vfact*vfact &
+ 0.25d0*fbc11(nnhc-Inhc)*vfact*dt_ys
enddo
! /* update the particle velocities */
pvfact = dexp(-0.5d0*vbc11(1)*dt_ys)
scale = scale*pvfact
! /* update the force */
fbc11(1)=(scale*scale*skin - gnkt)/qmcent11(1)
! /* update the thermostat position */
do Inhc = 1, nnhc
rbc11(Inhc) = rbc11(Inhc) &
+ 0.5d0*vbc11(Inhc)*dt_ys
enddo
! /* update the thermostat velocities */
do Inhc = 1, nnhc-1
vfact = dexp(-0.125d0*vbc11(Inhc+1)*dt_ys)
vbc11(Inhc) = vbc11(Inhc)*vfact*vfact &
+ 0.25d0*fbc11(Inhc)*vfact*dt_ys
fbc11(Inhc+1) = (qmcent11(Inhc)*vbc11(Inhc)*vbc11(Inhc) &
- gkt)/qmcent11(Inhc+1)
enddo
vbc11(nnhc) = vbc11(nnhc) + 0.25d0*fbc11(nnhc)*dt_ys
enddo
! /* update the paricle velocities */
vur(:,:,1) = vur(:,:,1) * scale
Return
End Subroutine
Subroutine nhc_integrate_cent3
Use Parameters
Implicit None
Double Precision :: skin, scale, dt_ys, vfact, pvfact
real(8) :: vrfact(3), pvrfact(3), dkinr(3)!, scaler(3)
integer :: i, iys, iatom, Inhc
! /* start multiple time step integration */
do iys = 1, Nys
! /* set time increment at this loop */
dt_ys = dt*ysweight(iys)
do iatom = 1, Natom
! /* calculate kinetic energy of centroid coordiates */
! /* for each atom */
dkinr(:) = fictmass(iatom,1)*vur(:,iatom,1)*vur(:,iatom,1)
! /* update the force */
frbc31(:,iatom,1) = (dkinr(:) - gkt)/qmcent31(1)
do Inhc = 2, Nnhc
frbc31(:,iatom,Inhc) &
= (qmcent31(Inhc-1)*vrbc31(:,iatom,Inhc-1)*vrbc31(:,iatom,Inhc-1) - gkt)/qmcent31(Inhc)
enddo
! /* update the thermostat velocities */
vrbc31(:,iatom,Nnhc) = vrbc31(:,iatom,Nnhc) + 0.25d0*frbc31(:,iatom,Nnhc)*dt_ys
do Inhc = 1, Nnhc-1
vrfact = dexp( -0.125d0*vrbc31(:,iatom,Nnhc-Inhc+1)*dt_ys )
vrbc31(:,iatom,Nnhc-Inhc) &
= vrbc31(:,iatom,Nnhc-Inhc)*vrfact*vrfact + 0.25d0*frbc31(:,iatom,Nnhc-Inhc)*vrfact*dt_ys
enddo
! /* update the particle velocities */
pvrfact(:) = dexp(-0.5d0*vrbc31(:,iatom,1)*dt_ys)
! /* update the force */
frbc31(:,iatom,1)=(pvrfact(:)*pvrfact(:)*dkinr(:) - gkt)/qmcent31(1)
! /* update the thermostat position */
do Inhc = 1, Nnhc
rbc31(:,iatom,Inhc) = rbc31(:,iatom,Inhc) + 0.5d0*vrbc31(:,iatom,Inhc)*dt_ys
enddo
! /* update the thermostat velocities */
do Inhc = 1, Nnhc-1
vrfact(:) = dexp(-0.125d0*vrbc31(:,iatom,Inhc+1)*dt_ys)
vrbc31(:,iatom,Inhc) &
= vrbc31(:,iatom,Inhc)*vrfact(:)*vrfact(:) + 0.25d0*frbc31(:,iatom,Inhc)*vrfact*dt_ys
frbc31(:,iatom,Inhc+1) &
= (qmcent31(Inhc)*vrbc31(:,iatom,Inhc)*vrbc31(:,iatom,Inhc) - gkt)/qmcent31(Inhc+1)
enddo
vrbc31(:,iatom,Nnhc) = vrbc31(:,iatom,Nnhc) + 0.25d0*frbc31(:,iatom,Nnhc)*dt_ys
! /* update the paricle velocities */
vur(:,iatom,1) = vur(:,iatom,1)*pvrfact(:)
enddo
enddo
return
End Subroutine