--- /dev/null
+c
+c
+c ###################################################
+c ## COPYRIGHT (C) 1996 by Jay William Ponder ##
+c ## All Rights Reserved ##
+c ###################################################
+c
+c ################################################################
+c ## ##
+c ## subroutine surfatom -- exposed surface area of an atom ##
+c ## ##
+c ################################################################
+c
+c
+c "surfatom" performs an analytical computation of the surface
+c area of a specified atom; a simplified version of "surface"
+c
+c literature references:
+c
+c T. J. Richmond, "Solvent Accessible Surface Area and
+c Excluded Volume in Proteins", Journal of Molecular Biology,
+c 178, 63-89 (1984)
+c
+c L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
+c Applied to Molecular Dynamics of Proteins in Solution",
+c Protein Science, 1, 227-235 (1992)
+c
+c variables and parameters:
+c
+c ir number of atom for which area is desired
+c area accessible surface area of the atom
+c radius radii of each of the individual atoms
+c
+c
+ subroutine surfatom (ir,area,radius)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'sizes.i'
+ include 'COMMON.GEO'
+ include 'COMMON.IOUNITS'
+ integer nres,nsup,nstart_sup
+ double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
+ common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
+ & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
+ & dc_work(MAXRES6),nres,nres0
+ integer maxarc
+ parameter (maxarc=300)
+ integer i,j,k,m
+ integer ii,ib,jb
+ integer io,ir
+ integer mi,ni,narc
+ integer key(maxarc)
+ integer intag(maxarc)
+ integer intag1(maxarc)
+ real*8 area,arcsum
+ real*8 arclen,exang
+ real*8 delta,delta2
+ real*8 eps,rmove
+ real*8 xr,yr,zr
+ real*8 rr,rrsq
+ real*8 rplus,rminus
+ real*8 axx,axy,axz
+ real*8 ayx,ayy
+ real*8 azx,azy,azz
+ real*8 uxj,uyj,uzj
+ real*8 tx,ty,tz
+ real*8 txb,tyb,td
+ real*8 tr2,tr,txr,tyr
+ real*8 tk1,tk2
+ real*8 thec,the,t,tb
+ real*8 txk,tyk,tzk
+ real*8 t1,ti,tf,tt
+ real*8 txj,tyj,tzj
+ real*8 ccsq,cc,xysq
+ real*8 bsqk,bk,cosine
+ real*8 dsqj,gi,pix2
+ real*8 therk,dk,gk
+ real*8 risqk,rik
+ real*8 radius(maxatm)
+ real*8 ri(maxarc),risq(maxarc)
+ real*8 ux(maxarc),uy(maxarc),uz(maxarc)
+ real*8 xc(maxarc),yc(maxarc),zc(maxarc)
+ real*8 xc1(maxarc),yc1(maxarc),zc1(maxarc)
+ real*8 dsq(maxarc),bsq(maxarc)
+ real*8 dsq1(maxarc),bsq1(maxarc)
+ real*8 arci(maxarc),arcf(maxarc)
+ real*8 ex(maxarc),lt(maxarc),gr(maxarc)
+ real*8 b(maxarc),b1(maxarc),bg(maxarc)
+ real*8 kent(maxarc),kout(maxarc)
+ real*8 ther(maxarc)
+ logical moved,top
+ logical omit(maxarc)
+c
+c
+c zero out the surface area for the sphere of interest
+c
+ area = 0.0d0
+c write (2,*) "ir",ir," radius",radius(ir)
+ if (radius(ir) .eq. 0.0d0) return
+c
+c set the overlap significance and connectivity shift
+c
+ pix2 = 2.0d0 * pi
+ delta = 1.0d-8
+ delta2 = delta * delta
+ eps = 1.0d-8
+ moved = .false.
+ rmove = 1.0d-8
+c
+c store coordinates and radius of the sphere of interest
+c
+ xr = c(1,ir)
+ yr = c(2,ir)
+ zr = c(3,ir)
+ rr = radius(ir)
+ rrsq = rr * rr
+c
+c initialize values of some counters and summations
+c
+ 10 continue
+ io = 0
+ jb = 0
+ ib = 0
+ arclen = 0.0d0
+ exang = 0.0d0
+c
+c test each sphere to see if it overlaps the sphere of interest
+c
+ do i = 1, 2*nres
+ if (i.eq.ir .or. radius(i).eq.0.0d0) goto 30
+ rplus = rr + radius(i)
+ tx = c(1,i) - xr
+ if (abs(tx) .ge. rplus) goto 30
+ ty = c(2,i) - yr
+ if (abs(ty) .ge. rplus) goto 30
+ tz = c(3,i) - zr
+ if (abs(tz) .ge. rplus) goto 30
+c
+c check for sphere overlap by testing distance against radii
+c
+ xysq = tx*tx + ty*ty
+ if (xysq .lt. delta2) then
+ tx = delta
+ ty = 0.0d0
+ xysq = delta2
+ end if
+ ccsq = xysq + tz*tz
+ cc = sqrt(ccsq)
+ if (rplus-cc .le. delta) goto 30
+ rminus = rr - radius(i)
+c
+c check to see if sphere of interest is completely buried
+c
+ if (cc-abs(rminus) .le. delta) then
+ if (rminus .le. 0.0d0) goto 170
+ goto 30
+ end if
+c
+c check for too many overlaps with sphere of interest
+c
+ if (io .ge. maxarc) then
+ write (iout,20)
+ 20 format (/,' SURFATOM -- Increase the Value of MAXARC')
+ stop
+ end if
+c
+c get overlap between current sphere and sphere of interest
+c
+ io = io + 1
+ xc1(io) = tx
+ yc1(io) = ty
+ zc1(io) = tz
+ dsq1(io) = xysq
+ bsq1(io) = ccsq
+ b1(io) = cc
+ gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
+ intag1(io) = i
+ omit(io) = .false.
+ 30 continue
+ end do
+c
+c case where no other spheres overlap the sphere of interest
+c
+ if (io .eq. 0) then
+ area = 4.0d0 * pi * rrsq
+ return
+ end if
+c
+c case where only one sphere overlaps the sphere of interest
+c
+ if (io .eq. 1) then
+ area = pix2 * (1.0d0 + gr(1))
+ area = mod(area,4.0d0*pi) * rrsq
+ return
+ end if
+c
+c case where many spheres intersect the sphere of interest;
+c sort the intersecting spheres by their degree of overlap
+c
+ call sort2 (io,gr,key)
+ do i = 1, io
+ k = key(i)
+ intag(i) = intag1(k)
+ xc(i) = xc1(k)
+ yc(i) = yc1(k)
+ zc(i) = zc1(k)
+ dsq(i) = dsq1(k)
+ b(i) = b1(k)
+ bsq(i) = bsq1(k)
+ end do
+c
+c get radius of each overlap circle on surface of the sphere
+c
+ do i = 1, io
+ gi = gr(i) * rr
+ bg(i) = b(i) * gi
+ risq(i) = rrsq - gi*gi
+ ri(i) = sqrt(risq(i))
+ ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
+ end do
+c
+c find boundary of inaccessible area on sphere of interest
+c
+ do k = 1, io-1
+ if (.not. omit(k)) then
+ txk = xc(k)
+ tyk = yc(k)
+ tzk = zc(k)
+ bk = b(k)
+ therk = ther(k)
+c
+c check to see if J circle is intersecting K circle;
+c get distance between circle centers and sum of radii
+c
+ do j = k+1, io
+ if (omit(j)) goto 60
+ cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
+ cc = acos(min(1.0d0,max(-1.0d0,cc)))
+ td = therk + ther(j)
+c
+c check to see if circles enclose separate regions
+c
+ if (cc .ge. td) goto 60
+c
+c check for circle J completely inside circle K
+c
+ if (cc+ther(j) .lt. therk) goto 40
+c
+c check for circles that are essentially parallel
+c
+ if (cc .gt. delta) goto 50
+ 40 continue
+ omit(j) = .true.
+ goto 60
+c
+c check to see if sphere of interest is completely buried
+c
+ 50 continue
+ if (pix2-cc .le. td) goto 170
+ 60 continue
+ end do
+ end if
+ end do
+c
+c find T value of circle intersections
+c
+ do k = 1, io
+ if (omit(k)) goto 110
+ omit(k) = .true.
+ narc = 0
+ top = .false.
+ txk = xc(k)
+ tyk = yc(k)
+ tzk = zc(k)
+ dk = sqrt(dsq(k))
+ bsqk = bsq(k)
+ bk = b(k)
+ gk = gr(k) * rr
+ risqk = risq(k)
+ rik = ri(k)
+ therk = ther(k)
+c
+c rotation matrix elements
+c
+ t1 = tzk / (bk*dk)
+ axx = txk * t1
+ axy = tyk * t1
+ axz = dk / bk
+ ayx = tyk / dk
+ ayy = txk / dk
+ azx = txk / bk
+ azy = tyk / bk
+ azz = tzk / bk
+ do j = 1, io
+ if (.not. omit(j)) then
+ txj = xc(j)
+ tyj = yc(j)
+ tzj = zc(j)
+c
+c rotate spheres so K vector colinear with z-axis
+c
+ uxj = txj*axx + tyj*axy - tzj*axz
+ uyj = tyj*ayy - txj*ayx
+ uzj = txj*azx + tyj*azy + tzj*azz
+ cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
+ if (acos(cosine) .lt. therk+ther(j)) then
+ dsqj = uxj*uxj + uyj*uyj
+ tb = uzj*gk - bg(j)
+ txb = uxj * tb
+ tyb = uyj * tb
+ td = rik * dsqj
+ tr2 = risqk*dsqj - tb*tb
+ tr2 = max(eps,tr2)
+ tr = sqrt(tr2)
+ txr = uxj * tr
+ tyr = uyj * tr
+c
+c get T values of intersection for K circle
+c
+ tb = (txb+tyr) / td
+ tb = min(1.0d0,max(-1.0d0,tb))
+ tk1 = acos(tb)
+ if (tyb-txr .lt. 0.0d0) tk1 = pix2 - tk1
+ tb = (txb-tyr) / td
+ tb = min(1.0d0,max(-1.0d0,tb))
+ tk2 = acos(tb)
+ if (tyb+txr .lt. 0.0d0) tk2 = pix2 - tk2
+ thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
+ if (abs(thec) .lt. 1.0d0) then
+ the = -acos(thec)
+ else if (thec .ge. 1.0d0) then
+ the = 0.0d0
+ else if (thec .le. -1.0d0) then
+ the = -pi
+ end if
+c
+c see if "tk1" is entry or exit point; check t=0 point;
+c "ti" is exit point, "tf" is entry point
+c
+ cosine = min(1.0d0,max(-1.0d0,
+ & (uzj*gk-uxj*rik)/(b(j)*rr)))
+ if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
+ ti = tk2
+ tf = tk1
+ else
+ ti = tk2
+ tf = tk1
+ end if
+ narc = narc + 1
+ if (narc .ge. maxarc) then
+ write (iout,70)
+ 70 format (/,' SURFATOM -- Increase the Value',
+ & ' of MAXARC')
+ stop
+ end if
+ if (tf .le. ti) then
+ arcf(narc) = tf
+ arci(narc) = 0.0d0
+ tf = pix2
+ lt(narc) = j
+ ex(narc) = the
+ top = .true.
+ narc = narc + 1
+ end if
+ arcf(narc) = tf
+ arci(narc) = ti
+ lt(narc) = j
+ ex(narc) = the
+ ux(j) = uxj
+ uy(j) = uyj
+ uz(j) = uzj
+ end if
+ end if
+ end do
+ omit(k) = .false.
+c
+c special case; K circle without intersections
+c
+ if (narc .le. 0) goto 90
+c
+c general case; sum up arclength and set connectivity code
+c
+ call sort2 (narc,arci,key)
+ arcsum = arci(1)
+ mi = key(1)
+ t = arcf(mi)
+ ni = mi
+ if (narc .gt. 1) then
+ do j = 2, narc
+ m = key(j)
+ if (t .lt. arci(j)) then
+ arcsum = arcsum + arci(j) - t
+ exang = exang + ex(ni)
+ jb = jb + 1
+ if (jb .ge. maxarc) then
+ write (iout,80)
+ 80 format (/,' SURFATOM -- Increase the Value',
+ & ' of MAXARC')
+ stop
+ end if
+ i = lt(ni)
+ kent(jb) = maxarc*i + k
+ i = lt(m)
+ kout(jb) = maxarc*k + i
+ end if
+ tt = arcf(m)
+ if (tt .ge. t) then
+ t = tt
+ ni = m
+ end if
+ end do
+ end if
+ arcsum = arcsum + pix2 - t
+ if (.not. top) then
+ exang = exang + ex(ni)
+ jb = jb + 1
+ i = lt(ni)
+ kent(jb) = maxarc*i + k
+ i = lt(mi)
+ kout(jb) = maxarc*k + i
+ end if
+ goto 100
+ 90 continue
+ arcsum = pix2
+ ib = ib + 1
+ 100 continue
+ arclen = arclen + gr(k)*arcsum
+ 110 continue
+ end do
+ if (arclen .eq. 0.0d0) goto 170
+ if (jb .eq. 0) goto 150
+c
+c find number of independent boundaries and check connectivity
+c
+ j = 0
+ do k = 1, jb
+ if (kout(k) .ne. 0) then
+ i = k
+ 120 continue
+ m = kout(i)
+ kout(i) = 0
+ j = j + 1
+ do ii = 1, jb
+ if (m .eq. kent(ii)) then
+ if (ii .eq. k) then
+ ib = ib + 1
+ if (j .eq. jb) goto 150
+ goto 130
+ end if
+ i = ii
+ goto 120
+ end if
+ end do
+ 130 continue
+ end if
+ end do
+ ib = ib + 1
+c
+c attempt to fix connectivity error by moving atom slightly
+c
+ if (moved) then
+ write (iout,140) ir
+ 140 format (/,' SURFATOM -- Connectivity Error at Atom',i6)
+ else
+ moved = .true.
+ xr = xr + rmove
+ yr = yr + rmove
+ zr = zr + rmove
+ goto 10
+ end if
+c
+c compute the exposed surface area for the sphere of interest
+c
+ 150 continue
+ area = ib*pix2 + exang + arclen
+ area = mod(area,4.0d0*pi) * rrsq
+c
+c attempt to fix negative area by moving atom slightly
+c
+ if (area .lt. 0.0d0) then
+ if (moved) then
+ write (iout,160) ir
+ 160 format (/,' SURFATOM -- Negative Area at Atom',i6)
+ else
+ moved = .true.
+ xr = xr + rmove
+ yr = yr + rmove
+ zr = zr + rmove
+ goto 10
+ end if
+ end if
+ 170 continue
+ return
+ end
projects / unres.git / blobdiff