c----------------------------------------------------------------------
      subroutine rot_mol(co,nuatom,angxd,angyd,angzd,center)
c----------------------------------------------------------------------
c
c        Rotation is performed about cartesian coordinate axes,
c     using the geometrical center of the molecule as origin.
c     ANGXD, ANGYD and ANGZD are assumed to be in degrees.
c
      logical center
      logical wireframe,ball_and_stick,need_sort,label_num,
     &                label_cha,need_lim,perspective,round_bonds,
     &                autoscl,showhb,havehb
      dimension co(3,nuatom)
c
      real max_z,min_z
c
      common/frame/scale,xori,yori,max_z,min_z
      common/draw_log/wireframe,ball_and_stick,need_sort,label_num,
     &                label_cha,need_lim,perspective,round_bonds,
     &                autoscl,showhb,havehb
      common/perspective/eye_to_scr,scr_to_top
c
      if (center) then
c
c     Find the center coordinates :
c
           xmax=co(1,1)
           xmin=co(1,1)
           ymax=co(2,1)
           ymin=co(2,1)
           zmax=co(3,1)
           zmin=co(3,1)
           do i=2,nuatom
                xmax=max(xmax,co(1,i))
                xmin=min(xmin,co(1,i))
                ymax=max(ymax,co(2,i))
                ymin=min(ymin,co(2,i))
                zmax=max(zmax,co(3,i))
                zmin=min(zmin,co(3,i))
           enddo
           xorigi=0.5*(xmax-xmin)
           yorigi=0.5*(ymax-ymin)
           zorigi=0.5*(zmax-zmin)
      endif
      call rot_3d(co,nuatom,angxd,angyd,angzd)
c
c     if the rotation involved the X or Y axis the limits on Z are
c     changed.  Get the new values :
c
      if (angx.ne.0.0e0 .or. angy.ne.0.0e0) then
           call get_limits(3,min_z,max_z)
      endif
      scr_to_top = 0.0e0
      eye_to_scr = 2.0e0*(max_z - min_z)
      return
      end
c-----------------------------------------------------------------------
      subroutine rot_3d(cart,number,angxd,angyd,angzd)
c-----------------------------------------------------------------------
c
c        Rotation is performed about cartesian coordinate axes,
c     ANGXD, ANGYD and ANGZD are assumed to be in degrees.
c
      dimension cart(3,number)
c
c     Set up the transformation matrix :
c
      if(angxd .eq. 0.0e0) then
         c1 = 1.0e0
         s1 = 0.0e0
      else
         c1=cosd(angxd)
         s1=sind(angxd)
      endif
      if(angyd .eq. 0.0e0) then
         c2= 1.0e0
         s2= 0.0e0
      else
         c2=cosd(angyd)
         s2=sind(angyd)
      endif
      if(angzd .eq. 0.0e0) then
         c3= 1.0e0
         s3= 0.0e0
      else
         c3=cosd(angzd)
         s3=sind(angzd)
      endif
c
      rot11=c2*c3
      rot21=c2*s3
      rot32=s1*c2
      rot33=c1*c2

      rot12=-s1*s2*c3-c1*s3
      rot13=s1*s3-c1*s2*c3
      rot22=c1*c3-s1*s2*s3
      rot23=-c1*s2*s3-s1*c3

      rot31=s2
c
c     Do the transformation :
c
      do ia=1,number
           coxold=cart(1,ia)
           coyold=cart(2,ia)
           cozold=cart(3,ia)
           cart(1,ia)=rot11*coxold+rot12*coyold+rot13*cozold
           cart(2,ia)=rot21*coxold+rot22*coyold+rot23*cozold
           cart(3,ia)=rot31*coxold+rot32*coyold+rot33*cozold
      enddo
      return
      end
c----------------------------------------------------------------------
      subroutine rank_seq(n)
c----------------------------------------------------------------------
c
c          This is an auxiliary subroutine for sorting.  Given a real
c     array CO of (3 X N) elements, it will create two complementary 
c     pointer lists (ISEQ and IRANK) that will aid in either sorting or 
c     ranking of the input array.  
c            
c          The ISEQ array contains pointers to the positions of the 
c     elements in CO in order of increasing magnitude, i. e. CO(3,SEQ(1))
c     is the smallest element of CO, X(3,SEQ(2)) is the second, etc.
c        
c          The RANK array contains the relative ranking of its 
c     corresponding element of CO, again in order of increasing 
c     magnitude, i. e. CO(3,1) is the RANK(1)'th smallest element of CO,
c     CO(3,i) is the smallest element of CO if RANK(i)=1, etc.
c     
c          Neither RANK nor SEQ need be initialized before calling this
c     subroutine, since they will be overwritten; CO and N will be 
c     returned unchanged.
c     
c          The lists are generated via a modified double-ended search/sort
c     algorithm, in which the logical array USED keeps track of the elements
c     already sequenced/ranked.
c     
c                                                   Julian Tirado-Rives
c                                                    Purdue University
c                                                       March, 1987
c     
c----------------------------------------------------------------------
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
c     logical used,foundmin,foundmax
      logical foundmin,foundmax
c
c     dimension used(natom)
c
      common/sortseq/iseq(natom),irank(natom)
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
c
c     Initializing the arrays :
c
      do j=1,n
c          used(j)  = .false.
           iseq(j)  = 0
           irank(j) = 0
      enddo
c
c     Initializing pointers :
c
c          IFOR  : current forward position in ISEQ
c          IBACK : current backward position in ISEQ
c          IMIN  : position of current minimum in X
c          IMAX  : position of current maximum in X
c
      ifor  = 1
      iback = n
      imin  = ifor
      imax  = iback
c
c     since this is a double-ended search, loop until the backward counter
c     is smaller than the forward.  At this point the fills in both
c     directions have already met in the middle.
c
      do while (iback.ge.ifor)
c
c          search for current minimum and maximum :
c
           foundmin=.false.
           foundmax=.false.
           do j=1,n
c             if (.not.used(j)) then
              if (irank(j).eq.0) then
                 if (j.ne.imin .and. co(3,j).le.co(3,imin)) then
                    imin = j
                    foundmin=.true.
                 endif
                 if (j.ne.imax .and. co(3,j).gt.co(3,imax)) then
                    imax = j
                    foundmax=.true.
                 endif
              endif
           enddo
c          if (.not.foundmin) foundmin = .not.used(imin)
c          if (.not.foundmax) foundmax = .not.used(imax)
           if (.not.foundmin) foundmin = irank(imin).eq.0
           if (.not.foundmax) foundmax = irank(imax).eq.0
           
           if ((.not.foundmin) .or. (.not.foundmax)) then
              do j=1,n
c                if (.not.used(j)) then
                 if (irank(j).eq.0) then
                    if (.not.foundmin) then
                       foundmin=.true.
                       imin=j
                    else if (.not.foundmax) then
                       foundmax=.true.
                       imax=j
                    endif
                 endif
              enddo
           endif
c
c          flag the current minimum and maximum as USED :
c
c          used(imin)=.true.
c          used(imax)=.true.
c
c          fill the ISEQ array with the positions of current minimum and
c          maximum :
c
           iseq(ifor)=imin
           iseq(iback)=imax
c
c          set the ranks of these elements :
c
           irank(imin)=ifor
           irank(imax)=iback
c
c          increment the pointers :
c
           ifor=ifor+1
           iback=iback-1
c
c          swap minimum and maximum to insure finding correct values in
c          the next iteration :
c
           itemp=imin
           imin=imax
           imax=itemp
      enddo
      return
      end
c-----------------------------------------------------------------------
      subroutine calc_slb(atomi,atomj,slbi,slbj)
c-----------------------------------------------------------------------
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      common/draw_log/wireframe,ball_and_stick,need_sort,label_num,
     &                label_cha,need_lim,perspective,round_bonds,
     &                autoscl,showhb,havehb
      common/per_tbl_num/radii(105)
      common/draw_spec/rad_fact,bond_fact,top_z,bottom_z
c
      integer atomi,atomj
      logical wireframe,ball_and_stick,need_sort,label_num,
     &                label_cha,need_lim,perspective,round_bonds,
     &                autoscl,showhb,havehb
c
      dimension slbi(3),slbj(3)
c
      xi = co(1,atomi)
      yi = co(2,atomi)
      zi = co(3,atomi)
      xj = co(1,atomj)
      yj = co(2,atomj)
      zj = co(3,atomj)
      ri  = atrad(atomi)
      rj  = atrad(atomj)
      
      if (perspective) then
           xi = persp(xi,co(3,atomi))
           xj = persp(xj,co(3,atomj))
           yi = persp(yi,co(3,atomi))
           yj = persp(yj,co(3,atomj))
           ri = persp(ri,co(3,atomi))
           rj = persp(rj,co(3,atomj))
      endif
c
      w=radii(1)*bond_fact*bond_fact*rad_fact
      wi=w
      wj=w
      if (perspective) then
           wi=persp(wi,co(3,atomi))*persp(1.0,co(3,atomi))
           wj=persp(wj,co(3,atomj))*persp(1.0,co(3,atomj))
      endif
c      
      rnorm = 1.0e0/((xi-xj)**2+(yi-yj)**2+(zi-zj)**2)
      ridis = sqrt((ri*ri-wi*wi)*rnorm)
      rjdis = sqrt((rj*rj-wj*wj)*rnorm)
      slbi(1)=xi+(xj-xi)*ridis
      slbi(2)=yi+(yj-yi)*ridis
      slbi(3)=zi+(zj-zi)*ridis
      slbj(1)=xj+(xi-xj)*rjdis
      slbj(2)=yj+(yi-yj)*rjdis
      slbj(3)=zj+(zi-zj)*rjdis
c

      return
      end
c----------------------------------------------------------------------
      subroutine calc_bond(bond,xb,yb,numpbond)
c----------------------------------------------------------------------
c
c          This subroutine calculates the verteces of the polygon that
c     represent the bonds.  Briefly, it takes the 3D-single-line bond,
c     calculates an iso-Z line orthogonal to it at each
c     end and selects two equidistant points 0.3*R away from it. This
c     gives the required four points.  If rounded bonds are required
c     the "cap" corresponding to the lower atom is expanded
c
c     points 1 and 4 are the iso-z endpoints in the bottom atom,
c     points 2 and 3  "   "    "       "      "  "  top     "
c
c     the lines 1-2 and 3-4 are the sides of the bond, while
c     the lines 2-3 and 4-?-1 correspond to the "caps" of the bond .
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      integer bond
      logical vertical
      logical wireframe,ball_and_stick,need_sort,label_num,
     &                label_cha,need_lim,perspective,round_bonds,
     &                autoscl,showhb,havehb
c
      common/bonds/numbonds,ibond(nbond),jbond(nbond),nbonds(nbond),
     &             lastbond(natom),numpoint
      common/draw_spec/rad_fact,bond_fact,top_z,bottom_z
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      common/per_tbl_num/radii(105)
      common/draw_log/wireframe,ball_and_stick,need_sort,label_num,
     &                label_cha,need_lim,perspective,round_bonds,
     &                autoscl,showhb,havehb
      common/sortseq/iseq(natom),irank(natom)
c
      dimension xb(20),yb(20),xcap(20),ycap(20)
      dimension slbi(3),slbj(3)
c
      if (irank(ibond(bond)).lt.irank(jbond(bond))) then
           i=ibond(bond)
           j=jbond(bond)
      else
           j=ibond(bond)
           i=jbond(bond)
      endif
      ibnd=2*(bond-1)+1
      jbnd=ibnd+1
c
      call calc_slb(i,j,slbi,slbj)
      x1=slbi(1)
      y1=slbi(2)
      z1=slbi(3)
      x2=slbj(1)
      y2=slbj(2)
      z2=slbj(3)
c
c     building a trapeze :
c
      call line_eq(x1,y1,x2,y2,slope,yint,vertical)
      w=radii(1)*bond_fact*bond_fact*rad_fact
      add1=w
      add2=w
      if (perspective) then
           add1=persp(add1,co(3,i))*persp(1.0,co(3,i))
           add2=persp(add2,co(3,j))*persp(1.0,co(3,j))
      endif
c      
      if (vertical) then
           xb(1)=x1-add1
           yb(1)=y1
           xb(2)=x2-add2
           yb(2)=y2
           xb(3)=x2+add2
           yb(3)=y2
           xb(4)=x1+add1
           yb(4)=y1
      else
           theta = atan(slope)
           costheta = cos(theta)
           sintheta = sin(theta)
c
           xb(1)=x1-add1*sintheta
           yb(1)=y1+add1*costheta
           xb(4)=x1+add1*sintheta
           yb(4)=y1-add1*costheta
           xb(2)=x2-add2*sintheta
           yb(2)=y2+add2*costheta
           xb(3)=x2+add2*sintheta
           yb(3)=y2-add2*costheta
      endif
c
c     if rounded bonds are needed go and calculate them, add the
c     calculated points at the end of the xb and yb arrays.  The
c     bonds where the two atoms have the same Z don't have to be
c     calculated.
c
      if (round_bonds .and. (co(3,j).ne.co(3,i)) ) then
           call calc_bond_cap(slbi,slbj,add1,numpcap,xcap,ycap)
           numpbond=4+numpcap
           do k=1,numpcap
                xb(k+4)=xcap(k)
                yb(k+4)=ycap(k)
           enddo
      else
           numpbond=4
      endif
      return
      end
c----------------------------------------------------------------------
      subroutine calc_bond_cap(bottom,top,width,cap_points,xcap,ycap)
c----------------------------------------------------------------------
c
c     Calculation of the rounded "caps" of bonds :
c
c     input  :  BOTTOM     contains the coordinates of the center of
c                          the bond cap for the atom in the bottom
c               TOP        contains the coordinates of the center of
c                          the bond cap for the atom on top
c               WIDTH      the half width of the bond
c
c     output :  CAP_POINTS the number of points to be added to the
c                          bond-polygon
c               XCAP, YCAP 2-D coordinates of these points.
c
c     Calculation :
c
c           a 3-D unitary cap (values are sines and cosines) centered
c              at the origin in the Y-Z plane is used as a template
c           this template is then multiplied by the wanted width
c           the cap is then rotated in the y axis according to the
c           angle that the central bond line makes with the X-Y plane
c           after that it is rotated in the z axis to match the
c           orthogonal to the slope of the line in the X-Y plane
c           3-D translated to coincide with the center of the bond, and
c           finally, the 2-D values extracted.
c
      integer cap_points
c
      dimension unity(9),unitz(9),xcap(20),ycap(20),cap(3,20)
      dimension bottom(3),top(3)
c
      data unity/-0.95106,-0.80902,-0.58779,-0.30902, 0.0e0,
     &            0.30902, 0.58779, 0.80902, 0.95106/
      data unitz/0.30902, 0.58779, 0.80902,0.95106, 1.0e0,
     &           0.95106, 0.80902, 0.58779, 0.30902/
      data rad_to_deg /57.295779515e+00/ 
c
c     scale the cap according to the width :
c
      do i=1,9
           cap(1,i)=0.0e0
           cap(2,i)=unity(i)*width
           cap(3,i)=unitz(i)*width
      enddo
c
c     do the rotations :
c
      phi=atan((top(3)-bottom(3))/sqrt((top(1)-bottom(1))**2+(top(2)
     &         -bottom(2))**2))*rad_to_deg
      if (top(1).gt.bottom(1)) phi = - phi
      if (top(1).eq.bottom(1)) then
         theta=90.0
         if (top(2).gt.bottom(2)) phi = phi+180.0
      else
           theta=atan((top(2)-bottom(2))/(top(1)-bottom(1)))*
     &           rad_to_deg 
      endif
      call rot_3d(cap,9,0.0e0,-phi,0.0e0,.false.)
      call rot_3d(cap,9,0.0e0,0.0e0,theta,.false.)
c
c     translate
c
      do i=1,9
           xcap(i)=cap(1,i)+bottom(1)
           ycap(i)=cap(2,i)+bottom(2)
      enddo
      cap_points=9
      return
      end
c----------------------------------------------------------------------
      subroutine line_eq(x1,y1,x2,y2,slope,yint,vertical)
c----------------------------------------------------------------------
c 
c          Calculation of the equation of a line (slope & Y-intercept).
c
      logical vertical
c
      if (x1.eq.x2) then
                vertical=.true.
                slope=0.0e0
                yint=0.0e0
           else
                vertical=.false.
                slope=(y1-y2)/(x1-x2)
                yint=y1-(slope*x1)
      endif
      return
      end
c-f--------------------------------------------------------------------
      logical function in_sphere(atom,radius,Ncenters,centers)
c-f--------------------------------------------------------------------
c
c     Will return .TRUE. if ATOM is within a sphere of RADIUS Angstroms
c     of any one of the CENTERS :
c
      integer atom,centers,test_center
      logical show
c
      dimension centers(Ncenters)
c
      show = .false.
      do i=1,Ncenters
           test_center=centers(i)
           if (.not. show) show=(atom.eq.test_center .or. 
     &                           (at_dist(atom,test_center).le.radius))
      enddo
      in_sphere=show
      return
      end
c-f--------------------------------------------------------------------
      real function at_dist(atom1,atom2)
c-f--------------------------------------------------------------------
c
c     Determine the distance between two atoms
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      integer atom1,atom2
c
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
c
      at_dist  = sqrt( (co(1,atom1)-co(1,atom2))**2 + 
     &                 (co(2,atom1)-co(2,atom2))**2 +
     &                 (co(3,atom1)-co(3,atom2))**2  )
      return
      end
c-f--------------------------------------------------------------------
      real function angle(atom1,atom2,atom3)
c-f--------------------------------------------------------------------
c
c          Will determine the angle formed by three atoms using the
c     law of cosines.
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      integer atom1,atom2,atom3
c
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
      data rad_to_deg /57.295779515e+00/
c
      a = at_dist(atom1,atom2)
      b = at_dist(atom2,atom3)
      c = at_dist(atom1,atom3)
c
      costht = (b**2 + a**2 - c**2)/(2 * a * b)
      angle = (acos(costht))*rad_to_deg
      return
      end
c-f--------------------------------------------------------------------
      real function dihedral(atom1,atom2,atom3,atom4)
c-f--------------------------------------------------------------------
c
c     will determine the dihedral angle formed by four atoms using
c     the directional cosines of the vectors normal to the two
c     planes.
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      real labc,lbcd,mabc,mbcd,nabc,nbcd,mgnabc,mgnbcd
      integer atom1,atom2,atom3,atom4
c
      common/const/pi,twopi,degtorad
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
c
      xa = co(1,atom1)-co(1,atom2)
      xc = co(1,atom3)-co(1,atom2)
      xd = co(1,atom4)-co(1,atom2)
      ya = co(2,atom1)-co(2,atom2)
      yc = co(2,atom3)-co(2,atom2)
      yd = co(2,atom4)-co(2,atom2)
      za = co(3,atom1)-co(3,atom2)
      zc = co(3,atom3)-co(3,atom2)
      zd = co(3,atom4)-co(3,atom2)
c
      labc = (ya*zc) - (yc*za)
      mabc = (za*xc) - (zc*xa)
      nabc = (xa*yc) - (xc*ya)
      lbcd = (yd*zc) - (yc*zd)
      mbcd = (zd*xc) - (zc*xd)
      nbcd = (xd*yc) - (xc*yd)
c
      mgnabc = sqrt(labc**2 + mabc**2 + nabc**2)
      mgnbcd = sqrt(lbcd**2 + mbcd**2 + nbcd**2)
c
      cosang = (labc*lbcd+mabc*mbcd+nabc*nbcd)/(mgnabc*mgnbcd)
c
      if (abs(cosang).gt.1.0) cosang = sign (1.0,cosang)
c
c     the variable S is the sign, it can have values of -1, 0, +1 
c
      sg = (xd-xc)*labc + (yd-yc)*mabc + (zd-zc)*nabc
      if (abs(sg).le.0.00001) then 
           s=1.0
      else  
           s= (-1.0)*sign(1.0,sg)
      endif
c
      if (cosang.eq.0.) then
           ang = s * pi/2.0
      else
           ang = s * acos(cosang)
      endif
      dihedral = ang*180.0/pi 
c
      return
      end
c-----------------------------------------------------------------------
      subroutine get_limits(index,minimum,maximum)
c-----------------------------------------------------------------------
c
c          Will return the limits (MINIMUM and MAXIMUM) in the coordinate
c     indicated by INDEX ( Index = 1, 2, and 3 for X, Y, and Z respectively)
c
      parameter (natom = 10000)
      parameter (nbond = 50000)
c
      real minimum,maximum
      integer atom
c
      common/mol_num/numat,num(natom),co(3,natom),atrad(natom)
c
      minimum = co(index,1)
      maximum = minimum
      do atom = 2, numat
c          maximum = max(coordinate,co(index,atom))
c          minimum = min(coordinate,co(index,atom))
           maximum = max(maximum,co(index,atom))
           minimum = min(minimum,co(index,atom))
      enddo
      return
      end
c-f---------------------------------------------------------------------
      real function persp(var,z)
c-f---------------------------------------------------------------------
c
      real max_z,min_z
c
      common/frame/scale,xori,yori,max_z,min_z
      common/perspective/eye_to_scr,scr_to_top
c
      persp = eye_to_scr * scale * var / (eye_to_scr * scale +
     &        scr_to_top * scale + max_z*scale - z*scale)
      return
      end
c-f---------------------------------------------------------------------
      real function persp_x(x,z)
c-f---------------------------------------------------------------------
c
      real max_z,min_z
c
      common/frame/scale,xori,yori,max_z,min_z
      common/perspective/eye_to_scr,scr_to_top
c
      persp_x = persp(x,z)
      return
      end
c-f---------------------------------------------------------------------
      real function persp_y(y,z)
c-f---------------------------------------------------------------------
c
      real max_z,min_z
c
      common/frame/scale,xori,yori,max_z,min_z
      common/perspective/eye_to_scr,scr_to_top
c
      persp_y=persp(y,z)
      return
      end