-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmapverlet.f90
102 lines (78 loc) · 2 KB
/
mapverlet.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
subroutine mapverlet(hel,x0,p0)
use global
implicit none
real*8 xrhs(nstate),prhs(nstate)
real*8 force(nstate)
real*8 xsumdum,psumdum,ddt
real*8 hel(nstate,nstate)
real*8 x0(nstate),p0(nstate)
integer i,n,m,nlit
nlit=20
ddt=0.5*dt/nlit
do i=1,nlit
do n=1,nstate
xrhs(n)=hel(n,n)*p0(n)
prhs(n)=-hel(n,n)*x0(n)
end do
do n=1,nstate
xsumdum=0.
psumdum=0.
do m=1,nstate
if(m.ne.n) then
xsumdum=xsumdum+hel(n,m)*p0(m)
psumdum=psumdum+hel(n,m)*x0(m)
end if
end do
xrhs(n)=xrhs(n)+xsumdum
prhs(n)=prhs(n)-psumdum
end do
! Generate current second derivatives
do n=1,nstate
force(n)=hel(n,n)*prhs(n)
end do
do n=1,nstate
xsumdum=0.
do m=1,nstate
if(m.ne.n) then
xsumdum=xsumdum+hel(n,m)*prhs(m)
end if
end do
force(n)=force(n)+xsumdum
end do
! other mapping force from gaussian overlap
! do n=1,nstate
! xrhs(n,i)=xrhs(n,i)+p0(n,i)*real(nb/beta)
! prhs(n,i)=prhs(n,i)-x0(n,i)*real(nb/beta)
! force(n,i)=force(n,i)-x0(n,i)*(nb/beta)**2
! end do
! Advance Ver step
do n=1,nstate
x0(n)=x0(n)+xrhs(n)*ddt+0.5*force(n)*ddt*ddt
p0(n)=p0(n)+0.5*prhs(n)*ddt
end do
! Compute new first derivatives
do n=1,nstate
prhs(n)=-hel(n,n)*x0(n)
end do
do n=1,nstate
psumdum=0.
do m=1,nstate
if(m.ne.n) then
psumdum=psumdum+hel(n,m)*x0(m)
end if
end do
prhs(n)=prhs(n)-psumdum
end do
! other mapping force from gaussian
! do n=1,nstate
! prhs(n,i)=prhs(n,i)-x0(n,i)*real(nb/beta)
! force(n,i)=force(n,i)-x0(n,i)*(nb/beta)**2
! end do
! Advance let step
do n=1,nstate
p0(n)=p0(n)+0.5*prhs(n)*ddt
end do
! end loop of the mapping propagation
end do
return
end subroutine mapverlet