X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fsrc_MD_DFA%2Fsurfatom.f;fp=source%2Funres%2Fsrc_MD_DFA%2Fsurfatom.f;h=99748421def07e9b8d92c7571c73d9b38cdebcdc;hb=d21d93a57773de0c363db8764a376202f3a0552d;hp=0000000000000000000000000000000000000000;hpb=b8f5b4118dc932c21ad2bdbba905850f8297a6a2;p=unres.git

diff --git a/source/unres/src_MD_DFA/surfatom.f b/source/unres/src_MD_DFA/surfatom.f
new file mode 100644
index 0000000..9974842
--- /dev/null
+++ b/source/unres/src_MD_DFA/surfatom.f
@@ -0,0 +1,494 @@
+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