+++ /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