getclu.f

最经典的分子对结软件

**********************************************************************
*---------------------------getclu-----------------------------------*
* this sub. sorts the spheres by radii, puting the spheres < then    *
* mxrad1 in spcors and those smaller in the various x arrays.        *
* bks, 08/89                                                         *
**********************************************************************
      subroutine getclu(mxrad1,mxrad2)
      include 'cluster.h'
      
      real mxrad1,mxrad2
ccc       mxrad1: user defined max-radius.
ccc       mxrad2: secondary max radius (spheres can be included, but 
ccc               cannot lead to further spheres being added.
      integer i, j, k, l
ccc         loop variables

ccc   initialise variables
      k=1
      l=1
      do 120 i = 1, nmax
           iatoms(k)=tiatom(i)
           label(k)=i
           do 53 j=1,3
                spcors(k,j)=tspcor(i,j)
   53      continue
           rads(k)=trads(i)
           jatoms(k)=tjatom(i)
           scolor(k)=tcolor(i)
           flaged(k)=tflagd(i)
ccc        only spheres smaller than the user defined radius.
           if(rads(k).le.mxrad1)then
                k=k+1
                if (k .gt. redsph) then
                  write (6,74) redsph
   74             format('cluster(getclu): reduced sphere limit ',i6)
                  write (6,84) 
   84             format('exceeded, program stops')
                  stop
                endif
           elseif(rads(k).le.mxrad2)then
                if (l .gt. redsph) then
                  write (6,74) redsph
                  write (6,84) 
                  stop
                endif
                xiatom(l)=iatoms(k)      
                xjatom(l)=jatoms(k) 
                xcolor(l)=scolor(k) 
                do 95 j=1,3         
                    xspcor(l,j)=spcors(k,j)
   95           continue                         
                xtrads(l)=rads(k)
                xflagd(l)=tflagd(k)
                l=l+1                                     
           endif      
  120 continue                                 
      nsph=k-1
      xsph=l-1
      write(6,175),nsph,xsph,mxrad1,mxrad2
  175 format('nsph,xsph,mxrad1,mxrad2 are:',2i8,2f10.3)

      return
      end