1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ########################################################
9c     ##                                                    ##
10c     ##  subroutine eimprop  --  improper dihedral energy  ##
11c     ##                                                    ##
12c     ########################################################
13c
14c
15c     "eimprop" calculates the improper dihedral potential energy
16c
17c
18      subroutine eimprop
19      use atoms
20      use bound
21      use energi
22      use group
23      use improp
24      use math
25      use torpot
26      use usage
27      implicit none
28      integer i,ia,ib,ic,id
29      real*8 e,eps,dt,fgrp
30      real*8 ideal,force
31      real*8 cosine,sine
32      real*8 rcb,angle
33      real*8 xt,yt,zt
34      real*8 xu,yu,zu
35      real*8 xtu,ytu,ztu
36      real*8 rt2,ru2,rtru
37      real*8 xia,yia,zia
38      real*8 xib,yib,zib
39      real*8 xic,yic,zic
40      real*8 xid,yid,zid
41      real*8 xba,yba,zba
42      real*8 xcb,ycb,zcb
43      real*8 xdc,ydc,zdc
44      logical proceed
45c
46c
47c     zero out improper dihedral energy
48c
49      eid = 0.0d0
50      if (niprop .eq. 0)  return
51c
52c     set tolerance for minimum distance and angle values
53c
54      eps = 0.0001d0
55c
56c     OpenMP directives for the major loop structure
57c
58!$OMP PARALLEL default(private) shared(niprop,iiprop,use,
59!$OMP& x,y,z,kprop,vprop,idihunit,eps,use_group,use_polymer)
60!$OMP& shared(eid)
61!$OMP DO reduction(+:eid) schedule(guided)
62c
63c     calculate the improper dihedral angle energy term
64c
65      do i = 1, niprop
66         ia = iiprop(1,i)
67         ib = iiprop(2,i)
68         ic = iiprop(3,i)
69         id = iiprop(4,i)
70c
71c     decide whether to compute the current interaction
72c
73         proceed = .true.
74         if (use_group)  call groups (proceed,fgrp,ia,ib,ic,id,0,0)
75         if (proceed)  proceed = (use(ia) .or. use(ib) .or.
76     &                              use(ic) .or. use(id))
77c
78c     compute the value of the improper dihedral angle
79c
80         if (proceed) then
81            xia = x(ia)
82            yia = y(ia)
83            zia = z(ia)
84            xib = x(ib)
85            yib = y(ib)
86            zib = z(ib)
87            xic = x(ic)
88            yic = y(ic)
89            zic = z(ic)
90            xid = x(id)
91            yid = y(id)
92            zid = z(id)
93            xba = xib - xia
94            yba = yib - yia
95            zba = zib - zia
96            xcb = xic - xib
97            ycb = yic - yib
98            zcb = zic - zib
99            xdc = xid - xic
100            ydc = yid - yic
101            zdc = zid - zic
102            if (use_polymer) then
103               call image (xba,yba,zba)
104               call image (xcb,ycb,zcb)
105               call image (xdc,ydc,zdc)
106            end if
107            rcb = sqrt(max(xcb*xcb+ycb*ycb+zcb*zcb,eps))
108            xt = yba*zcb - ycb*zba
109            yt = zba*xcb - zcb*xba
110            zt = xba*ycb - xcb*yba
111            xu = ycb*zdc - ydc*zcb
112            yu = zcb*xdc - zdc*xcb
113            zu = xcb*ydc - xdc*ycb
114            xtu = yt*zu - yu*zt
115            ytu = zt*xu - zu*xt
116            ztu = xt*yu - xu*yt
117            rt2 = max(xt*xt+yt*yt+zt*zt,eps)
118            ru2 = max(xu*xu+yu*yu+zu*zu,eps)
119            rtru = sqrt(rt2 * ru2)
120            cosine = (xt*xu + yt*yu + zt*zu) / rtru
121            sine = (xcb*xtu + ycb*ytu + zcb*ztu) / (rcb*rtru)
122            cosine = min(1.0d0,max(-1.0d0,cosine))
123            angle = radian * acos(cosine)
124            if (sine .lt. 0.0d0)  angle = -angle
125c
126c     set the improper dihedral parameters for this angle
127c
128            ideal = vprop(i)
129            force = kprop(i)
130            if (abs(angle+ideal) .lt. abs(angle-ideal))
131     &         ideal = -ideal
132            dt = angle - ideal
133            do while (dt .gt. 180.0d0)
134               dt = dt - 360.0d0
135            end do
136            do while (dt .lt. -180.0d0)
137               dt = dt + 360.0d0
138            end do
139c
140c     calculate the improper dihedral energy
141c
142            e = idihunit * force * dt**2
143c
144c     scale the interaction based on its group membership
145c
146            if (use_group)  e = e * fgrp
147c
148c     increment the total improper dihedral energy
149c
150            eid = eid + e
151         end if
152      end do
153c
154c     OpenMP directives for the major loop structure
155c
156!$OMP END DO
157!$OMP END PARALLEL
158      return
159      end
160