1c
2c
3c     ###################################################
4c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
5c     ##              All Rights Reserved              ##
6c     ###################################################
7c
8c     ################################################################
9c     ##                                                            ##
10c     ##  subroutine xyzatm  --  single atom internal to Cartesian  ##
11c     ##                                                            ##
12c     ################################################################
13c
14c
15c     "xyzatm" computes the Cartesian coordinates of a single
16c     atom from its defining internal coordinate values
17c
18c
19      subroutine xyzatm (i,ia,bond,ib,angle1,ic,angle2,chiral)
20      use atoms
21      use inform
22      use iounit
23      use math
24      implicit none
25      integer i,ia,ib,ic,chiral
26      real*8 bond,angle1,angle2
27      real*8 eps,rad1,rad2
28      real*8 sin1,cos1,sin2,cos2
29      real*8 cosine,sine,sine2
30      real*8 xab,yab,zab,rab
31      real*8 xba,yba,zba,rba
32      real*8 xbc,ybc,zbc,rbc
33      real*8 xac,yac,zac,rac
34      real*8 xt,yt,zt,xu,yu,zu
35      real*8 cosb,sinb,cosg,sing
36      real*8 xtmp,ztmp,a,b,c
37c
38c
39c     convert angles to radians, and get their sines and cosines
40c
41      eps = 0.00000001d0
42      rad1 = angle1 / radian
43      rad2 = angle2 / radian
44      sin1 = sin(rad1)
45      cos1 = cos(rad1)
46      sin2 = sin(rad2)
47      cos2 = cos(rad2)
48c
49c     if no second site given, place the atom at the origin
50c
51      if (ia .eq. 0) then
52         x(i) = 0.0d0
53         y(i) = 0.0d0
54         z(i) = 0.0d0
55c
56c     if no third site given, place the atom along the z-axis
57c
58      else if (ib .eq. 0) then
59         x(i) = x(ia)
60         y(i) = y(ia)
61         z(i) = z(ia) + bond
62c
63c     if no fourth site given, place the atom in the x,z-plane
64c
65      else if (ic .eq. 0) then
66         xab = x(ia) - x(ib)
67         yab = y(ia) - y(ib)
68         zab = z(ia) - z(ib)
69         rab = sqrt(xab**2 + yab**2 + zab**2)
70         xab = xab / rab
71         yab = yab / rab
72         zab = zab / rab
73         cosb = zab
74         sinb = sqrt(xab**2 + yab**2)
75         if (sinb .eq. 0.0d0) then
76            cosg = 1.0d0
77            sing = 0.0d0
78         else
79            cosg = yab / sinb
80            sing = xab / sinb
81         end if
82         xtmp = bond*sin1
83         ztmp = rab - bond*cos1
84         x(i) = x(ib) + xtmp*cosg + ztmp*sing*sinb
85         y(i) = y(ib) - xtmp*sing + ztmp*cosg*sinb
86         z(i) = z(ib) + ztmp*cosb
87c
88c     general case where the second angle is a dihedral angle
89c
90      else if (chiral .eq. 0) then
91         xab = x(ia) - x(ib)
92         yab = y(ia) - y(ib)
93         zab = z(ia) - z(ib)
94         rab = sqrt(xab**2 + yab**2 + zab**2)
95         xab = xab / rab
96         yab = yab / rab
97         zab = zab / rab
98         xbc = x(ib) - x(ic)
99         ybc = y(ib) - y(ic)
100         zbc = z(ib) - z(ic)
101         rbc = sqrt(xbc**2 + ybc**2 + zbc**2)
102         xbc = xbc / rbc
103         ybc = ybc / rbc
104         zbc = zbc / rbc
105         xt = zab*ybc - yab*zbc
106         yt = xab*zbc - zab*xbc
107         zt = yab*xbc - xab*ybc
108         cosine = xab*xbc + yab*ybc + zab*zbc
109         sine = sqrt(max(1.0d0-cosine**2,eps))
110         xt = xt / sine
111         yt = yt / sine
112         zt = zt / sine
113         xu = yt*zab - zt*yab
114         yu = zt*xab - xt*zab
115         zu = xt*yab - yt*xab
116         x(i) = x(ia) + bond * (xu*sin1*cos2 + xt*sin1*sin2 - xab*cos1)
117         y(i) = y(ia) + bond * (yu*sin1*cos2 + yt*sin1*sin2 - yab*cos1)
118         z(i) = z(ia) + bond * (zu*sin1*cos2 + zt*sin1*sin2 - zab*cos1)
119         if (abs(cosine) .ge. 1.0d0) then
120            cosb = zab
121            sinb = sqrt(xab**2 + yab**2)
122            if (sinb .eq. 0.0d0) then
123               cosg = 1.0d0
124               sing = 0.0d0
125            else
126               cosg = yab / sinb
127               sing = xab / sinb
128            end if
129            xtmp = bond*sin1
130            ztmp = rab - bond*cos1
131            x(i) = x(ib) + xtmp*cosg + ztmp*sing*sinb
132            y(i) = y(ib) - xtmp*sing + ztmp*cosg*sinb
133            z(i) = z(ib) + ztmp*cosb
134            write (iout,10)  i
135   10       format (/,' XYZATM  --  Warning, Undefined Dihedral',
136     &                 ' Angle at Atom',i6)
137         end if
138c
139c     general case where the second angle is a bond angle
140c
141      else if (abs(chiral) .eq. 1) then
142         xba = x(ib) - x(ia)
143         yba = y(ib) - y(ia)
144         zba = z(ib) - z(ia)
145         rba = sqrt(xba**2 + yba**2 + zba**2)
146         xba = xba / rba
147         yba = yba / rba
148         zba = zba / rba
149         xac = x(ia) - x(ic)
150         yac = y(ia) - y(ic)
151         zac = z(ia) - z(ic)
152         rac = sqrt(xac**2 + yac**2 + zac**2)
153         xac = xac / rac
154         yac = yac / rac
155         zac = zac / rac
156         xt = zba*yac - yba*zac
157         yt = xba*zac - zba*xac
158         zt = yba*xac - xba*yac
159         cosine = xba*xac + yba*yac + zba*zac
160         sine2 = max(1.0d0-cosine**2,eps)
161         if (abs(cosine) .ge. 1.0d0) then
162            write (iout,20)  i
163   20       format (/,' XYZATM  --  Warning, Collinear Defining',
164     &                 ' Atoms at Atom',i6)
165         end if
166         a = (-cos2 - cosine*cos1) / sine2
167         b = (cos1 + cosine*cos2) / sine2
168         c = (1.0d0 + a*cos2 - b*cos1) / sine2
169         if (c .gt. eps) then
170            c = chiral * sqrt(c)
171         else if (c .lt. -eps) then
172            c = sqrt((a*xac+b*xba)**2 + (a*yac+b*yba)**2
173     &                       + (a*zac+b*zba)**2)
174            a = a / c
175            b = b / c
176            c = 0.0d0
177            if (debug) then
178               write (iout,30)  ia
179   30          format (/,' XYZATM  --  Warning, Sum of Bond Angles',
180     &                    ' Too Large at Atom',i6)
181            end if
182         else
183            c = 0.0d0
184         end if
185         x(i) = x(ia) + bond * (a*xac + b*xba + c*xt)
186         y(i) = y(ia) + bond * (a*yac + b*yba + c*yt)
187         z(i) = z(ia) + bond * (a*zac + b*zba + c*zt)
188      end if
189      return
190      end
191