1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     #############################################################
9c     ##                                                         ##
10c     ##  subroutine eurey3  --  Urey-Bradley energy & analysis  ##
11c     ##                                                         ##
12c     #############################################################
13c
14c
15c     "eurey3" calculates the Urey-Bradley energy; also
16c     partitions the energy among the atoms
17c
18c
19      subroutine eurey3
20      use action
21      use analyz
22      use atomid
23      use atoms
24      use bound
25      use energi
26      use group
27      use inform
28      use iounit
29      use urey
30      use urypot
31      use usage
32      implicit none
33      integer i,ia,ib,ic
34      real*8 e,ideal,force
35      real*8 dt,dt2,fgrp
36      real*8 xac,yac,zac,rac
37      logical proceed
38      logical header,huge
39c
40c
41c     zero out the Urey-Bradley energy and partitioning terms
42c
43      neub = 0
44      eub = 0.0d0
45      do i = 1, n
46         aeub(i) = 0.0d0
47      end do
48      if (nurey .eq. 0)  return
49c
50c     print header information if debug output was requested
51c
52      header = .true.
53      if (debug .and. nurey.ne.0) then
54         header = .false.
55         write (iout,10)
56   10    format (/,' Individual Urey-Bradley Interactions :',
57     &           //,' Type',18x,'Atom Names',18x,'Ideal',
58     &              4x,'Actual',6x,'Energy',/)
59      end if
60c
61c     OpenMP directives for the major loop structure
62c
63!$OMP PARALLEL default(private) shared(nurey,iury,ul,uk,
64!$OMP& use,x,y,z,cury,qury,ureyunit,use_group,use_polymer,
65!$OMP& name,verbose,debug,header,iout)
66!$OMP& shared(eub,neub,aeub)
67!$OMP DO reduction(+:eub,neub,aeub) schedule(guided)
68c
69c     calculate the Urey-Bradley 1-3 energy term
70c
71      do i = 1, nurey
72         ia = iury(1,i)
73         ib = iury(2,i)
74         ic = iury(3,i)
75         ideal = ul(i)
76         force = uk(i)
77c
78c     decide whether to compute the current interaction
79c
80         proceed = .true.
81         if (use_group)  call groups (proceed,fgrp,ia,ic,0,0,0,0)
82         if (proceed)  proceed = (use(ia) .or. use(ic))
83c
84c     compute the value of the 1-3 distance deviation
85c
86         if (proceed) then
87            xac = x(ia) - x(ic)
88            yac = y(ia) - y(ic)
89            zac = z(ia) - z(ic)
90            if (use_polymer)  call image (xac,yac,zac)
91            rac = sqrt(xac*xac + yac*yac + zac*zac)
92            dt = rac - ideal
93            dt2 = dt * dt
94c
95c     calculate the Urey-Bradley energy for this interaction
96c
97            e = ureyunit * force * dt2 * (1.0d0+cury*dt+qury*dt2)
98c
99c     scale the interaction based on its group membership
100c
101            if (use_group)  e = e * fgrp
102c
103c     increment the total Urey-Bradley energy
104c
105            neub = neub + 1
106            eub = eub + e
107            aeub(ia) = aeub(ia) + 0.5d0*e
108            aeub(ic) = aeub(ic) + 0.5d0*e
109c
110c     print a message if the energy of this interaction is large
111c
112            huge = (e .gt. 5.0d0)
113            if (debug .or. (verbose.and.huge)) then
114               if (header) then
115                  header = .false.
116                  write (iout,20)
117   20             format (/,' Individual Urey-Bradley Interactions :',
118     &                    //,' Type',18x,'Atom Names',18x,'Ideal',
119     &                       4x,'Actual',6x,'Energy',/)
120               end if
121               write (iout,30)  ia,name(ia),ib,name(ib),
122     &                          ic,name(ic),ideal,rac,e
123   30          format (' UreyBrad',2x,i7,'-',a3,i7,'-',a3,
124     &                    i7,'-',a3,2x,2f10.4,f12.4)
125            end if
126         end if
127      end do
128c
129c     OpenMP directives for the major loop structure
130c
131!$OMP END DO
132!$OMP END PARALLEL
133      return
134      end
135