1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ###############################################################
9c     ##                                                           ##
10c     ##  subroutine estrbnd3  --  stretch-bend energy & analysis  ##
11c     ##                                                           ##
12c     ###############################################################
13c
14c
15c     "estrbnd3" calculates the stretch-bend potential energy;
16c     also partitions the energy among the atoms
17c
18c
19      subroutine estrbnd3
20      use action
21      use analyz
22      use angbnd
23      use angpot
24      use atomid
25      use atoms
26      use bndstr
27      use bound
28      use energi
29      use group
30      use inform
31      use iounit
32      use math
33      use strbnd
34      use usage
35      implicit none
36      integer i,j,k
37      integer istrbnd
38      integer ia,ib,ic
39      real*8 e,eps,dt
40      real*8 dr1,dr2
41      real*8 fgrp,angle
42      real*8 force1,force2
43      real*8 dot,cosine
44      real*8 xia,yia,zia
45      real*8 xib,yib,zib
46      real*8 xic,yic,zic
47      real*8 xab,yab,zab
48      real*8 xcb,ycb,zcb
49      real*8 rab,rcb
50      logical proceed
51      logical header,huge
52c
53c
54c     zero out the energy component and partitioning terms
55c
56      neba = 0
57      eba = 0.0d0
58      do i = 1, n
59         aeba(i) = 0.0d0
60      end do
61      if (nstrbnd .eq. 0)  return
62c
63c     set tolerance for minimum distance and angle values
64c
65      eps = 0.0001d0
66c
67c     print header information if debug output was requested
68c
69      header = .true.
70      if (debug .and. nstrbnd.ne.0) then
71         header = .false.
72         write (iout,10)
73   10    format (/,' Individual Stretch-Bend Interactions :',
74     &           //,' Type',18x,'Atom Names',18x,'dSB 1',
75     &              5x,'dSB 2',6x,'Energy',/)
76      end if
77c
78c     OpenMP directives for the major loop structure
79c
80!$OMP PARALLEL default(private) shared(nstrbnd,isb,iang,sbk,
81!$OMP& anat,bl,bk,use,x,y,z,stbnunit,eps,use_group,use_polymer,
82!$OMP& name,verbose,debug,header,iout)
83!$OMP& shared(eba,neba,aeba)
84!$OMP DO reduction(+:eba,neba,aeba) schedule(guided)
85c
86c     calculate the stretch-bend interaction energy term
87c
88      do istrbnd = 1, nstrbnd
89         i = isb(1,istrbnd)
90         ia = iang(1,i)
91         ib = iang(2,i)
92         ic = iang(3,i)
93         force1 = sbk(1,istrbnd)
94         force2 = sbk(2,istrbnd)
95c
96c     decide whether to compute the current interaction
97c
98         proceed = .true.
99         if (use_group)  call groups (proceed,fgrp,ia,ib,ic,0,0,0)
100         if (proceed)  proceed = (use(ia) .or. use(ib) .or. use(ic))
101c
102c     get the coordinates of the atoms in the angle
103c
104         if (proceed) then
105            xia = x(ia)
106            yia = y(ia)
107            zia = z(ia)
108            xib = x(ib)
109            yib = y(ib)
110            zib = z(ib)
111            xic = x(ic)
112            yic = y(ic)
113            zic = z(ic)
114c
115c     compute the value of the bond angle
116c
117            xab = xia - xib
118            yab = yia - yib
119            zab = zia - zib
120            xcb = xic - xib
121            ycb = yic - yib
122            zcb = zic - zib
123            if (use_polymer) then
124               call image (xab,yab,zab)
125               call image (xcb,ycb,zcb)
126            end if
127            rab = sqrt(max(xab*xab+yab*yab+zab*zab,eps))
128            rcb = sqrt(max(xcb*xcb+ycb*ycb+zcb*zcb,eps))
129            dot = xab*xcb + yab*ycb + zab*zcb
130            cosine = dot / (rab*rcb)
131            cosine = min(1.0d0,max(-1.0d0,cosine))
132            angle = radian * acos(cosine)
133            dt = angle - anat(i)
134c
135c     get the stretch-bend interaction energy
136c
137            j = isb(2,istrbnd)
138            k = isb(3,istrbnd)
139            dr1 = rab - bl(j)
140            dr2 = rcb - bl(k)
141            e = stbnunit * (force1*dr1+force2*dr2) * dt
142c
143c     scale the interaction based on its group membership
144c
145            if (use_group)  e = e * fgrp
146c
147c     increment the total stretch-bend energy
148c
149            neba = neba + 1
150            eba = eba + e
151            aeba(ib) = aeba(ib) + e
152c
153c     print a message if the energy of this interaction is large
154c
155            huge = (abs(e) .gt. 5.0d0)
156            if (debug .or. (verbose.and.huge)) then
157               if (header) then
158                  header = .false.
159                  write (iout,20)
160   20             format (/,' Individual Stretch-Bend',
161     &                       ' Interactions :',
162     &                    //,' Type',18x,'Atom Names',18x,'dSB 1',
163     &                       5x,'dSB 2',6x,'Energy',/)
164               end if
165               write (iout,30)  ia,name(ia),ib,name(ib),
166     &                          ic,name(ic),dr1*dt,dr2*dt,e
167   30          format (' StrBend',3x,3(i7,'-',a3),2x,2f10.4,f12.4)
168            end if
169         end if
170      end do
171c
172c     OpenMP directives for the major loop structure
173c
174!$OMP END DO
175!$OMP END PARALLEL
176      return
177      end
178