avrsph.f

最经典的分子对结软件

***********************************************************************
* this routine of the main recluster averages spheres whose centers   *
* and radii are close to one anothers.  'close' for radii and centers *
* are user defined variable parameters.                               *
***********************************************************************
      subroutine avrsph
      include 'cluster.h'
      integer i,j,l,k,sphcnt,kk
ccc       loop variables
      integer numavr, totavr
ccc       numavr: number of spheres averaged, per-residue.
ccc       totavr: total number of spheres averaged
      real dist
ccc         dist: the distance btwn a sphere center and and average
ccc               center.
      real arad, crad, disrad
ccc         arad: temp variable for radii.
ccc         crad: radius of average sphere.
ccc       disrad: difference btwn two radii.
      real temp(3)
ccc         temp: temp variable for average sphere.
      real center(3)
ccc       center: average sphere for residue.
      logical logavr
ccc       logavr: whether any spheres have been averaged.

ccc   run thru spheres until no more averaging occurs.
   25 continue
ccc       initialize variables.
*         do 10 k=1,3
*              temp(k)=0.0
*              center(k)=0.0
*  10     continue
          do 20 j=1,redsph 
               do 18 i=1,5
                    do 15 k=1,7 
                         clsxyz(j,i,k)=0.0
   15               continue
   18          continue
   20     continue
          arad=0.0
          i=0
          sphcnt=1
          totavr=0
          logavr=.false.
ccc       for all spheres.
   50     continue
               do 60 k=1,3
                    temp(k)=0.0
                    center(k)=0.0
   60          continue
               j=0
               numavr=0
ccc            for all spheres.
   75          continue
                    arad=numavr*crad+trads(sphcnt+j)
                    arad=arad/(numavr+1)
                    disrad=arad-trads(sphcnt+j)
                    dist=0.0
                    do 100 k=1,3
                         temp(k)=numavr*center(k)+tspcor(sphcnt+j,k)
                         temp(k)=temp(k)/(numavr+1)
                         dist=dist+(temp(k)-tspcor(sphcnt+j,k))**2
  100               continue
                    dist=sqrt(dist)
                
ccc                 if the distance and radii criteria are fulfilled,
ccc                 .and. neither of the spheres are flagged.
                    if(dist.le.nearnb.and.disrad.le.nearad
     &              .and. .not.tflagd(sphcnt+j) .and. .not.
     &              tflagd(sphcnt))then
                         numavr=numavr+1

ccc                      keep track of closest spheres to avrsph.
                         if(numavr.eq.2)then
                             logavr=.true.
                             do 110 k=1,3
                                 clsxyz(sphcnt,2,k)=tspcor(sphcnt+j,k)
                                 clsxyz(sphcnt,1,k)=tspcor(sphcnt,k)
  110                        continue
                             clsxyz(sphcnt,2,4)=trads(sphcnt+j)
                             clsxyz(sphcnt,1,4)=trads(sphcnt)
                             clsxyz(sphcnt,1,5)=dist
                             clsxyz(sphcnt,2,5)=dist
                             clsxyz(sphcnt,1,6)=tiatom(sphcnt)
                             clsxyz(sphcnt,2,6)=tiatom(sphcnt+j)
                             clsxyz(sphcnt,1,7)=tjatom(sphcnt)
                             clsxyz(sphcnt,2,7)=tjatom(sphcnt+j)
                             clsxyz(sphcnt,1,8)=tcolor(sphcnt)
                             clsxyz(sphcnt,2,8)=tcolor(sphcnt+j)
                         elseif(numavr.le.5)then
                             do 115 k=1,3
                                  clsxyz(sphcnt,numavr,k)=
     &                               tspcor(sphcnt+j,k)
  115                        continue
                             clsxyz(sphcnt,numavr,4)=trads(sphcnt+j)
                             clsxyz(sphcnt,numavr,5)=dist
                             clsxyz(sphcnt,numavr,6)=tiatom(sphcnt+j)
                             clsxyz(sphcnt,numavr,7)=tjatom(sphcnt+j)
                             clsxyz(sphcnt,numavr,8)=tcolor(sphcnt+j)
                         else
                             call reords(sphcnt)
                             k=0
  125                        continue
                                 k=k+1
                                 if(dist.lt.clsxyz(sphcnt,k,5))then
ccc                                 push the distance list down.
                                    do 140 l=5,k+1,-1
                                         do 130 kk=1,5
                                              clsxyz(sphcnt,l,kk)=
     &                                        clsxyz(sphcnt,l-1,kk)
  130                                    continue
  140                               continue
                                    do 145 kk=1,3
                                         clsxyz(sphcnt,k,kk)=
     &                                     tspcor(sphcnt+j,kk)
  145                               continue
                                    clsxyz(sphcnt,k,4)=trads(sphcnt+j)
                                    clsxyz(sphcnt,k,5)=dist
ccc                                 loop breakout condition.
                                    k=6
                                 endif
                             if(k.lt.5)goto 125
                         endif

                         if(numavr.gt.1)then
                              do 200 k=sphcnt+j,nmax-1
                                   tiatom(k)=tiatom(k+1)
                                   trads(k)=trads(k+1)
                                   tjatom(k)=tjatom(k+1)
                                   tcolor(k)=tcolor(k+1)
                                   tflagd(k)=tflagd(k+1)
                                   do 150 l=1,3
                                        tspcor(k,l)=tspcor(k+1,l)
  150                              continue
  200                         continue
                              j=j-1
                              nmax=nmax-1
                         endif
                         do 240 k=1,3
                            center(k)=temp(k)
  240                    continue
                         crad=arad
                    endif

                    j=j+1
              if(sphcnt+j.lt.nmax)goto 75

  350         format('for residue', i5,1x,i5,1x,'spheres averaged')
  375         format('center coordinates',4f10.5)
              if(numavr.gt.1)then
                   totavr=totavr+numavr

                   if(.not.tflagd(sphcnt))then
                        do 390 k=1,3
                             tspcor(sphcnt,k)=clsxyz(sphcnt,1,k)
  390                   continue
                        trads(sphcnt)=clsxyz(sphcnt,1,4)
                        tiatom(sphcnt)=clsxyz(sphcnt,1,6)
                        tjatom(sphcnt)=clsxyz(sphcnt,1,7)
                        tcolor(sphcnt)=clsxyz(sphcnt,1,8)
                        tflagd(sphcnt)=.false.
                   endif
              endif
              sphcnt = sphcnt + 1
          if(sphcnt.lt.nmax)goto 50
      if(logavr)goto 25
      write(6,*)'number of spheres is -',nmax,' after averaging'

      return
      end