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