C
C  James O. Jensen
C  jojensen@crdec.apgea.army.mil
C  Chemical Research, Development and Engineering Center
C  SMCCR-RSP-C
C  Aberdeen Proving Ground
C  Maryland   21010-5423
C
C  Janet L. Jensen
C  jljensen@crdec6.apgea.army.mil
C  Chemical Research, Development and Engineering Center
C  SMCCR-DDT
C  Aberdeen Proving Ground
C  Maryland   21010-5423
C
      subroutine zcon(a,natoms)
C
C  THIS ROUTINE CONVERTS A GENERAL SET OF CARTESIAN COORDINATES
C   INTO Z-MATRIX ORIENTATION THAT IS COMPATIBLE WITH
C   GAUSSIAN90.
C
C  Atom #1 will end up at the origin.
C
C  Atom #2 will end up on the positive Z-axis.
C
C  Atom #3 will end up on the XZ-Plane and have positive 
C   x coordinate. 
C
      implicit double precision (a-h,o-z)
C
      dimension a(3,100)
C
      pid4 = datan(1.0d0)
      pid2 = pid4 + pid4
      pi = pid2 + pid2 
C
C  PUT ATOM #1 AT THE ORIGIN 
C    all other atoms also move
C
      scrx = a(1,1)
      scry = a(2,1)
      scrz = a(3,1)
C
        do 300 i=1,natoms
        a(1,i)= a(1,i) - scrx 
        a(2,i)= a(2,i) - scry 
        a(3,i)= a(3,i) - scrz 
  300   continue
C
C  ROTATE ATOM #2 TO THE YZ-PLANE 
C
      if (a(3,2).ne.0.0d0) then
        theta1 = datan(a(1,2)/a(3,2))
      else
        theta1 = pid2
      endif 
C
      cost1 = dcos(theta1)
      sint1 = dsin(theta1)
C
        do 400 i=2,natoms
        scrx = a(1,i)
        scrz = a(3,i) 
        a(1,i) = cost1*scrx - sint1*scrz
        a(3,i) = sint1*scrx + cost1*scrz  
  400   continue 
C
C  ROTATE ATOM #2 TO THE Z-AXIS 
C 
      if (a(3,2).ne.0.0d0) then
        theta2 = datan(a(2,2)/a(3,2))
      else
        theta2 = pid2
      endif 
C
      cost2 = dcos(theta2)
      sint2 = dsin(theta2)
C
        do 500 i=2,natoms
        scry = a(2,i)
        scrz = a(3,i) 
        a(2,i) = cost2*scry - sint2*scrz
        a(3,i) = sint2*scry + cost2*scrz
  500   continue 
C
C  ROTATE ATOM #3 TO THE XZ-PLANE
C
      if (a(1,3).ne.0.0d0) then
        theta3 = datan(a(2,3)/a(1,3))
      else
        theta3 = pid2
      endif 
C
C  Make sure atom #3 is pointing in the positive x direction
C   relative to atom #1 and atom #2.
C
      if (a(1,3).lt.0.0d0) then 
        sinn = -1.0d0
      else 
        sinn = 1.0d0
      endif  
C
      cost3 = dcos(theta3)
      sint3 = dsin(theta3)
C
        do 600 i=3,natoms
        scry = a(2,i)
        scrx = a(1,i) 
        a(2,i) = (cost3*scry - sint3*scrx) * sinn
        a(1,i) = (sint3*scry + cost3*scrx) * sinn 
  600   continue 
C
C
C  Make sure that atom #2 ends up pointing in the positive Z direction
C    relative to atom #1.  Otherwise rotate 180 degrees.  
C
      if (a(3,2).lt.0.0d0) then 
        sinn = -1.0d0
      else 
        sinn = 1.0d0
      endif  
C
        do 700 i=1,natoms
        a(2,i) = a(2,i) * sinn 
        a(3,i) = a(3,i) * sinn  
  700   continue 
C
C  Make sure atom #3 is pointing in the positive x direction
C   relative to atom #1 and atom #2.  Otherwise rotate 180 degrees.
C
      if (a(1,3).lt.0.0d0) then 
        sinn = -1.0d0
      else 
        sinn = 1.0d0
      endif  
C
        do 800 i=1,natoms
        a(2,i) = a(2,i) * sinn
        a(1,i) = a(1,i) * sinn 
  800   continue 
C
      return
      end