PROGRAM hepta1
      IMPLICIT NONE
C Version 1.0
C August 1996
C
C by Feliu Maseras (Univ. Montpellier II)  
C maseras@lsd.univ-montp2.fr              
C 
C Fits the best possible coordination polyhedron 
C (pentagonal bipyramid, capped octahedron, capped trigonal prism) 
C to a given seven-coordinate geometry. 
C
C Input: Cartesian coordinates of the 8 atomic centers implied in 
C 7 coordination, starting with the central atom. 
C Each line starts with an arbitrary integer number, followed by 
C 3 real numbers corresponding to the cartesian coordinates. Input 
C is free format. 
C
C Example of input file:
C 1           3.10125   1.74646   3.22288
C 2           3.10125    .88272    .51223
C 3           3.10125   3.44358   5.44623
C 4           3.10125   -.23984   3.10382
C 5           1.49604   1.12256   4.29579
C 6           4.50674   3.18476   2.41439
C 7          4.70646   1.12256   4.29579
C 8          1.69576   3.18476   2.41439
C
C Example of compilation of the program in an IBM RISC workstation:
C prompt> xlf -o hepta1.exe hepta1.f
C
C Example of execution of the program in a unix system:
C (output to screen):
C prompt> hepta1.exe < example.hee 
C
C Example of execution of the program in a unix system:
C (output to file):
C prompt> hepta1.exe < example.hee > example.hes
C
      REAL*8 coord(3,8),vedir(3,7),chose(7,7),sum,pi
      INTEGER i1,i2,i3,iprint,iaux
C Writes header
      WRITE(6,5)
    5 FORMAT(2X,'*********************************************',/,
     1       2X,'*  HEPTA  version 1.0 (August 1996)         *',/,
     2       2X,'*                                           *',/,
     3       2X,'*  Users are kindly requested to refer the  *',/,
     4       2X,'*  following paper:                         *',/,
     5       2X,'*  F. Maseras, O. Eisenstein; New J. Chem., *',/,
     6       2X,'*  21, 961-967 (1997).                      *',/,
     7       2X,'*********************************************',/,
     8       2X,/,
     9       72('*'))
C Initialize pi
      pi= 4.*ATAN(1.)
C Initializes iprint
      iprint= 5
C
C Reads coordinates
C
      DO i1= 1,8
        READ(5,*) iaux,(coord(i2,i1),i2=1,3)
        END DO 
C Debug printing
      IF(iprint.GE.8) THEN
        WRITE(6,110)
  110 FORMAT(2X,'Input cartesian coordinates')
        DO i1= 1,8
          WRITE(6,120) i1,(coord(i2,i1),i2=1,3)
  120 FORMAT(2X,I2,3F10.6)
          END DO 
        END IF
C 
C Fills the matrix of directional vectors vedir, and normalizes it
C
      DO i1= 2,8
        i3= i1-1
        DO i2= 1,3
          vedir(i2,i3)= coord(i2,i1) - coord(i2,1)
          END DO 
        END DO 
C Debug printing
      IF(iprint.GE.9) THEN
        WRITE(6,310)
  310 FORMAT(2X,'Non-normalized matrix of directional vectors')
        DO i1= 1,8
          WRITE(6,120) i1,(vedir(i2,i1),i2=1,3)
          END DO 
        END IF
C Normalizes
      DO i1= 1,7
        sum= 0.
        DO i2= 1,3
          sum= sum + vedir(i2,i1)*vedir(i2,i1)
          END DO 
        sum= DSQRT(sum)
C
        IF(sum.LT.1.D-05) THEN
          i3= i1+1
          WRITE(6,10) 1,i3,sum 
   10 FORMAT(2X,'Distance between atoms',I3,' and',I3,' is only',D10.4)
          STOP
          END IF
C
        DO i2= 1,3
          vedir(i2,i1)= vedir(i2,i1)/sum
          END DO 
        END DO 
C Debug printing
      IF(iprint.GE.8) THEN
        WRITE(6,210)
  210 FORMAT(2X,'Normalized matrix of directional vectors')
        DO i1= 1,8
          WRITE(6,120) i1,(vedir(i2,i1),i2=1,3)
          END DO 
        END IF
C
C Fills up the chose square matrix, containing the angles 
C between each pair of ligands. 
C
      CALL cang(vedir,chose)
C Debug printing
      IF(iprint.GE.8) THEN
        WRITE(6,410)
  410 FORMAT(2X,'Ligand-ligand angles in input coordinates (radians)')
        WRITE(6,420) (i1,i1=1,7)
  420 FORMAT(2X,4X,7I8)
        DO i1= 1,7
          WRITE(6,430) i1,(chose(i2,i1),i2=1,7)
  430 FORMAT(2X,i4,7F8.4)
          END DO 
        END IF
      IF(iprint.GE.2) THEN
        WRITE(6,510)
  510 FORMAT(2X,'Ligand-ligand angles in input coordinates (degrees)')
        WRITE(6,420) (i1,i1=1,7)
        DO i1= 1,7
          WRITE(6,530) i1,(((180.*chose(i2,i1))/(pi)),i2=1,7)
  530 FORMAT(2X,i4,7F8.2)
          END DO 
        END IF
C
      CALL pbp(chose,iprint)
      CALL co(chose,iprint)
      CALL ctp(chose,iprint,vedir)
      CALL cofi(chose,iprint)
      CALL ctpfi(chose,iprint,vedir)
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE pbp(chose,iprint)
      IMPLICIT NONE
C Finds the best possible pentagonal bipyramid
C 
C pbpser(i)   Error associated to the choice of atom i as the top 
C             axial
C pbpord(j,i) Optimal ordering of the ligands when i is axial.
C             In increasing order of angle.
C             pbpord(1,i) is always i.
C             pbpord(7,i) is the other axial.
C             pbpord(2,i) to pbpord(6,i) are the equatorials
C             pbpord(2,i) is closest to 1 (in principle arbitrary)
C             pbpord(3,i) is the closest between 3-6 to 2.
C             pbpord(4,i) is the closest between 4-6 to 3
C             pbpord(5,i) is the closest between 5-6 to 4.
C chase(k,l)  Auxiliary matrix to be compared with chose. It contains 
C             the ideal values of ligand-ligand angles for a given 
C             choice of axial ligand. 
C chese(k,l)  Auxiliary matrix used in the definition of chase. 
C             ligands are placed in their ideal positions in order to 
C             compute their ideal angles.
C
      REAL*8 chose(7,7),pbpser(7)
      REAL*8 bpval,pi,chase(7,7),chese(3,7),pi25,difaux
      INTEGER i1,i2,i3,iprint,i3min,pbpord(7,7)
      INTEGER bpbp
C
      WRITE(6,1234)
C
      WRITE(6,123) 
  123 FORMAT(2x,'D5H PENTAGONAL BIPYRAMID')
C
C Initialize pi
      pi= 4*ATAN(1.)
      pi25= (2*pi)/5.
C
C Large loop pbp1
C
      DO i1=1,7
C i1 is the candidate for the top axial position
C
C Sorts the ligands in the pbpord(?,i1) vector
C
C Initialize pbpord(?,i1)
        DO i2=1,7
          pbpord(i2,i1)= i2
          END DO 
C
        CALL luck(pbpord,chose,i1,i1,1,7)
        CALL luck(pbpord,chose,i1,2,3,6)
        CALL luck(pbpord,chose,i1,3,4,6)
        CALL luck(pbpord,chose,i1,4,5,6)
C
C Fills the matrix with the optimal values for positions
C (in the precise order of ligands)
C
        CALL pbpide(pbpord,i1,chese)
C Debug printing
      IF(iprint.GE.8) THEN
        WRITE(6,210) i1
  210 FORMAT(2X,'Ideal coordinates for PBP with axial',I3)
        DO i2= 1,7
          WRITE(6,120) i2,(chese(i3,i2),i3=1,3)
  120 FORMAT(2X,I2,3F10.6)
          END DO 
        END IF
C 
C Fills the matrix chase with the optimal values for the angles 
C
      CALL cang(chese,chase)
C Debug printing
      IF(iprint.GE.10) THEN
        WRITE(6,410) i1
  410 FORMAT(2X,'Ideal angles (radians) for PBP with axial',i3) 
        WRITE(6,420) (i2,i2=1,7)
  420 FORMAT(2X,4X,7I8)
        DO i2= 1,7
          WRITE(6,430) i2,(chase(i2,i3),i3=1,7)
  430 FORMAT(2X,i4,7F8.4)
          END DO 
        END IF
C Evaluates deviation between real and optimal
        pbpser(i1)= 0.
        DO i2=1,7
          i3min= i2+1
          DO i3=i3min,7
            difaux= chose(i3,i2) - chase(i3,i2)
            pbpser(i1)= pbpser(i1) + difaux*difaux
            END DO 
          END DO 
C Scales pbpser to standard deviation in degrees
        pbpser(i1)= (180./pi)*(sqrt(pbpser(i1)/21.))
        END DO
C
C Decides the best pbp
C
      bpbp= 1
      bpval= pbpser(1)
      DO i1=1,7
        IF(pbpser(i1).LT.bpval) THEN
          bpval= pbpser(i1)
          bpbp= i1
          END IF
        END DO 
C 
C Output
C
      IF(iprint.GE.0) THEN
        WRITE(6,10) bpbp,bpval
   10 FORMAT(2X,'Best pentagonal bipyramid with ligand',I2,
     1 ' in axial position',/,2X,'Deviation associated:',F10.4)
        END IF
C
      IF(iprint.GE.2) THEN
        WRITE(6,110) 
  110 FORMAT(2X,'Axial ligands:')
        WRITE(6,125) pbpord(1,bpbp),pbpord(7,bpbp)
  125 FORMAT(2X,2I5)
        WRITE(6,130) 
  130 FORMAT(2X,'Equatorial ligands:')
        WRITE(6,140) (pbpord(i2,bpbp),i2=2,6)
  140 FORMAT(2X,5I5)
        END IF
C
      IF(iprint.GE.2) THEN
        CALL pbpide(pbpord,bpbp,chese)
        CALL cang(chese,chase)
        WRITE(6,510)
  510 FORMAT(2X,'Ideal ligand-ligand angles (degrees)')
        WRITE(6,520) (i1,i1=1,7)
  520 FORMAT(2X,4X,7I8)
        DO i1= 1,7
          WRITE(6,530) i1,(((180.*chase(i2,i1))/(pi)),i2=1,7)
  530 FORMAT(2X,i4,7F8.2)
          END DO
        END IF
C
      IF(iprint.GE.5) THEN
        WRITE(6,20)
   20 FORMAT(2X,'Deviation associated to each possible PBP')
        DO i1= 1,7
          WRITE(6,30) i1,pbpser(i1)
   30 FORMAT(2X,'Ligand in axial position:',I2,5X,
     1 'Deviation associated:',F8.4)
          END DO
        END IF
      RETURN
C
 1234 FORMAT(72('*'))
C
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE co(chose,iprint)
      IMPLICIT NONE
C axial=capping
C
C Finds the best possible capped octahedron
C 
C cocser(i)   Error associated to the choice of atom i as the top 
C             axial
C cocord(j,i) Optimal ordering of the ligands when i is axial.
C             cocord(1,i) is always i.
C             cocord(2,i) to cocord(4,i) are the capped face.
C             cocord(5,i) to cocord(7,i) are the non capped face.
C             cocord(2,i) is closest to 1 (in principle arbitrary)
C             cocord(5,i) is the closest between 5-7 to 2.
C             cocord(3,i) is the closest between 3-4 to 5.
C             cocord(6,i) is the closest between 6-7 to 3
C             Ideal value of capping-center-capped is average 
C               between 1-2, 1-3, 1-4
C             Ideal value of capping-center-noncapped is average 
C               between 1-5, 1-6, 1-7
C chase(k,l)  Auxiliary matrix to be compared with chose. It contains 
C             the ideal values of ligand-ligand angles for a given 
C             choice of axial ligand. 
C chese(k,l)  Auxiliary matrix used in the definition of chase. 
C             ligands are placed in their ideal positions in order to 
C             compute their ideal angles.
C
      REAL*8 chose(7,7),cocser(7)
      REAL*8 cocval,pi,chase(7,7),chese(3,7),pi25,difaux
      REAL*8 axca(7),axnca(7),sum
      INTEGER i1,i2,i3,iprint,i3min,cocord(7,7),iaux
      INTEGER bco
C
      WRITE(6,1234)
C
      WRITE(6,123) 
  123 FORMAT(2x,'C3V CAPPED OCTAHEDRON')
C
C Initialize pi
      pi= 4*ATAN(1.)
      pi25= (2*pi)/5.
C
C Large loop coc1
C
      DO i1=1,7
C i1 is the candidate for the top axial position
C
C Sorts the ligands in the cocord(?,i1) vector
C
C Initialize cocord(?,i1)
        DO i2=1,7
          cocord(i2,i1)= i2
          END DO 
C
C Orders the ligands in cocord(?,i1)
C
        CALL luck(cocord,chose,i1,i1,1,7)
        CALL luck(cocord,chose,i1,2,5,7)
        CALL luck(cocord,chose,i1,5,3,4)
        CALL luck(cocord,chose,i1,3,6,7)
C Computes the ideal angles with axial axca(i1), axnca(i1)
        sum= 0.
        DO i3= 2,4
          iaux= cocord(i3,i1)
          sum= sum + chose(iaux,i1)
          END DO
        axca(i1)= sum/3.
C
        sum= 0.
        DO i3= 5,7
          iaux= cocord(i3,i1)
          sum= sum + chose(iaux,i1)
          END DO
        axnca(i1)= sum/3.
C
C Fills the matrix with the optimal values for positions
C (in the precise order of ligands)
C
        CALL cocide(cocord,i1,chese,axca,axnca)
C Debug printing
      IF(iprint.GE.8) THEN
        WRITE(6,210) i1
  210 FORMAT(2X,'Ideal coordinates for CO with capping',I3)
        WRITE(6,220) (360./(2.*pi))*axca(i1),(360./(2.*pi))*axnca(i1)
  220 FORMAT(2X,'Ideal values for capping-capped and capping-noncapped:',
     1 /,2F10.3)
        DO i2= 1,7
          WRITE(6,120) i2,(chese(i3,i2),i3=1,3)
  120 FORMAT(2X,I2,3F10.6)
          END DO 
        END IF
C 
C Fills the matrix chase with the optimal values for the angles 
C
      CALL cang(chese,chase)
C Debug printing
      IF(iprint.GE.10) THEN
        WRITE(6,410) i1
  410 FORMAT(2X,'Ideal angles (radians) for CO with capping',i3) 
        WRITE(6,420) (i2,i2=1,7)
  420 FORMAT(2X,4X,7I8)
        DO i2= 1,7
          WRITE(6,430) i2,(chase(i2,i3),i3=1,7)
  430 FORMAT(2X,i4,7F8.4)
          END DO 
        END IF
C Evaluates deviation between real and optimal
        cocser(i1)= 0.
        DO i2=1,7
          i3min= i2+1
          DO i3=i3min,7
            difaux= chose(i3,i2) - chase(i3,i2)
            cocser(i1)= cocser(i1) + difaux*difaux
            END DO 
          END DO 
C Scales cocser to standard deviation in degrees
        cocser(i1)= (180./pi)*(sqrt(cocser(i1)/21.))
        END DO
C
C Decides the best pbp
C
      bco= 1
      cocval= cocser(1)
      DO i1=1,7
        IF(cocser(i1).LT.cocval) THEN
          cocval= cocser(i1)
          bco= i1
          END IF
        END DO 
C 
C Output
C
      IF(iprint.GE.0) THEN
        WRITE(6,10) bco,cocval
   10 FORMAT(2X,'Best capped octahedron with ligand',I2,
     1 ' in capping position',/,2X,'Deviation associated:',F10.4)
        END IF
C
      IF(iprint.GE.2) THEN
        WRITE(6,110) 
  110 FORMAT(2X,'Capping ligand:')
        WRITE(6,125) cocord(1,bco)
  125 FORMAT(2X,2I5)
        WRITE(6,130) 
  130 FORMAT(2X,'Cappped face ligands:')
        WRITE(6,140) (cocord(i2,bco),i2=2,4)
  140 FORMAT(2X,3I5)
        WRITE(6,150) 
  150 FORMAT(2X,'Non-cappped face ligands:')
        WRITE(6,140) (cocord(i2,bco),i2=5,7)
        WRITE(6,160) (360./(2.*pi))*axca(bco),(360./(2.*pi))*axnca(bco)
  160 FORMAT(2X,'Values for capping-capped and capping-noncap angles:',
     1 /,2F10.3)
        END IF
C
      IF(iprint.GE.2) THEN
        CALL cocide(cocord,bco,chese,axca,axnca)
        CALL cang(chese,chase)
        WRITE(6,510)
  510 FORMAT(2X,'Ideal ligand-ligand angles (degrees)')
        WRITE(6,520) (i1,i1=1,7)
  520 FORMAT(2X,4X,7I8)
        DO i1= 1,7
          WRITE(6,530) i1,(((180.*chase(i2,i1))/(pi)),i2=1,7)
  530 FORMAT(2X,i4,7F8.2)
          END DO
        END IF
C
      IF(iprint.GE.5) THEN
        WRITE(6,20)
   20 FORMAT(2X,'Deviation associated to each possible CO')
        DO i1= 1,7
          WRITE(6,30) i1,cocser(i1)
   30 FORMAT(2X,'Ligand in capping position:',I2,5X,
     1 'Deviation associated:',F8.4)
          END DO
        END IF
      RETURN
C
 1234 FORMAT(72('*'))
C
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE ctp(chose,iprint,vedir)
      IMPLICIT NONE
C (Kind of a mess concerning terminology:
C    axial=capping
C    equatorial=quadrilateral face
C    basal=edge)
C
C Finds the best possible capped trigonal prism
C 
C ctpser(i)   Error associated to the choice of atom i as the top 
C             axial
C ctpord(j,i) Optimal ordering of the ligands when i is axial.
C             ctpord(1,i) is always i.
C             ctpord(2,i) to ctpord(5,i) are the equatorial ligands.
C             ctpord(6,i) to ctpord(7,i) are the basal ligands.
C             ctpord(2,i) is closest to 1 (in principle arbitrary)
C             ctpord(6,i) is the closest between 6-7 to 2.
C             ctpord(3,i) is the closest between 3-5 to 6.
C             ctpord(4,i) is the closest between 4-5 to 3
C             Ideal value axeq(i1) of axial-center-equatorial is average 
C               between 1-2, 1-3, 1-4,1-5
C             Ideal value axbas(i1) of axial-center-basal is average 
C               between 1-6, 1-7
C             Ideal value dieq(i1) of basal-center-axial-basal 
C              dihedral angle is more complicated to compute 
C              (see subroutine cdieq).
C chase(k,l)  Auxiliary matrix to be compared with chose. It contains 
C             the ideal values of ligand-ligand angles for a given 
C             choice of axial ligand. 
C chese(k,l)  Auxiliary matrix used in the definition of chase. 
C             ligands are placed in their ideal positions in order to 
C             compute their ideal angles.
C vedir is passed only because it must be fed to cdieq
C
      REAL*8 chose(7,7),ctpser(7)
      REAL*8 ctpval,pi,chase(7,7),chese(3,7),pi25,difaux
      REAL*8 vedir(3,7)
      REAL*8 axeq(7),axbas(7),dieq(7),sum
      INTEGER i1,i2,i3,iprint,i3min,ctpord(7,7),iaux
      INTEGER bctp
C
      WRITE(6,1234)
C
      WRITE(6,123) 
  123 FORMAT(2x,'C2V CAPPED TRIGONAL PRISM')
C
C Initialize pi
      pi= 4*ATAN(1.)
      pi25= (2*pi)/5.
C
C Large loop ctp1
C
      DO i1=1,7
C i1 is the candidate for the top axial position
C
C Sorts the ligands in the ctpord(?,i1) vector
C
C Initialize ctpord(?,i1)
        DO i2=1,7
          ctpord(i2,i1)= i2
          END DO 
C
C Orders the ligands in ctpord(?,i1)
C
        CALL luck(ctpord,chose,i1,i1,1,7)
        CALL luck(ctpord,chose,i1,2,6,7)
        CALL luck(ctpord,chose,i1,6,3,5)
        CALL luck(ctpord,chose,i1,3,4,5)
C Computes the ideal angles with axial axeq(i1), axbas(i1)
        sum= 0.
        DO i3= 2,5
          iaux= ctpord(i3,i1)
          sum= sum + chose(iaux,i1)
          END DO
        axeq(i1)= sum/4.
C
        sum= 0.
        DO i3= 6,7
          iaux= ctpord(i3,i1)
          sum= sum + chose(iaux,i1)
          END DO
        axbas(i1)= sum/2.
C Computes the ideal dihedral angle dieq(i1)
        CALL cdieq(chose,ctpord,i1,iprint,vedir,dieq)
C
C Fills the matrix with the optimal values for positions
C (in the precise order of ligands)
C
        CALL ctpide(ctpord,i1,chese,axeq,axbas,dieq)
C Debug printing
      IF(iprint.GE.8) THEN
        WRITE(6,210) i1
  210 FORMAT(2X,'Ideal coordinates for CTP with capping',I3)
        WRITE(6,220) (360./(2.*pi))*axeq(i1),(360./(2.*pi))*axbas(i1)
  220 FORMAT(2X,'Ideal values for capping-quad.face and capping-edge:',
     1 /,2F10.3)
        WRITE(6,230) (360./(2.*pi))*dieq(i1)
  230 FORMAT(2X,'Ideal value for quadrilateral face dihedral:',F10.3)
        DO i2= 1,7
          WRITE(6,120) i2,(chese(i3,i2),i3=1,3)
  120 FORMAT(2X,I2,3F10.6)
          END DO 
        END IF
C 
C Fills the matrix chase with the optimal values for the angles 
C
      CALL cang(chese,chase)
C Debug printing
      IF(iprint.GE.10) THEN
        WRITE(6,410) i1
  410 FORMAT(2X,'Ideal angles (radians) for CTP with capping',i3) 
        WRITE(6,420) (i2,i2=1,7)
  420 FORMAT(2X,4X,7I8)
        DO i2= 1,7
          WRITE(6,430) i2,(chase(i2,i3),i3=1,7)
  430 FORMAT(2X,i4,7F8.4)
          END DO 
        END IF
C Evaluates deviation between real and optimal
        ctpser(i1)= 0.
        DO i2=1,7
          i3min= i2+1
          DO i3=i3min,7
            difaux= chose(i3,i2) - chase(i3,i2)
            ctpser(i1)= ctpser(i1) + difaux*difaux
            END DO 
          END DO 
C Scales ctpser to standard deviation in degrees
        ctpser(i1)= (180./pi)*(sqrt(ctpser(i1)/21.))
        END DO
C
C Decides the best ctp
C
      bctp= 1
      ctpval= ctpser(1)
      DO i1=1,7
        IF(ctpser(i1).LT.ctpval) THEN
          ctpval= ctpser(i1)
          bctp= i1
          END IF
        END DO 
C 
C Output
C
      IF(iprint.GE.0) THEN
        WRITE(6,10) bctp,ctpval
   10 FORMAT(2X,'Best capped trigonal prism with ligand',I2,
     1 ' in capping position',/,2X,'Deviation associated:',F10.4)
        END IF
C
      IF(iprint.GE.2) THEN
        WRITE(6,110) 
  110 FORMAT(2X,'Capping ligand:')
        WRITE(6,125) ctpord(1,bctp)
  125 FORMAT(2X,2I5)
        WRITE(6,130) 
  130 FORMAT(2X,'Quad.face ligands:')
        WRITE(6,140) (ctpord(i2,bctp),i2=2,5)
  140 FORMAT(2X,4I5)
        WRITE(6,150) 
  150 FORMAT(2X,'Edge ligands:')
        WRITE(6,140) (ctpord(i2,bctp),i2=6,7)
      WRITE(6,160) (360./(2.*pi))*axeq(bctp),(360./(2.*pi))*axbas(bctp)
  160 FORMAT(2X,'Values for capping-quad.face and capping-edge angles:',
     1 /,2F10.3)
        WRITE(6,170) (360./(2.*pi))*dieq(bctp)
  170 FORMAT(2X,'Value for quadrilateral face dihedral angle:',F10.3) 
        END IF
C
      IF(iprint.GE.2) THEN
        CALL ctpide(ctpord,bctp,chese,axeq,axbas,dieq)
        CALL cang(chese,chase)
        WRITE(6,510)
  510 FORMAT(2X,'Ideal ligand-ligand angles (degrees)')
        WRITE(6,520) (i1,i1=1,7)
  520 FORMAT(2X,4X,7I8)
        DO i1= 1,7
          WRITE(6,530) i1,(((180.*chase(i2,i1))/(pi)),i2=1,7)
  530 FORMAT(2X,i4,7F8.2)
          END DO
        END IF
C
      IF(iprint.GE.5) THEN
        WRITE(6,20)
   20 FORMAT(2X,'Deviation associated to each possible CTP')
        DO i1= 1,7
          WRITE(6,30) i1,ctpser(i1)
   30 FORMAT(2X,'Ligand in capping position:',I2,5X,
     1 'Deviation associated:',F8.4)
          END DO
        END IF
      RETURN
C
 1234 FORMAT(72('*'))
C
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE cofi(chose,iprint)
      IMPLICIT NONE
C axial=capping
C
C Finds the best possible capped octahedron
C WITH FIXED STANDARD VALUES FOR THE FREE ANGLES
C 
C cocser(i)   Error associated to the choice of atom i as the top 
C             axial
C cocord(j,i) Optimal ordering of the ligands when i is axial.
C             cocord(1,i) is always i.
C             cocord(2,i) to cocord(4,i) are the capped face.
C             cocord(5,i) to cocord(7,i) are the non capped face.
C             cocord(2,i) is closest to 1 (in principle arbitrary)
C             cocord(5,i) is the closest between 5-7 to 2.
C             cocord(3,i) is the closest between 3-4 to 5.
C             cocord(6,i) is the closest between 6-7 to 3
C             Ideal value of capping-center-capped is average 
C               between 1-2, 1-3, 1-4
C             Ideal value of capping-center-noncapped is average 
C               between 1-5, 1-6, 1-7
C chase(k,l)  Auxiliary matrix to be compared with chose. It contains 
C             the ideal values of ligand-ligand angles for a given 
C             choice of axial ligand. 
C chese(k,l)  Auxiliary matrix used in the definition of chase. 
C             ligands are placed in their ideal positions in order to 
C             compute their ideal angles.
C
      REAL*8 chose(7,7),cocser(7)
      REAL*8 cocval,pi,chase(7,7),chese(3,7),pi25,difaux
      REAL*8 axca(7),axnca(7),sum,axcaf,axncaf
      INTEGER i1,i2,i3,iprint,i3min,cocord(7,7),iaux
      INTEGER bco
C
      WRITE(6,1234)
C
      WRITE(6,123) 
  123 FORMAT(2x,'C3V CAPPED OCTAHEDRON (STANDARD ANGLES)')
C
C Initialize pi
      pi= 4*ATAN(1.)
      pi25= (2*pi)/5.
C Initialize axcaf, axncaf
C (currently, Bartell & Thompson with n=1)
      axcaf= ((2.*pi)/360.)*73.3
      axncaf= ((2.*pi)/360.)*128.8
C
C Large loop coc1
C
      DO i1=1,7
C i1 is the candidate for the top axial position
C
C Sorts the ligands in the cocord(?,i1) vector
C
C Initialize cocord(?,i1)
        DO i2=1,7
          cocord(i2,i1)= i2
          END DO 
C
C Orders the ligands in cocord(?,i1)
C
        CALL luck(cocord,chose,i1,i1,1,7)
        CALL luck(cocord,chose,i1,2,5,7)
        CALL luck(cocord,chose,i1,5,3,4)
        CALL luck(cocord,chose,i1,3,6,7)
C Computes the ideal angles with axial axca(i1), axnca(i1)
C       sum= 0.
C       DO i3= 2,4
C         iaux= cocord(i3,i1)
C         sum= sum + chose(iaux,i1)
C         END DO
        axca(i1)= axcaf
C
C       sum= 0.
C       DO i3= 5,7
C         iaux= cocord(i3,i1)
C         sum= sum + chose(iaux,i1)
C         END DO
        axnca(i1)= axncaf
C
C Fills the matrix with the optimal values for positions
C (in the precise order of ligands) 
C
        CALL cocide(cocord,i1,chese,axca,axnca)
C Debug printing
      IF(iprint.GE.8) THEN
        WRITE(6,210) i1
  210 FORMAT(2X,'Ideal coordinates for CO with capping',I3)
        WRITE(6,220) (360./(2.*pi))*axca(i1),(360./(2.*pi))*axnca(i1)
  220 FORMAT(2X,'Fixed values for capping-capped and capping-noncapped:',
     1 /,2F10.3)
        DO i2= 1,7
          WRITE(6,120) i2,(chese(i3,i2),i3=1,3)
  120 FORMAT(2X,I2,3F10.6)
          END DO 
        END IF
C 
C Fills the matrix chase with the optimal values for the angles 
C
      CALL cang(chese,chase)
C Debug printing
      IF(iprint.GE.10) THEN
        WRITE(6,410) i1
  410 FORMAT(2X,'Fixed angles (radians) for CO with capping',i3) 
        WRITE(6,420) (i2,i2=1,7)
  420 FORMAT(2X,4X,7I8)
        DO i2= 1,7
          WRITE(6,430) i2,(chase(i2,i3),i3=1,7)
  430 FORMAT(2X,i4,7F8.4)
          END DO 
        END IF
C Evaluates deviation between real and optimal
        cocser(i1)= 0.
        DO i2=1,7
          i3min= i2+1
          DO i3=i3min,7
            difaux= chose(i3,i2) - chase(i3,i2)
            cocser(i1)= cocser(i1) + difaux*difaux
            END DO 
          END DO 
C Scales cocser to standard deviation in degrees
        cocser(i1)= (180./pi)*(sqrt(cocser(i1)/21.))
        END DO
C
C Decides the best pbp
C
      bco= 1
      cocval= cocser(1)
      DO i1=1,7
        IF(cocser(i1).LT.cocval) THEN
          cocval= cocser(i1)
          bco= i1
          END IF
        END DO 
C 
C Output
C
      IF(iprint.GE.0) THEN
        WRITE(6,10) bco,cocval
   10 FORMAT(2X,'Best capped octahedron with ligand',I2,
     1 ' in capping position',/,2X,'Deviation associated:',F10.4)
        END IF
C
      IF(iprint.GE.2) THEN
        WRITE(6,110) 
  110 FORMAT(2X,'Capping ligand:')
        WRITE(6,125) cocord(1,bco)
  125 FORMAT(2X,2I5)
        WRITE(6,130) 
  130 FORMAT(2X,'Cappped face ligands:')
        WRITE(6,140) (cocord(i2,bco),i2=2,4)
  140 FORMAT(2X,3I5)
        WRITE(6,150) 
  150 FORMAT(2X,'Non-cappped face ligands:')
        WRITE(6,140) (cocord(i2,bco),i2=5,7)
        WRITE(6,160) (360./(2.*pi))*axca(bco),(360./(2.*pi))*axnca(bco)
  160 FORMAT(2X,'Values for capping-capped and capping-noncap angles:',
     1 /,2F10.3)
        END IF
C
      IF(iprint.GE.2) THEN
        CALL cocide(cocord,bco,chese,axca,axnca)
        CALL cang(chese,chase)
        WRITE(6,510)
  510 FORMAT(2X,'Ideal ligand-ligand angles (degrees)')
        WRITE(6,520) (i1,i1=1,7)
  520 FORMAT(2X,4X,7I8)
        DO i1= 1,7
          WRITE(6,530) i1,(((180.*chase(i2,i1))/(pi)),i2=1,7)
  530 FORMAT(2X,i4,7F8.2)
          END DO
        END IF
C
      IF(iprint.GE.5) THEN
        WRITE(6,20)
   20 FORMAT(2X,'Deviation associated to each possible CO')
        DO i1= 1,7
          WRITE(6,30) i1,cocser(i1)
   30 FORMAT(2X,'Ligand in capping position:',I2,5X,
     1 'Deviation associated:',F8.4)
          END DO
        END IF
      RETURN
C
 1234 FORMAT(72('*'))
C
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE ctpfi(chose,iprint,vedir)
      IMPLICIT NONE
C (Kind of a mess concerning terminology:
C    axial=capping
C    equatorial=quadrilateral face
C    basal=edge)
C
C Finds the best possible capped trigonal prism 
C WITH FIXED STANDARD VALUES FOR FREE ANGLES
C 
C ctpser(i)   Error associated to the choice of atom i as the top 
C             axial
C ctpord(j,i) Optimal ordering of the ligands when i is axial.
C             ctpord(1,i) is always i.
C             ctpord(2,i) to ctpord(5,i) are the equatorial ligands.
C             ctpord(6,i) to ctpord(7,i) are the basal ligands.
C             ctpord(2,i) is closest to 1 (in principle arbitrary)
C             ctpord(6,i) is the closest between 6-7 to 2.
C             ctpord(3,i) is the closest between 3-5 to 6.
C             ctpord(4,i) is the closest between 4-5 to 3
C             Ideal value axeq(i1) of axial-center-equatorial is average 
C               between 1-2, 1-3, 1-4,1-5
C             Ideal value axbas(i1) of axial-center-basal is average 
C               between 1-6, 1-7
C             Ideal value dieq(i1) of basal-center-axial-basal 
C              dihedral angle is more complicated to compute 
C              (see subroutine cdieq).
C chase(k,l)  Auxiliary matrix to be compared with chose. It contains 
C             the ideal values of ligand-ligand angles for a given 
C             choice of axial ligand. 
C chese(k,l)  Auxiliary matrix used in the definition of chase. 
C             ligands are placed in their ideal positions in order to 
C             compute their ideal angles.
C vedir is passed only because it must be fed to cdieq
C
      REAL*8 chose(7,7),ctpser(7)
      REAL*8 ctpval,pi,chase(7,7),chese(3,7),pi25,difaux
      REAL*8 vedir(3,7)
      REAL*8 axeq(7),axbas(7),dieq(7),sum,axeqf,axbasf,dieqf
      INTEGER i1,i2,i3,iprint,i3min,ctpord(7,7),iaux
      INTEGER bctp
C
      WRITE(6,1234)
C
      WRITE(6,123) 
  123 FORMAT(2x,'C2V CAPPED TRIGONAL PRISM (STANDARD ANGLES)')
C
C Initialize pi
      pi= 4*ATAN(1.)
      pi25= (2*pi)/5.
C Initialize axeqf,axbasf,dieqf
C (currently, Thompson&Bartell, n=1)
      axeqf= ((2.*pi)/360.)*80.8
      axbasf= ((2.*pi)/360.)*144.2
      dieqf= ((2.*pi)/360.)*49.0
C
C Large loop ctp1
C
      DO i1=1,7
C i1 is the candidate for the top axial position
C
C Sorts the ligands in the ctpord(?,i1) vector
C
C Initialize ctpord(?,i1)
        DO i2=1,7
          ctpord(i2,i1)= i2
          END DO 
C
C Orders the ligands in ctpord(?,i1)
C
        CALL luck(ctpord,chose,i1,i1,1,7)
        CALL luck(ctpord,chose,i1,2,6,7)
        CALL luck(ctpord,chose,i1,6,3,5)
        CALL luck(ctpord,chose,i1,3,4,5)
C Computes the ideal angles with axial axeq(i1), axbas(i1)
C       sum= 0.
C       DO i3= 2,5
C         iaux= ctpord(i3,i1)
C         sum= sum + chose(iaux,i1)
C         END DO
        axeq(i1)= axeqf
C
C       sum= 0.
C       DO i3= 6,7
C         iaux= ctpord(i3,i1)
C         sum= sum + chose(iaux,i1)
C         END DO
        axbas(i1)= axbasf
C Computes the ideal dihedral angle dieq(i1)
C       CALL cdieq(chose,ctpord,i1,iprint,vedir,dieq)
        dieq(i1)= dieqf
C
C Fills the matrix with the optimal values for positions
C (in the precise order of ligands)
C
        CALL ctpide(ctpord,i1,chese,axeq,axbas,dieq)
C Debug printing
      IF(iprint.GE.8) THEN
        WRITE(6,210) i1
  210 FORMAT(2X,'Ideal coordinates for CTP with capping',I3)
        WRITE(6,220) (360./(2.*pi))*axeq(i1),(360./(2.*pi))*axbas(i1)
  220 FORMAT(2X,'Fixed values for capping-quad.face and capping-edge:',
     1 /,2F10.3)
        WRITE(6,230) (360./(2.*pi))*dieq(i1)
  230 FORMAT(2X,'Fixed value for quadrilateral face dihedral:',F10.3)
        DO i2= 1,7
          WRITE(6,120) i2,(chese(i3,i2),i3=1,3)
  120 FORMAT(2X,I2,3F10.6)
          END DO 
        END IF
C 
C Fills the matrix chase with the optimal values for the angles 
C
      CALL cang(chese,chase)
C Debug printing
      IF(iprint.GE.10) THEN
        WRITE(6,410) i1
  410 FORMAT(2X,'Ideal angles (radians) for CTP with capping',i3) 
        WRITE(6,420) (i2,i2=1,7)
  420 FORMAT(2X,4X,7I8)
        DO i2= 1,7
          WRITE(6,430) i2,(chase(i2,i3),i3=1,7)
  430 FORMAT(2X,i4,7F8.4)
          END DO 
        END IF
C Evaluates deviation between real and optimal
        ctpser(i1)= 0.
        DO i2=1,7
          i3min= i2+1
          DO i3=i3min,7
            difaux= chose(i3,i2) - chase(i3,i2)
            ctpser(i1)= ctpser(i1) + difaux*difaux
            END DO 
          END DO 
C Scales ctpser to standard deviation in degrees
        ctpser(i1)= (180./pi)*(sqrt(ctpser(i1)/21.))
        END DO
C
C Decides the best ctp
C
      bctp= 1
      ctpval= ctpser(1)
      DO i1=1,7
        IF(ctpser(i1).LT.ctpval) THEN
          ctpval= ctpser(i1)
          bctp= i1
          END IF
        END DO 
C 
C Output
C
      IF(iprint.GE.0) THEN
        WRITE(6,10) bctp,ctpval
   10 FORMAT(2X,'Best capped trigonal prism with ligand',I2,
     1 ' in capping position',/,2X,'Deviation associated:',F10.4)
        END IF
C
      IF(iprint.GE.2) THEN
        WRITE(6,110) 
  110 FORMAT(2X,'Capping ligand:')
        WRITE(6,125) ctpord(1,bctp)
  125 FORMAT(2X,2I5)
        WRITE(6,130) 
  130 FORMAT(2X,'Quad.face ligands:')
        WRITE(6,140) (ctpord(i2,bctp),i2=2,5)
  140 FORMAT(2X,4I5)
        WRITE(6,150) 
  150 FORMAT(2X,'Edge ligands:')
        WRITE(6,140) (ctpord(i2,bctp),i2=6,7)
      WRITE(6,160) (360./(2.*pi))*axeq(bctp),(360./(2.*pi))*axbas(bctp)
  160 FORMAT(2X,'Values for capping-quad.face and capping-edge angles:',
     1 /,2F10.3)
        WRITE(6,170) (360./(2.*pi))*dieq(bctp)
  170 FORMAT(2X,'Value for quadrilateral face dihedral angle:',F10.3) 
        END IF
C
      IF(iprint.GE.2) THEN
        CALL ctpide(ctpord,bctp,chese,axeq,axbas,dieq)
        CALL cang(chese,chase)
        WRITE(6,510)
  510 FORMAT(2X,'Ideal ligand-ligand angles (degrees)')
        WRITE(6,520) (i1,i1=1,7)
  520 FORMAT(2X,4X,7I8)
        DO i1= 1,7
          WRITE(6,530) i1,(((180.*chase(i2,i1))/(pi)),i2=1,7)
  530 FORMAT(2X,i4,7F8.2)
          END DO
        END IF
C
      IF(iprint.GE.5) THEN
        WRITE(6,20)
   20 FORMAT(2X,'Deviation associated to each possible CTP')
        DO i1= 1,7
          WRITE(6,30) i1,ctpser(i1)
   30 FORMAT(2X,'Ligand in capping position:',I2,5X,
     1 'Deviation associated:',F8.4)
          END DO
        END IF
      RETURN
C
 1234 FORMAT(72('*'))
C
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE luck(ord,chose,icol,iref,ista,ifin)
      IMPLICIT NONE
C
C Sorts in column icol of matrix ord the ligands from ista to ifin
C according to increasing angle with iref.
C Angles are taken from chose.
C
      REAL*8 chose(7,7)
      INTEGER i2,i3,j2,j3,i3min,ord(7,7),icol,iref,ista,ifin
      INTEGER aux,iaux
C Bubble 
      iaux= ord(iref,icol)
      DO i2=ista,ifin
        i3min= i2+1
        DO i3=i3min,ifin
          j2= ord(i2,icol)
          j3= ord(i3,icol)
          IF(chose(j2,iaux).GT.chose(j3,iaux)) THEN
            aux= ord(i2,icol)
            ord(i2,icol)= ord(i3,icol)
            ord(i3,icol)= aux
            END IF
          END DO 
        END DO 
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE pbpide(pbpord,icol,chese)
      IMPLICIT NONE
C
C Fills the chese matrix with the ideal PBP coordinates.
C (in the precise order of ligands indicated in the icol column 
C of pbpord)
C
      REAL*8 pi,chese(3,7),pi25
      INTEGER pbpord(7,7),iaux,icol
C Initialize pi
      pi= 4*ATAN(1.)
      pi25= (2*pi)/5.
C Fills the matrix with the optimal values for positions
C Axial 1
      iaux= pbpord(1,icol)
      chese(1,iaux)= 0.0
      chese(2,iaux)= 0.0
      chese(3,iaux)= 1.0
C Axial 2
      iaux= pbpord(7,icol)
      chese(1,iaux)= 0.0
      chese(2,iaux)= 0.0
      chese(3,iaux)= -1.0
C Equatorial 1
      iaux= pbpord(2,icol)
      chese(1,iaux)= 1.0
      chese(2,iaux)= 0.0
      chese(3,iaux)= 0.0
C Equatorial 2
      iaux= pbpord(3,icol)
      chese(1,iaux)= COS(pi25)
      chese(2,iaux)= SIN(pi25)
      chese(3,iaux)= 0.0
C Equatorial 3
      iaux= pbpord(4,icol)
      chese(1,iaux)= COS(2*pi25)
      chese(2,iaux)= SIN(2*pi25)
      chese(3,iaux)= 0.0
C Equatorial 4
      iaux= pbpord(5,icol)
      chese(1,iaux)= COS(3*pi25)
      chese(2,iaux)= SIN(3*pi25)
      chese(3,iaux)= 0.0
C Equatorial 5
      iaux= pbpord(6,icol)
      chese(1,iaux)= COS(4*pi25)
      chese(2,iaux)= SIN(4*pi25)
      chese(3,iaux)= 0.0
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE cocide(cocord,icol,chese,axca,axnca)
      IMPLICIT NONE
C
C Fills the chese matrix with the ideal PBP coordinates.
C (in the precise order of ligands indicated in the icol column 
C of cocord)
C
      REAL*8 pi,chese(3,7),pi23,pi13,axca(7),axnca(7)
      INTEGER cocord(7,7),iaux,icol
C Initialize pi
      pi= 4*ATAN(1.)
      pi23= (2*pi)/3.
      pi13= (2*pi)/6.
C Fills the matrix with the optimal values for positions
C Axial 1
      iaux= cocord(1,icol)
      chese(1,iaux)= 0.0
      chese(2,iaux)= 0.0
      chese(3,iaux)= 1.0
C Capped face 1 (theta=axca(icol), phi=0.)
      iaux= cocord(2,icol)
      chese(1,iaux)= SIN(axca(icol))
      chese(2,iaux)= 0.0
      chese(3,iaux)= COS(axca(icol))
C Capped face 2 (theta=axca(icol), phi=(2*pi)/3)
      iaux= cocord(3,icol)
      chese(1,iaux)= SIN(axca(icol))*COS(pi23)
      chese(2,iaux)= SIN(axca(icol))*SIN(pi23)
      chese(3,iaux)= COS(axca(icol))
C Capped face 3 (theta=axca(icol), phi=(4*pi)/3)
      iaux= cocord(4,icol)
      chese(1,iaux)= SIN(axca(icol))*COS(2.*pi23)
      chese(2,iaux)= SIN(axca(icol))*SIN(2.*pi23)
      chese(3,iaux)= COS(axca(icol))
C Non-capped face 1 (theta=axnca(icol), phi=(2*pi)/6.)
      iaux= cocord(5,icol)
      chese(1,iaux)= SIN(axnca(icol))*COS(pi13)
      chese(2,iaux)= SIN(axnca(icol))*SIN(pi13)
      chese(3,iaux)= COS(axnca(icol))
C Non-capped face 2 (theta=axnca(icol), phi=(pi))
      iaux= cocord(6,icol)
      chese(1,iaux)= SIN(axnca(icol))*COS(pi)
      chese(2,iaux)= SIN(axnca(icol))*SIN(pi)
      chese(3,iaux)= COS(axnca(icol))
C Non-capped face 3 (theta=axnca(icol), phi=(10*pi)/6)
      iaux= cocord(7,icol)
      chese(1,iaux)= SIN(axnca(icol))*COS(5.*pi13)
      chese(2,iaux)= SIN(axnca(icol))*SIN(5.*pi13)
      chese(3,iaux)= COS(axnca(icol))
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE ctpide(ctpord,icol,chese,axeq,axbas,dieq)
      IMPLICIT NONE
C
C Fills the chese matrix with the ideal CTP coordinates.
C (in the precise order of ligands indicated in the icol column 
C of ctpord)
C
      REAL*8 pi,chese(3,7),axeq(7),axbas(7),dieq(7)
      INTEGER ctpord(7,7),iaux,icol
C Initialize pi
      pi= 4*ATAN(1.)
C Fills the matrix with the optimal values for positions
C Axial 1
      iaux= ctpord(1,icol)
      chese(1,iaux)= 0.0
      chese(2,iaux)= 0.0
      chese(3,iaux)= 1.0
C Equatorial 1 (theta=axeq(icol), phi=dieq(icol))
      iaux= ctpord(2,icol)
      chese(1,iaux)= SIN(axeq(icol))*COS(dieq(icol))
      chese(2,iaux)= SIN(axeq(icol))*SIN(dieq(icol))
      chese(3,iaux)= COS(axeq(icol))
C Equatorial 2 (theta=axeq(icol), phi=-dieq(icol))
      iaux= ctpord(3,icol)
      chese(1,iaux)= SIN(axeq(icol))*COS(-dieq(icol))
      chese(2,iaux)= SIN(axeq(icol))*SIN(-dieq(icol))
      chese(3,iaux)= COS(axeq(icol))
C Equatorial 3 (theta=axeq(icol), phi=pi+dieq(icol))
      iaux= ctpord(4,icol)
      chese(1,iaux)= SIN(axeq(icol))*COS(pi+dieq(icol))
      chese(2,iaux)= SIN(axeq(icol))*SIN(pi+dieq(icol))
      chese(3,iaux)= COS(axeq(icol))
C Equatorial 4 (theta=axeq(icol), phi=(pi-dieq(icol))
      iaux= ctpord(5,icol)
      chese(1,iaux)= SIN(axeq(icol))*COS(pi-dieq(icol))
      chese(2,iaux)= SIN(axeq(icol))*SIN(pi-dieq(icol))
      chese(3,iaux)= COS(axeq(icol))
C Basal 1 (theta=axbas(icol), phi=0.)
      iaux= ctpord(6,icol)
      chese(1,iaux)= SIN(axbas(icol))
      chese(2,iaux)= 0.
      chese(3,iaux)= COS(axbas(icol))
C Basal 1 (theta=axbas(icol), phi=pi)
      iaux= ctpord(7,icol)
      chese(1,iaux)= -SIN(axbas(icol))
      chese(2,iaux)= 0.
      chese(3,iaux)= COS(axbas(icol))
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE cang(coord,angles)
      IMPLICIT NONE
C 
C Fills the angles square matrix containing the ligand-metal-ligand 
C angles from the normalized direction vectors contained in the 
C coord matrix 
C First it computes the cosinus as the scalar product, then the angle.
C The matrix is triangular, so only the side below the diagonal needs 
C to be computed
C
      REAL*8 angles(7,7),coord(3,7),sum
      INTEGER i2,j1,j2,j2min
C
      DO j1= 1,7
        j2min= j1+1
        DO j2= j2min,7
          sum= 0.
          DO i2= 1,3
            sum= sum + coord(i2,j1)*coord(i2,j2)
            END DO 
          angles(j2,j1)= ACOS(sum)
C Elements above the diagonal are symmetric
          angles(j1,j2)= angles(j2,j1)
          END DO 
C Elements in diagonal are 0
C (cosinus=1., ligand forms an angle of zero with itself)
        angles(j1,j1)= 0.
        END DO 
      RETURN 
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE cdieq(chose,ctpord,icol,iprint,vedir,dieq)
      IMPLICIT NONE
C
C Finds the optimal value for the equatorial-center-axial-equatorial 
C dihedral angle in a CTP
C The algorithm consists in projecting the axial component out of 
C each of the direction vectors of the four equatorial ligands. 
C Afterwards, they remain in a plane, and dihedral angles can be 
C computed directly  as if they were bond angles through scalar 
C products.
C
      REAL*8 chose(7,7),chise(7,7),chuse(3,7)
      REAL*8 dieq(7),sum
      REAL*8 vedir(3,7)
      INTEGER i1,i2,i3,i4,iprint,ctpord(7,7),j1,icol
C
C Fills auxiliary matrix chuse
C
      DO i1=1,7
        DO i2=1,3
          chuse(i2,i1)= vedir(i2,i1)
          END DO
        END DO
C
C Projects out axial component
C
      DO i1=2,5
        j1= ctpord(i1,icol)
        DO i2=1,3
          chuse(i2,j1)= chuse(i2,j1) - 
     1     COS(chose(j1,icol))*chuse(i2,icol)
          END DO 
C Debug printing
        IF(iprint.GE.10) THEN
          WRITE(6,110)
  110     FORMAT(2X,'Basal direction vector after projecting out ',
     1 'the axial component',/,2X,'(not normalized)')
          WRITE(6,120) j1,(chuse(i2,j1),i2=1,3)
  120     FORMAT(2X,i3,3F8.4)
          END IF
C Normalizes the resulting vector
        sum= 0.
        DO i2=1,3
          sum= sum + chuse(i2,j1)*chuse(i2,j1)
          END DO
        sum= DSQRT(sum)
C
        DO i2= 1,3
          chuse(i2,j1)= chuse(i2,j1)/sum
          END DO 
        END DO 
C Debug printing
        IF(iprint.GE.10) THEN
          WRITE(6,210)
  210     FORMAT(2X,'Basal direction vectors after projecting out ',
     1 'the axial component',/,2X,'(normalized)')
          DO i4=2,5
            i3= ctpord(i4,icol)
            WRITE(6,120) i3,(chuse(i2,i3),i2=1,3)
            END DO
          END IF
C
C Computes angles
C
      CALL cang(chuse,chise)
C
C Picks up the results and evaluates dieq
C
      i1= ctpord(2,icol)
      i2= ctpord(3,icol)
      i3= ctpord(4,icol)
      i4= ctpord(5,icol)
      dieq(icol)= (abs(chise(i1,i2))+abs(chise(i3,i4))) / 4.
      RETURN
      END