1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ##############################################################
9c     ##                                                          ##
10c     ##  subroutine estrbnd  --  stretch-bend cross term energy  ##
11c     ##                                                          ##
12c     ##############################################################
13c
14c
15c     "estrbnd" calculates the stretch-bend potential energy
16c
17c
18      subroutine estrbnd
19      use angbnd
20      use angpot
21      use atoms
22      use bndstr
23      use bound
24      use energi
25      use group
26      use math
27      use strbnd
28      use usage
29      implicit none
30      integer i,j,k
31      integer istrbnd
32      integer ia,ib,ic
33      real*8 e,eps,dt
34      real*8 dr1,dr2
35      real*8 fgrp,angle
36      real*8 force1,force2
37      real*8 dot,cosine
38      real*8 xia,yia,zia
39      real*8 xib,yib,zib
40      real*8 xic,yic,zic
41      real*8 xab,yab,zab
42      real*8 xcb,ycb,zcb
43      real*8 rab,rcb
44      logical proceed
45c
46c
47c     zero out the stretch-bend cross term energy
48c
49      eba = 0.0d0
50      if (nstrbnd .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(nstrbnd,isb,iang,sbk,
59!$OMP& anat,bl,bk,use,x,y,z,stbnunit,eps,use_group,use_polymer)
60!$OMP& shared(eba)
61!$OMP DO reduction(+:eba) schedule(guided)
62c
63c     calculate the stretch-bend interaction energy term
64c
65      do istrbnd = 1, nstrbnd
66         i = isb(1,istrbnd)
67         ia = iang(1,i)
68         ib = iang(2,i)
69         ic = iang(3,i)
70         force1 = sbk(1,istrbnd)
71         force2 = sbk(2,istrbnd)
72c
73c     decide whether to compute the current interaction
74c
75         proceed = .true.
76         if (use_group)  call groups (proceed,fgrp,ia,ib,ic,0,0,0)
77         if (proceed)  proceed = (use(ia) .or. use(ib) .or. use(ic))
78c
79c     get the coordinates of the atoms in the angle
80c
81         if (proceed) then
82            xia = x(ia)
83            yia = y(ia)
84            zia = z(ia)
85            xib = x(ib)
86            yib = y(ib)
87            zib = z(ib)
88            xic = x(ic)
89            yic = y(ic)
90            zic = z(ic)
91c
92c     compute the value of the bond angle
93c
94            xab = xia - xib
95            yab = yia - yib
96            zab = zia - zib
97            xcb = xic - xib
98            ycb = yic - yib
99            zcb = zic - zib
100            if (use_polymer) then
101               call image (xab,yab,zab)
102               call image (xcb,ycb,zcb)
103            end if
104            rab = sqrt(max(xab*xab+yab*yab+zab*zab,eps))
105            rcb = sqrt(max(xcb*xcb+ycb*ycb+zcb*zcb,eps))
106            dot = xab*xcb + yab*ycb + zab*zcb
107            cosine = dot / (rab*rcb)
108            cosine = min(1.0d0,max(-1.0d0,cosine))
109            angle = radian * acos(cosine)
110            dt = angle - anat(i)
111c
112c     get the stretch-bend interaction energy
113c
114            j = isb(2,istrbnd)
115            k = isb(3,istrbnd)
116            dr1 = rab - bl(j)
117            dr2 = rcb - bl(k)
118            e = stbnunit * (force1*dr1+force2*dr2) * dt
119c
120c     scale the interaction based on its group membership
121c
122            if (use_group)  e = e * fgrp
123c
124c     increment the total stretch-bend energy
125c
126            eba = eba + e
127         end if
128      end do
129c
130c     OpenMP directives for the major loop structure
131c
132!$OMP END DO
133!$OMP END PARALLEL
134      return
135      end
136