1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ################################################################
9c     ##                                                            ##
10c     ##  subroutine eurey2  --  atom-by-atom Urey-Bradley Hessian  ##
11c     ##                                                            ##
12c     ################################################################
13c
14c
15c     "eurey2" calculates second derivatives of the Urey-Bradley
16c     interaction energy for a single atom at a time
17c
18c
19      subroutine eurey2 (i)
20      use atoms
21      use bound
22      use couple
23      use group
24      use hessn
25      use urey
26      use urypot
27      implicit none
28      integer i,j,ia,ic,iurey
29      real*8 ideal,force,fgrp
30      real*8 xac,yac,zac
31      real*8 rac,rac2
32      real*8 dt,dt2,term
33      real*8 termx,termy,termz
34      real*8 de,deddt,d2eddt2
35      real*8 d2e(3,3)
36      logical proceed
37c
38c
39c     calculate the Urey-Bradley interaction Hessian elements
40c
41      do iurey = 1, nurey
42         ia = iury(1,iurey)
43         ic = iury(3,iurey)
44         ideal = ul(iurey)
45         force = uk(iurey)
46c
47c     decide whether to compute the current interaction
48c
49         proceed = (i.eq.ia .or. i.eq.ic)
50         if (proceed .and. use_group)
51     &      call groups (proceed,fgrp,ia,ic,0,0,0,0)
52c
53c     compute the value of the 1-3 distance deviation
54c
55         if (proceed) then
56            if (i .eq. ic) then
57               ic = ia
58               ia = i
59            end if
60            xac = x(ia) - x(ic)
61            yac = y(ia) - y(ic)
62            zac = z(ia) - z(ic)
63            if (use_polymer)  call image (xac,yac,zac)
64            rac2 = xac*xac + yac*yac + zac*zac
65            rac = sqrt(rac2)
66            dt = rac - ideal
67            dt2 = dt * dt
68            deddt = 2.0d0 * ureyunit * force * dt
69     &                 * (1.0d0+1.5d0*cury*dt+2.0d0*qury*dt2)
70            d2eddt2 = 2.0d0 * ureyunit * force
71     &                   * (1.0d0+3.0d0*cury*dt+6.0d0*qury*dt2)
72c
73c     scale the interaction based on its group membership
74c
75            if (use_group) then
76               deddt = deddt * fgrp
77               d2eddt2 = d2eddt2 * fgrp
78            end if
79c
80c     set the chain rule terms for the Hessian elements
81c
82            de = deddt / rac
83            term = (d2eddt2-de) / rac2
84            termx = term * xac
85            termy = term * yac
86            termz = term * zac
87            d2e(1,1) = termx*xac + de
88            d2e(1,2) = termx*yac
89            d2e(1,3) = termx*zac
90            d2e(2,1) = d2e(1,2)
91            d2e(2,2) = termy*yac + de
92            d2e(2,3) = termy*zac
93            d2e(3,1) = d2e(1,3)
94            d2e(3,2) = d2e(2,3)
95            d2e(3,3) = termz*zac + de
96c
97c     increment diagonal and non-diagonal Hessian elements
98c
99            do j = 1, 3
100               hessx(j,ia) = hessx(j,ia) + d2e(1,j)
101               hessy(j,ia) = hessy(j,ia) + d2e(2,j)
102               hessz(j,ia) = hessz(j,ia) + d2e(3,j)
103               hessx(j,ic) = hessx(j,ic) - d2e(1,j)
104               hessy(j,ic) = hessy(j,ic) - d2e(2,j)
105               hessz(j,ic) = hessz(j,ic) - d2e(3,j)
106            end do
107         end if
108      end do
109      return
110      end
111