1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ################################################################
9c     ##                                                            ##
10c     ##  subroutine eurey1  --  bond stretch energy & derivatives  ##
11c     ##                                                            ##
12c     ################################################################
13c
14c
15c     "eurey1" calculates the Urey-Bradley interaction energy and
16c     its first derivatives with respect to Cartesian coordinates
17c
18c
19      subroutine eurey1
20      use atoms
21      use bound
22      use deriv
23      use energi
24      use group
25      use urey
26      use urypot
27      use usage
28      use virial
29      implicit none
30      integer i,ia,ic
31      real*8 e,de
32      real*8 ideal,force
33      real*8 dt,dt2,deddt,fgrp
34      real*8 dedx,dedy,dedz
35      real*8 xac,yac,zac,rac
36      real*8 vxx,vyy,vzz
37      real*8 vyx,vzx,vzy
38      logical proceed
39c
40c
41c     zero out the Urey-Bradley energy and first derivatives
42c
43      eub = 0.0d0
44      do i = 1, n
45         deub(1,i) = 0.0d0
46         deub(2,i) = 0.0d0
47         deub(3,i) = 0.0d0
48      end do
49      if (nurey .eq. 0)  return
50c
51c     OpenMP directives for the major loop structure
52c
53!$OMP PARALLEL default(private) shared(nurey,iury,ul,uk,
54!$OMP& use,x,y,z,cury,qury,ureyunit,use_group,use_polymer)
55!$OMP& shared(eub,deub,vir)
56!$OMP DO reduction(+:eub,deub,vir) schedule(guided)
57c
58c     calculate the Urey-Bradley 1-3 energy and first derivatives
59c
60      do i = 1, nurey
61         ia = iury(1,i)
62         ic = iury(3,i)
63         ideal = ul(i)
64         force = uk(i)
65c
66c     decide whether to compute the current interaction
67c
68         proceed = .true.
69         if (use_group)  call groups (proceed,fgrp,ia,ic,0,0,0,0)
70         if (proceed)  proceed = (use(ia) .or. use(ic))
71c
72c     compute the value of the 1-3 distance deviation
73c
74         if (proceed) then
75            xac = x(ia) - x(ic)
76            yac = y(ia) - y(ic)
77            zac = z(ia) - z(ic)
78            if (use_polymer)  call image (xac,yac,zac)
79            rac = sqrt(xac*xac + yac*yac + zac*zac)
80            dt = rac - ideal
81            dt2 = dt * dt
82            e = ureyunit * force * dt2 * (1.0d0+cury*dt+qury*dt2)
83            deddt = 2.0d0 * ureyunit * force * dt
84     &                 * (1.0d0+1.5d0*cury*dt+2.0d0*qury*dt2)
85c
86c     scale the interaction based on its group membership
87c
88            if (use_group) then
89               e = e * fgrp
90               deddt = deddt * fgrp
91            end if
92c
93c     compute chain rule terms needed for derivatives
94c
95            de = deddt / rac
96            dedx = de * xac
97            dedy = de * yac
98            dedz = de * zac
99c
100c     increment the total Urey-Bradley energy and first derivatives
101c
102            eub = eub + e
103            deub(1,ia) = deub(1,ia) + dedx
104            deub(2,ia) = deub(2,ia) + dedy
105            deub(3,ia) = deub(3,ia) + dedz
106            deub(1,ic) = deub(1,ic) - dedx
107            deub(2,ic) = deub(2,ic) - dedy
108            deub(3,ic) = deub(3,ic) - dedz
109c
110c     increment the internal virial tensor components
111c
112            vxx = xac * dedx
113            vyx = yac * dedx
114            vzx = zac * dedx
115            vyy = yac * dedy
116            vzy = zac * dedy
117            vzz = zac * dedz
118            vir(1,1) = vir(1,1) + vxx
119            vir(2,1) = vir(2,1) + vyx
120            vir(3,1) = vir(3,1) + vzx
121            vir(1,2) = vir(1,2) + vyx
122            vir(2,2) = vir(2,2) + vyy
123            vir(3,2) = vir(3,2) + vzy
124            vir(1,3) = vir(1,3) + vzx
125            vir(2,3) = vir(2,3) + vzy
126            vir(3,3) = vir(3,3) + vzz
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