1c
2c
3c     #############################################################
4c     ##  COPYRIGHT (C) 2003 by Pengyu Ren & Jay William Ponder  ##
5c     ##                   All Rights Reserved                   ##
6c     #############################################################
7c
8c     #################################################################
9c     ##                                                             ##
10c     ##  subroutine etortor  --  torsion-torsion cross term energy  ##
11c     ##                                                             ##
12c     #################################################################
13c
14c
15c     "etortor" calculates the torsion-torsion potential energy
16c
17c
18      subroutine etortor
19      use atoms
20      use bitor
21      use bound
22      use energi
23      use group
24      use ktrtor
25      use math
26      use torpot
27      use tortor
28      use usage
29      implicit none
30      integer i,k,itortor
31      integer pos1,pos2
32      integer ia,ib,ic,id,ie
33      integer nlo,nhi,nt
34      integer xlo,ylo
35      real*8 e,fgrp,sign
36      real*8 angle1,angle2
37      real*8 value1,value2
38      real*8 cosine1,cosine2
39      real*8 xt,yt,zt,rt2
40      real*8 xu,yu,zu,ru2
41      real*8 xv,yv,zv,rv2
42      real*8 rtru,rurv
43      real*8 x1l,x1u
44      real*8 y1l,y1u
45      real*8 xia,yia,zia
46      real*8 xib,yib,zib
47      real*8 xic,yic,zic
48      real*8 xid,yid,zid
49      real*8 xie,yie,zie
50      real*8 xba,yba,zba
51      real*8 xdc,ydc,zdc
52      real*8 xcb,ycb,zcb
53      real*8 xed,yed,zed
54      real*8 ftt(4),ft12(4)
55      real*8 ft1(4),ft2(4)
56      logical proceed
57c
58c
59c     zero out the torsion-torsion energy
60c
61      ett = 0.0d0
62      if (ntortor .eq. 0)  return
63c
64c     OpenMP directives for the major loop structure
65c
66!$OMP PARALLEL default(private) shared(ntortor,itt,ibitor,
67!$OMP& use,x,y,z,tnx,ttx,tny,tty,tbf,tbx,tby,tbxy,ttorunit,
68!$OMP& use_group,use_polymer)
69!$OMP& shared(ett)
70!$OMP DO reduction(+:ett) schedule(guided)
71c
72c     calculate the torsion-torsion interaction energy term
73c
74      do itortor = 1, ntortor
75         i = itt(1,itortor)
76         k = itt(2,itortor)
77         if (itt(3,itortor) .eq. 1) then
78            ia = ibitor(1,i)
79            ib = ibitor(2,i)
80            ic = ibitor(3,i)
81            id = ibitor(4,i)
82            ie = ibitor(5,i)
83         else
84            ia = ibitor(5,i)
85            ib = ibitor(4,i)
86            ic = ibitor(3,i)
87            id = ibitor(2,i)
88            ie = ibitor(1,i)
89         end if
90c
91c     decide whether to compute the current interaction
92c
93         proceed = .true.
94         if (use_group)  call groups (proceed,fgrp,ia,ib,ic,id,ie,0)
95         if (proceed)  proceed = (use(ia) .or. use(ib) .or. use(ic)
96     &                               .or. use(id) .or. use(ie))
97c
98c     compute the values of the torsional angles
99c
100         if (proceed) then
101            xia = x(ia)
102            yia = y(ia)
103            zia = z(ia)
104            xib = x(ib)
105            yib = y(ib)
106            zib = z(ib)
107            xic = x(ic)
108            yic = y(ic)
109            zic = z(ic)
110            xid = x(id)
111            yid = y(id)
112            zid = z(id)
113            xie = x(ie)
114            yie = y(ie)
115            zie = z(ie)
116            xba = xib - xia
117            yba = yib - yia
118            zba = zib - zia
119            xcb = xic - xib
120            ycb = yic - yib
121            zcb = zic - zib
122            xdc = xid - xic
123            ydc = yid - yic
124            zdc = zid - zic
125            xed = xie - xid
126            yed = yie - yid
127            zed = zie - zid
128            if (use_polymer) then
129               call image (xba,yba,zba)
130               call image (xcb,ycb,zcb)
131               call image (xdc,ydc,zdc)
132               call image (xed,yed,zed)
133            end if
134            xt = yba*zcb - ycb*zba
135            yt = zba*xcb - zcb*xba
136            zt = xba*ycb - xcb*yba
137            xu = ycb*zdc - ydc*zcb
138            yu = zcb*xdc - zdc*xcb
139            zu = xcb*ydc - xdc*ycb
140            rt2 = xt*xt + yt*yt + zt*zt
141            ru2 = xu*xu + yu*yu + zu*zu
142            rtru = sqrt(rt2 * ru2)
143            xv = ydc*zed - yed*zdc
144            yv = zdc*xed - zed*xdc
145            zv = xdc*yed - xed*ydc
146            rv2 = xv*xv + yv*yv + zv*zv
147            rurv = sqrt(ru2 * rv2)
148            if (rtru.ne.0.0d0 .and. rurv.ne.0.0d0) then
149               cosine1 = (xt*xu + yt*yu + zt*zu) / rtru
150               cosine1 = min(1.0d0,max(-1.0d0,cosine1))
151               angle1 = radian * acos(cosine1)
152               sign = xba*xu + yba*yu + zba*zu
153               if (sign .lt. 0.0d0)  angle1 = -angle1
154               value1 = angle1
155               cosine2 = (xu*xv + yu*yv + zu*zv) / rurv
156               cosine2 = min(1.0d0,max(-1.0d0,cosine2))
157               angle2 = radian * acos(cosine2)
158               sign = xcb*xv + ycb*yv + zcb*zv
159               if (sign .lt. 0.0d0)  angle2 = -angle2
160               value2 = angle2
161c
162c     check for inverted chirality at the central atom
163c
164               call chkttor (ib,ic,id,sign,value1,value2)
165c
166c     use bicubic interpolation to compute spline values
167c
168               nlo = 1
169               nhi = tnx(k)
170               do while (nhi-nlo .gt. 1)
171                  nt = (nhi+nlo) / 2
172                  if (ttx(nt,k) .gt. value1) then
173                     nhi = nt
174                  else
175                     nlo = nt
176                  end if
177               end do
178               xlo = nlo
179               nlo = 1
180               nhi = tny(k)
181               do while (nhi-nlo .gt. 1)
182                  nt = (nhi + nlo)/2
183                  if (tty(nt,k) .gt. value2) then
184                     nhi = nt
185                  else
186                     nlo = nt
187                  end if
188               end do
189               ylo = nlo
190               x1l = ttx(xlo,k)
191               x1u = ttx(xlo+1,k)
192               y1l = tty(ylo,k)
193               y1u = tty(ylo+1,k)
194               pos2 = ylo*tnx(k) + xlo
195               pos1 = pos2 - tnx(k)
196               ftt(1) = tbf(pos1,k)
197               ftt(2) = tbf(pos1+1,k)
198               ftt(3) = tbf(pos2+1,k)
199               ftt(4) = tbf(pos2,k)
200               ft1(1) = tbx(pos1,k)
201               ft1(2) = tbx(pos1+1,k)
202               ft1(3) = tbx(pos2+1,k)
203               ft1(4) = tbx(pos2,k)
204               ft2(1) = tby(pos1,k)
205               ft2(2) = tby(pos1+1,k)
206               ft2(3) = tby(pos2+1,k)
207               ft2(4) = tby(pos2,k)
208               ft12(1) = tbxy(pos1,k)
209               ft12(2) = tbxy(pos1+1,k)
210               ft12(3) = tbxy(pos2+1,k)
211               ft12(4) = tbxy(pos2,k)
212               call bcuint (ftt,ft1,ft2,ft12,x1l,x1u,
213     &                      y1l,y1u,value1,value2,e)
214               e = ttorunit * e
215c
216c     scale the interaction based on its group membership
217c
218               if (use_group)  e = e * fgrp
219c
220c     increment the total torsion-torsion energy
221c
222               ett = ett + e
223            end if
224         end if
225      end do
226c
227c     OpenMP directives for the major loop structure
228c
229!$OMP END DO
230!$OMP END PARALLEL
231      return
232      end
233c
234c
235c     ###############################################################
236c     ##                                                           ##
237c     ##  subroutine chkttor  --  check torsion-torsion chirality  ##
238c     ##                                                           ##
239c     ###############################################################
240c
241c
242c     "chkttor" tests the attached atoms at a torsion-torsion central
243c     site and changes the sign of the torsion values if the site has
244c     opposite chirality to that for the original parameter
245c
246c     note that the sign convention used in this version is correct
247c     for phi-psi torsion-torsion interactions defined for L-amino
248c     acids in a protein force field; the code may need to be altered
249c     for other chiral torsion-torsion situations, such as the sugar
250c     rings in nucleic acids
251c
252c
253      subroutine chkttor (ib,ic,id,sign,value1,value2)
254      use atomid
255      use atoms
256      use couple
257      implicit none
258      integer i,j,k,m
259      integer ia,ib,ic,id
260      real*8 sign
261      real*8 value1
262      real*8 value2
263      real*8 xac,yac,zac
264      real*8 xbc,ybc,zbc
265      real*8 xdc,ydc,zdc
266      real*8 c1,c2,c3,vol
267c
268c
269c     test for chirality at the central torsion-torsion site
270c
271      sign = 1.0d0
272      if (n12(ic) .eq. 4) then
273         j = 0
274         do i = 1, 4
275            m = i12(i,ic)
276            if (m.ne.ib .and. m.ne.id) then
277               if (j .eq. 0) then
278                  j = m
279               else
280                  k = m
281               end if
282            end if
283         end do
284         ia = 0
285         if (type(j) .gt. type(k))  ia = j
286         if (type(k) .gt. type(j))  ia = k
287         if (atomic(j) .gt. atomic(k))  ia = j
288         if (atomic(k) .gt. atomic(j))  ia = k
289c
290c     compute the signed parallelpiped volume at central site
291c
292         if (ia .ne. 0) then
293            xac = x(ia) - x(ic)
294            yac = y(ia) - y(ic)
295            zac = z(ia) - z(ic)
296            xbc = x(ib) - x(ic)
297            ybc = y(ib) - y(ic)
298            zbc = z(ib) - z(ic)
299            xdc = x(id) - x(ic)
300            ydc = y(id) - y(ic)
301            zdc = z(id) - z(ic)
302            c1 = ybc*zdc - zbc*ydc
303            c2 = ydc*zac - zdc*yac
304            c3 = yac*zbc - zac*ybc
305            vol = xac*c1 + xbc*c2 + xdc*c3
306c
307c     invert the angle values if chirality has an inverted sign
308c
309            if (vol .lt. 0.0d0) then
310               sign = -1.0d0
311               value1 = -value1
312               value2 = -value2
313            end if
314         end if
315      end if
316      return
317      end
318