some gfortran warning fixes
[unres.git] / 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