forked from TinkerTools/tinker
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy patheangang.f
155 lines (155 loc) · 4.2 KB
/
eangang.f
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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
c
c
c ###################################################
c ## COPYRIGHT (C) 1993 by Jay William Ponder ##
c ###################################################
c
c #############################################################
c ## ##
c ## subroutine eangang -- angle-angle cross term energy ##
c ## ##
c #############################################################
c
c
c "eangang" calculates the angle-angle potential energy
c
c
subroutine eangang
use angang
use angbnd
use angpot
use atoms
use bound
use energi
use group
use math
use usage
implicit none
integer i,k,iangang
integer ia,ib,ic,id,ie
real*8 e,angle
real*8 eps,fgrp
real*8 dt1,dt2
real*8 dot,cosine
real*8 xia,yia,zia
real*8 xib,yib,zib
real*8 xic,yic,zic
real*8 xid,yid,zid
real*8 xie,yie,zie
real*8 xab,yab,zab
real*8 xcb,ycb,zcb
real*8 xdb,ydb,zdb
real*8 xeb,yeb,zeb
real*8 rab2,rcb2
real*8 rdb2,reb2
logical proceed
c
c
c zero out the angle-angle cross term energy
c
eaa = 0.0d0
if (nangang .eq. 0) return
c
c set tolerance for minimum distance and angle values
c
eps = 0.0001d0
c
c OpenMP directives for the major loop structure
c
!$OMP PARALLEL default(private) shared(nangang,iaa,iang,
!$OMP& use,x,y,z,anat,kaa,aaunit,eps,use_group,use_polymer)
!$OMP& shared(eaa)
!$OMP DO reduction(+:eaa) schedule(guided)
c
c calculate the angle-angle interaction energy term
c
do iangang = 1, nangang
i = iaa(1,iangang)
k = iaa(2,iangang)
ia = iang(1,i)
ib = iang(2,i)
ic = iang(3,i)
id = iang(1,k)
ie = iang(3,k)
c
c decide whether to compute the current interaction
c
proceed = .true.
if (use_group) call groups (proceed,fgrp,ia,ib,ic,id,ie,0)
if (proceed) proceed = (use(ia) .or. use(ib) .or. use(ic)
& .or. use(id) .or. use(ie))
c
c get the coordinates of the atoms in the angle
c
if (proceed) then
xia = x(ia)
yia = y(ia)
zia = z(ia)
xib = x(ib)
yib = y(ib)
zib = z(ib)
xic = x(ic)
yic = y(ic)
zic = z(ic)
xid = x(id)
yid = y(id)
zid = z(id)
xie = x(ie)
yie = y(ie)
zie = z(ie)
c
c compute the values of the two bond angles
c
xab = xia - xib
yab = yia - yib
zab = zia - zib
xcb = xic - xib
ycb = yic - yib
zcb = zic - zib
xdb = xid - xib
ydb = yid - yib
zdb = zid - zib
xeb = xie - xib
yeb = yie - yib
zeb = zie - zib
if (use_polymer) then
call image (xab,yab,zab)
call image (xcb,ycb,zcb)
call image (xdb,ydb,zdb)
call image (xeb,yeb,zeb)
end if
rab2 = max(xab*xab+yab*yab+zab*zab,eps)
rcb2 = max(xcb*xcb+ycb*ycb+zcb*zcb,eps)
rdb2 = max(xdb*xdb+ydb*ydb+zdb*zdb,eps)
reb2 = max(xeb*xeb+yeb*yeb+zeb*zeb,eps)
dot = xab*xcb + yab*ycb + zab*zcb
cosine = dot / sqrt(rab2*rcb2)
cosine = min(1.0d0,max(-1.0d0,cosine))
angle = radian * acos(cosine)
dt1 = angle - anat(i)
dot = xdb*xeb + ydb*yeb + zdb*zeb
cosine = dot / sqrt(rdb2*reb2)
cosine = min(1.0d0,max(-1.0d0,cosine))
angle = radian * acos(cosine)
dt2 = angle - anat(k)
c
c get the angle-angle interaction energy
c
e = aaunit * kaa(iangang) * dt1 * dt2
c
c scale the interaction based on its group membership
c
if (use_group) e = e * fgrp
c
c increment the total angle-angle energy
c
eaa = eaa + e
end if
end do
c
c OpenMP directives for the major loop structure
c
!$OMP END DO
!$OMP END PARALLEL
return
end