source/unres/src_MD/surfatom.f

   1 c
   2 c
   3 c     ###################################################
   4 c     ##  COPYRIGHT (C)  1996  by  Jay William Ponder  ##
   5 c     ##              All Rights Reserved              ##
   6 c     ###################################################
   7 c
   8 c     ################################################################
   9 c     ##                                                            ##
  10 c     ##  subroutine surfatom  --  exposed surface area of an atom  ##
  11 c     ##                                                            ##
  12 c     ################################################################
  13 c
  14 c
  15 c     "surfatom" performs an analytical computation of the surface
  16 c     area of a specified atom; a simplified version of "surface"
  17 c
  18 c     literature references:
  19 c
  20 c     T. J. Richmond, "Solvent Accessible Surface Area and
  21 c     Excluded Volume in Proteins", Journal of Molecular Biology,
  22 c     178, 63-89 (1984)
  23 c
  24 c     L. Wesson and D. Eisenberg, "Atomic Solvation Parameters
  25 c     Applied to Molecular Dynamics of Proteins in Solution",
  26 c     Protein Science, 1, 227-235 (1992)
  27 c
  28 c     variables and parameters:
  29 c
  30 c     ir       number of atom for which area is desired
  31 c     area     accessible surface area of the atom
  32 c     radius   radii of each of the individual atoms
  33 c
  34 c
  35       subroutine surfatom (ir,area,radius)
  36       implicit real*8 (a-h,o-z)
  37       include 'DIMENSIONS'
  38       include 'sizes.i'
  39       include 'COMMON.GEO'
  40       include 'COMMON.IOUNITS'
  41       integer nres,nsup,nstart_sup
  42       double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm
  43       common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2),
  44      & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2),
  45      & dc_work(MAXRES6),nres,nres0
  46       integer maxarc
  47       parameter (maxarc=300)
  48       integer i,j,k,m
  49       integer ii,ib,jb
  50       integer io,ir
  51       integer mi,ni,narc
  52       integer key(maxarc)
  53       integer intag(maxarc)
  54       integer intag1(maxarc)
  55       real*8 area,arcsum
  56       real*8 arclen,exang
  57       real*8 delta,delta2
  58       real*8 eps,rmove
  59       real*8 xr,yr,zr
  60       real*8 rr,rrsq
  61       real*8 rplus,rminus
  62       real*8 axx,axy,axz
  63       real*8 ayx,ayy
  64       real*8 azx,azy,azz
  65       real*8 uxj,uyj,uzj
  66       real*8 tx,ty,tz
  67       real*8 txb,tyb,td
  68       real*8 tr2,tr,txr,tyr
  69       real*8 tk1,tk2
  70       real*8 thec,the,t,tb
  71       real*8 txk,tyk,tzk
  72       real*8 t1,ti,tf,tt
  73       real*8 txj,tyj,tzj
  74       real*8 ccsq,cc,xysq
  75       real*8 bsqk,bk,cosine
  76       real*8 dsqj,gi,pix2
  77       real*8 therk,dk,gk
  78       real*8 risqk,rik
  79       real*8 radius(maxatm)
  80       real*8 ri(maxarc),risq(maxarc)
  81       real*8 ux(maxarc),uy(maxarc),uz(maxarc)
  82       real*8 xc(maxarc),yc(maxarc),zc(maxarc)
  83       real*8 xc1(maxarc),yc1(maxarc),zc1(maxarc)
  84       real*8 dsq(maxarc),bsq(maxarc)
  85       real*8 dsq1(maxarc),bsq1(maxarc)
  86       real*8 arci(maxarc),arcf(maxarc)
  87       real*8 ex(maxarc),lt(maxarc),gr(maxarc)
  88       real*8 b(maxarc),b1(maxarc),bg(maxarc)
  89       real*8 kent(maxarc),kout(maxarc)
  90       real*8 ther(maxarc)
  91       logical moved,top
  92       logical omit(maxarc)
  93 c
  94 c
  95 c     zero out the surface area for the sphere of interest
  96 c
  97       area = 0.0d0
  98 c      write (2,*) "ir",ir," radius",radius(ir)
  99       if (radius(ir) .eq. 0.0d0)  return
 100 c
 101 c     set the overlap significance and connectivity shift
 102 c
 103       pix2 = 2.0d0 * pi
 104       delta = 1.0d-8
 105       delta2 = delta * delta
 106       eps = 1.0d-8
 107       moved = .false.
 108       rmove = 1.0d-8
 109 c
 110 c     store coordinates and radius of the sphere of interest
 111 c
 112       xr = c(1,ir)
 113       yr = c(2,ir)
 114       zr = c(3,ir)
 115       rr = radius(ir)
 116       rrsq = rr * rr
 117 c
 118 c     initialize values of some counters and summations
 119 c
 120    10 continue
 121       io = 0
 122       jb = 0
 123       ib = 0
 124       arclen = 0.0d0
 125       exang = 0.0d0
 126 c
 127 c     test each sphere to see if it overlaps the sphere of interest
 128 c
 129       do i = 1, 2*nres
 130          if (i.eq.ir .or. radius(i).eq.0.0d0)  goto 30
 131          rplus = rr + radius(i)
 132          tx = c(1,i) - xr
 133          if (abs(tx) .ge. rplus)  goto 30
 134          ty = c(2,i) - yr
 135          if (abs(ty) .ge. rplus)  goto 30
 136          tz = c(3,i) - zr
 137          if (abs(tz) .ge. rplus)  goto 30
 138 c
 139 c     check for sphere overlap by testing distance against radii
 140 c
 141          xysq = tx*tx + ty*ty
 142          if (xysq .lt. delta2) then
 143             tx = delta
 144             ty = 0.0d0
 145             xysq = delta2
 146          end if
 147          ccsq = xysq + tz*tz
 148          cc = sqrt(ccsq)
 149          if (rplus-cc .le. delta)  goto 30
 150          rminus = rr - radius(i)
 151 c
 152 c     check to see if sphere of interest is completely buried
 153 c
 154          if (cc-abs(rminus) .le. delta) then
 155             if (rminus .le. 0.0d0)  goto 170
 156             goto 30
 157          end if
 158 c
 159 c     check for too many overlaps with sphere of interest
 160 c
 161          if (io .ge. maxarc) then
 162             write (iout,20)
 163    20       format (/,' SURFATOM  --  Increase the Value of MAXARC')
 164             stop
 165          end if
 166 c
 167 c     get overlap between current sphere and sphere of interest
 168 c
 169          io = io + 1
 170          xc1(io) = tx
 171          yc1(io) = ty
 172          zc1(io) = tz
 173          dsq1(io) = xysq
 174          bsq1(io) = ccsq
 175          b1(io) = cc
 176          gr(io) = (ccsq+rplus*rminus) / (2.0d0*rr*b1(io))
 177          intag1(io) = i
 178          omit(io) = .false.
 179    30    continue
 180       end do
 181 c
 182 c     case where no other spheres overlap the sphere of interest
 183 c
 184       if (io .eq. 0) then
 185          area = 4.0d0 * pi * rrsq
 186          return
 187       end if
 188 c
 189 c     case where only one sphere overlaps the sphere of interest
 190 c
 191       if (io .eq. 1) then
 192          area = pix2 * (1.0d0 + gr(1))
 193          area = mod(area,4.0d0*pi) * rrsq
 194          return
 195       end if
 196 c
 197 c     case where many spheres intersect the sphere of interest;
 198 c     sort the intersecting spheres by their degree of overlap
 199 c
 200       call sort2 (io,gr,key)
 201       do i = 1, io
 202          k = key(i)
 203          intag(i) = intag1(k)
 204          xc(i) = xc1(k)
 205          yc(i) = yc1(k)
 206          zc(i) = zc1(k)
 207          dsq(i) = dsq1(k)
 208          b(i) = b1(k)
 209          bsq(i) = bsq1(k)
 210       end do
 211 c
 212 c     get radius of each overlap circle on surface of the sphere
 213 c
 214       do i = 1, io
 215          gi = gr(i) * rr
 216          bg(i) = b(i) * gi
 217          risq(i) = rrsq - gi*gi
 218          ri(i) = sqrt(risq(i))
 219          ther(i) = 0.5d0*pi - asin(min(1.0d0,max(-1.0d0,gr(i))))
 220       end do
 221 c
 222 c     find boundary of inaccessible area on sphere of interest
 223 c
 224       do k = 1, io-1
 225          if (.not. omit(k)) then
 226             txk = xc(k)
 227             tyk = yc(k)
 228             tzk = zc(k)
 229             bk = b(k)
 230             therk = ther(k)
 231 c
 232 c     check to see if J circle is intersecting K circle;
 233 c     get distance between circle centers and sum of radii
 234 c
 235             do j = k+1, io
 236                if (omit(j))  goto 60
 237                cc = (txk*xc(j)+tyk*yc(j)+tzk*zc(j))/(bk*b(j))
 238                cc = acos(min(1.0d0,max(-1.0d0,cc)))
 239                td = therk + ther(j)
 240 c
 241 c     check to see if circles enclose separate regions
 242 c
 243                if (cc .ge. td)  goto 60
 244 c
 245 c     check for circle J completely inside circle K
 246 c
 247                if (cc+ther(j) .lt. therk)  goto 40
 248 c
 249 c     check for circles that are essentially parallel
 250 c
 251                if (cc .gt. delta)  goto 50
 252    40          continue
 253                omit(j) = .true.
 254                goto 60
 255 c
 256 c     check to see if sphere of interest is completely buried
 257 c
 258    50          continue
 259                if (pix2-cc .le. td)  goto 170
 260    60          continue
 261             end do
 262          end if
 263       end do
 264 c
 265 c     find T value of circle intersections
 266 c
 267       do k = 1, io
 268          if (omit(k))  goto 110
 269          omit(k) = .true.
 270          narc = 0
 271          top = .false.
 272          txk = xc(k)
 273          tyk = yc(k)
 274          tzk = zc(k)
 275          dk = sqrt(dsq(k))
 276          bsqk = bsq(k)
 277          bk = b(k)
 278          gk = gr(k) * rr
 279          risqk = risq(k)
 280          rik = ri(k)
 281          therk = ther(k)
 282 c
 283 c     rotation matrix elements
 284 c
 285          t1 = tzk / (bk*dk)
 286          axx = txk * t1
 287          axy = tyk * t1
 288          axz = dk / bk
 289          ayx = tyk / dk
 290          ayy = txk / dk
 291          azx = txk / bk
 292          azy = tyk / bk
 293          azz = tzk / bk
 294          do j = 1, io
 295             if (.not. omit(j)) then
 296                txj = xc(j)
 297                tyj = yc(j)
 298                tzj = zc(j)
 299 c
 300 c     rotate spheres so K vector colinear with z-axis
 301 c
 302                uxj = txj*axx + tyj*axy - tzj*axz
 303                uyj = tyj*ayy - txj*ayx
 304                uzj = txj*azx + tyj*azy + tzj*azz
 305                cosine = min(1.0d0,max(-1.0d0,uzj/b(j)))
 306                if (acos(cosine) .lt. therk+ther(j)) then
 307                   dsqj = uxj*uxj + uyj*uyj
 308                   tb = uzj*gk - bg(j)
 309                   txb = uxj * tb
 310                   tyb = uyj * tb
 311                   td = rik * dsqj
 312                   tr2 = risqk*dsqj - tb*tb
 313                   tr2 = max(eps,tr2)
 314                   tr = sqrt(tr2)
 315                   txr = uxj * tr
 316                   tyr = uyj * tr
 317 c
 318 c     get T values of intersection for K circle
 319 c
 320                   tb = (txb+tyr) / td
 321                   tb = min(1.0d0,max(-1.0d0,tb))
 322                   tk1 = acos(tb)
 323                   if (tyb-txr .lt. 0.0d0)  tk1 = pix2 - tk1
 324                   tb = (txb-tyr) / td
 325                   tb = min(1.0d0,max(-1.0d0,tb))
 326                   tk2 = acos(tb)
 327                   if (tyb+txr .lt. 0.0d0)  tk2 = pix2 - tk2
 328                   thec = (rrsq*uzj-gk*bg(j)) / (rik*ri(j)*b(j))
 329                   if (abs(thec) .lt. 1.0d0) then
 330                      the = -acos(thec)
 331                   else if (thec .ge. 1.0d0) then
 332                      the = 0.0d0
 333                   else if (thec .le. -1.0d0) then
 334                      the = -pi
 335                   end if
 336 c
 337 c     see if "tk1" is entry or exit point; check t=0 point;
 338 c     "ti" is exit point, "tf" is entry point
 339 c
 340                   cosine = min(1.0d0,max(-1.0d0,
 341      &                            (uzj*gk-uxj*rik)/(b(j)*rr)))
 342                   if ((acos(cosine)-ther(j))*(tk2-tk1) .le. 0.0d0) then
 343                      ti = tk2
 344                      tf = tk1
 345                   else
 346                      ti = tk2
 347                      tf = tk1
 348                   end if
 349                   narc = narc + 1
 350                   if (narc .ge. maxarc) then
 351                      write (iout,70)
 352    70                format (/,' SURFATOM  --  Increase the Value',
 353      &                          ' of MAXARC')
 354                      stop
 355                   end if
 356                   if (tf .le. ti) then
 357                      arcf(narc) = tf
 358                      arci(narc) = 0.0d0
 359                      tf = pix2
 360                      lt(narc) = j
 361                      ex(narc) = the
 362                      top = .true.
 363                      narc = narc + 1
 364                   end if
 365                   arcf(narc) = tf
 366                   arci(narc) = ti
 367                   lt(narc) = j
 368                   ex(narc) = the
 369                   ux(j) = uxj
 370                   uy(j) = uyj
 371                   uz(j) = uzj
 372                end if
 373             end if
 374          end do
 375          omit(k) = .false.
 376 c
 377 c     special case; K circle without intersections
 378 c
 379          if (narc .le. 0)  goto 90
 380 c
 381 c     general case; sum up arclength and set connectivity code
 382 c
 383          call sort2 (narc,arci,key)
 384          arcsum = arci(1)
 385          mi = key(1)
 386          t = arcf(mi)
 387          ni = mi
 388          if (narc .gt. 1) then
 389             do j = 2, narc
 390                m = key(j)
 391                if (t .lt. arci(j)) then
 392                   arcsum = arcsum + arci(j) - t
 393                   exang = exang + ex(ni)
 394                   jb = jb + 1
 395                   if (jb .ge. maxarc) then
 396                      write (iout,80)
 397    80                format (/,' SURFATOM  --  Increase the Value',
 398      &                          ' of MAXARC')
 399                      stop
 400                   end if
 401                   i = lt(ni)
 402                   kent(jb) = maxarc*i + k
 403                   i = lt(m)
 404                   kout(jb) = maxarc*k + i
 405                end if
 406                tt = arcf(m)
 407                if (tt .ge. t) then
 408                   t = tt
 409                   ni = m
 410                end if
 411             end do
 412          end if
 413          arcsum = arcsum + pix2 - t
 414          if (.not. top) then
 415             exang = exang + ex(ni)
 416             jb = jb + 1
 417             i = lt(ni)
 418             kent(jb) = maxarc*i + k
 419             i = lt(mi)
 420             kout(jb) = maxarc*k + i
 421          end if
 422          goto 100
 423    90    continue
 424          arcsum = pix2
 425          ib = ib + 1
 426   100    continue
 427          arclen = arclen + gr(k)*arcsum
 428   110    continue
 429       end do
 430       if (arclen .eq. 0.0d0)  goto 170
 431       if (jb .eq. 0)  goto 150
 432 c
 433 c     find number of independent boundaries and check connectivity
 434 c
 435       j = 0
 436       do k = 1, jb
 437          if (kout(k) .ne. 0) then
 438             i = k
 439   120       continue
 440             m = kout(i)
 441             kout(i) = 0
 442             j = j + 1
 443             do ii = 1, jb
 444                if (m .eq. kent(ii)) then
 445                   if (ii .eq. k) then
 446                      ib = ib + 1
 447                      if (j .eq. jb)  goto 150
 448                      goto 130
 449                   end if
 450                   i = ii
 451                   goto 120
 452                end if
 453             end do
 454   130       continue
 455          end if
 456       end do
 457       ib = ib + 1
 458 c
 459 c     attempt to fix connectivity error by moving atom slightly
 460 c
 461       if (moved) then
 462          write (iout,140)  ir
 463   140    format (/,' SURFATOM  --  Connectivity Error at Atom',i6)
 464       else
 465          moved = .true.
 466          xr = xr + rmove
 467          yr = yr + rmove
 468          zr = zr + rmove
 469          goto 10
 470       end if
 471 c
 472 c     compute the exposed surface area for the sphere of interest
 473 c
 474   150 continue
 475       area = ib*pix2 + exang + arclen
 476       area = mod(area,4.0d0*pi) * rrsq
 477 c
 478 c     attempt to fix negative area by moving atom slightly
 479 c
 480       if (area .lt. 0.0d0) then
 481          if (moved) then
 482             write (iout,160)  ir
 483   160       format (/,' SURFATOM  --  Negative Area at Atom',i6)
 484          else
 485             moved = .true.
 486             xr = xr + rmove
 487             yr = yr + rmove
 488             zr = zr + rmove
 489             goto 10
 490          end if
 491       end if
 492   170 continue
 493       return
 494       end