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