1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1999  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ################################################################
9c     ##                                                            ##
10c     ##  function geometry  --  evaluate distance, angle, torsion  ##
11c     ##                                                            ##
12c     ################################################################
13c
14c
15c     "geometry" finds the value of the interatomic distance, angle
16c     or dihedral angle defined by two to four input atoms
17c
18c
19      function geometry (ia,ib,ic,id)
20      use atoms
21      use math
22      implicit none
23      integer ia,ib,ic,id
24      real*8 xab,yab,zab
25      real*8 xba,yba,zba
26      real*8 xcb,ycb,zcb
27      real*8 xdc,ydc,zdc
28      real*8 xt,yt,zt
29      real*8 xu,yu,zu
30      real*8 rab2,rcb2,rabc
31      real*8 rt2,ru2,rtru
32      real*8 cosine,sign
33      real*8 geometry
34c
35c
36c     set default in case atoms are coincident or colinear
37c
38      geometry = 0.0d0
39c
40c     compute the value of the distance in angstroms
41c
42      if (ic .eq. 0) then
43         xab = x(ia) - x(ib)
44         yab = y(ia) - y(ib)
45         zab = z(ia) - z(ib)
46         geometry = sqrt(xab*xab + yab*yab + zab*zab)
47c
48c     compute the value of the angle in degrees
49c
50      else if (id .eq. 0) then
51         xab = x(ia) - x(ib)
52         yab = y(ia) - y(ib)
53         zab = z(ia) - z(ib)
54         xcb = x(ic) - x(ib)
55         ycb = y(ic) - y(ib)
56         zcb = z(ic) - z(ib)
57         rab2 = max(xab*xab+yab*yab+zab*zab,0.0001d0)
58         rcb2 = max(xcb*xcb+ycb*ycb+zcb*zcb,0.0001d0)
59         rabc = sqrt(rab2*rcb2)
60         cosine = (xab*xcb + yab*ycb + zab*zcb) / rabc
61         cosine = min(1.0d0,max(-1.0d0,cosine))
62         geometry = radian * acos(cosine)
63c
64c     compute the value of the dihedral angle in degrees
65c
66      else
67         xba = x(ib) - x(ia)
68         yba = y(ib) - y(ia)
69         zba = z(ib) - z(ia)
70         xcb = x(ic) - x(ib)
71         ycb = y(ic) - y(ib)
72         zcb = z(ic) - z(ib)
73         xdc = x(id) - x(ic)
74         ydc = y(id) - y(ic)
75         zdc = z(id) - z(ic)
76         xt = yba*zcb - ycb*zba
77         yt = xcb*zba - xba*zcb
78         zt = xba*ycb - xcb*yba
79         xu = ycb*zdc - ydc*zcb
80         yu = xdc*zcb - xcb*zdc
81         zu = xcb*ydc - xdc*ycb
82         rt2 = max(xt*xt+yt*yt+zt*zt,0.0001d0)
83         ru2 = max(xu*xu+yu*yu+zu*zu,0.0001d0)
84         rtru = sqrt(rt2*ru2)
85         cosine = (xt*xu + yt*yu + zt*zu) / rtru
86         cosine = min(1.0d0,max(-1.0d0,cosine))
87         geometry = radian * acos(cosine)
88         sign = xba*xu + yba*yu + zba*zu
89         if (sign .lt. 0.0d0)  geometry = -geometry
90      end if
91      return
92      end
93