The Hartree-Fock program from Szabo-Ostlund

C********************************************************************
C
C MINIMAL BASIS STO-3G CALCULATION ON HEH+
C
C THIS IS A LITTLE DUMMY MAIN PROGRAM WHICH CALLS HFCALC
C
C APPENDIX B: TWO-ELECTRON SELF-CONSISTENT-FIELD PROGRAM
C OF MODERN QUANTUM CHEMISTRY by
C Attila Szabo and Neil S. Ostlund
C Ed. 2nd (1989) Dover Publications INC.
C
C Labourly Typed by Michael Zitolo (Feb., 2005)
C Edited and Compiled by Michael Zitolo and Xihua Chen
C
C Cleaned up and debugged again by Andrew Long (2012) 
C                   and Daniele (kalium) Dondi (2013)
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      IOP=1
      N=3
      R=1.4632D0
      ZETA1=2.0925D0
      ZETA2=1.24D0
      ZA=2.0D0
      ZB=1.0D0
      CALL HFCALC(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
      END

C*********************************************************************
      SUBROUTINE HFCALC(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C
C DOES A HARTREE-FOCK CALCULATION FOR A TWO-ELECTRON DIATOMIC
C USING THE 1S MINIMAL STO-NG BASIS SET
C MINIMAL BASIS SET HAS BASIS FUNCTIONS 1 AND 2 ON NUCLEI A AND B
C
C IOP=0 NO PRINTING WHATSOEVER (TO OPTIMIZE EXPONENTS, SAY)
C IOP=1 PRINT ONLY CONVERGED RESULTS
C IOP=2 PRINT EVERY ITERATION
C N STO-NG CALCULATION (N=1,2 OR 3)
C R BONDLENGTH (AU)
C ZETA1 SLATER ORBITAL EXPONENT (FUNCTION 1)
C ZETA2 SLATER ORBITAL EXPONENT (FUNCTION 2)
C ZA ATOMIC NUMBER (ATOM A)
C ZB ATOMIC NUMBER (ATOM B)
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      IF (IOP.EQ.0) GO TO 20
      PRINT 10,N,ZA,ZB
   10 FORMAT(' ',2X,'STO-',I1,'G FOR ATOMIC NUMBERS ',F5.2,' AND ',
     $ F5.2//)
   20 CONTINUE
C CALCULATE ALL THE ONE AND TWO ELECTRON INTEGRALS
      CALL INTGRL(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C BE INEFFICIENT AND PUT ALL INTEGRALS IN PRETTY ARRAYS
      CALL COLECT(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C PERFORM THE SCF CALCULATION
      CALL SCF(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
      RETURN
      END

C*********************************************************************
      SUBROUTINE INTGRL(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C
C CALCULATES ALL THE BASIC INTEGRALS NEEDED FOR SCF CALCULATION
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/INT/S12,T11,T12,T22,V11A,V12A,V22A,V11B,V12B,V22B,
     $ V1111,V2111,V2121,V2211,V2221,V2222
      DIMENSION COEF(3,3),EXPON(3,3),D1(3),A1(3),D2(3),A2(3)
      DATA PI/3.1415926535898D0/
C THESE ARE THE CONTRACTION COEFFICIENTS AND EXPONENTS FOR
C A NORMALIZED SLATER ORBITAL WITH EXPONENT 1.0 IN TERMS OF
C NORMALIZED 1S PRIMITIVE GAUSSIANS
      DATA COEF,EXPON/1.0D0,2*0.0D0,0.678914D0,0.430129D0,0.0D0,
     $ 0.444635D0,0.535328D0,0.154329D0,0.270950D0,2*0.0D0,0.151623D0,
     $ 0.851819D0,0.0D0,0.109818D0,0.405771D0,2.22766D0/
      R2=R*R
C SCALE THE EXPONENTS (A) OF PRIMITIVE GAUSSIANS
C INCLUDE NORMALIZATION IN CONTRACTION COEFFICIENTS (D)
      DO 10 I=1,N
      A1(I)=EXPON(I,N)*(ZETA1**2)
      D1(I)=COEF(I,N)*((2.0D0*A1(I)/PI)**0.75D0)
      A2(I)=EXPON(I,N)*(ZETA2**2)
      D2(I)=COEF(I,N)*((2.0D0*A2(I)/PI)**0.75D0)
   10 CONTINUE
C D AND A ARE NOW THE CONTRACTION COEFFICIENTS AND EXPONENTS
C IN TERMS OF UNNORMALIZED PRIMITIVE GAUSSIANS
      S12=0.0D0
      T11=0.0D0
      T12=0.0D0
      T22=0.0D0
      V11A=0.0D0
      V12A=0.0D0
      V22A=0.0D0
      V11B=0.0D0
      V12B=0.0D0
      V22B=0.0D0
      V1111=0.0D0
      V2111=0.0D0
      V2121=0.0D0
      V2211=0.0D0
      V2221=0.0D0
      V2222=0.0D0
C CALCULATE ONE-ELECTRON INTEGRALS
C CENTER A IS FIRST ATOM, CETER B IS SECOND ATOM
C ORIGIN IS ON CENTER A
C V12A = OFF-DIAGONAL NUCLEAR ATTRACTION TO CENTER A, ETC.
      DO 20 I=1,N
      DO 20 J=1,N
C RAP2 = SQUARED DISTANCE BETWEEN CENTER A AND CENTER P, ETC.
      RAP=A2(J)*R/(A1(I)+A2(J))
      RAP2=RAP**2
      RBP2=(R-RAP)**2
      S12=S12+S(A1(I),A2(J),R2)*D1(I)*D2(J)
      T11=T11+T(A1(I),A1(J),0.0D0)*D1(I)*D1(J)
      T12=T12+T(A1(I),A2(J),R2)*D1(I)*D2(J)
      T22=T22+T(A2(I),A2(J),0.0D0)*D2(I)*D2(J)
      V11A=V11A+V(A1(I),A1(J),0.0D0,0.0D0,ZA)*D1(I)*D1(J)
      V12A=V12A+V(A1(I),A2(J),R2,RAP2,ZA)*D1(I)*D2(J)
      V22A=V22A+V(A2(I),A2(J),0.0D0,R2,ZA)*D2(I)*D2(J)
      V11B=V11B+V(A1(I),A1(J),0.0D0,R2,ZB)*D1(I)*D1(J)
      V12B=V12B+V(A1(I),A2(J),R2,RBP2,ZB)*D1(I)*D2(J)
      V22B=V22B+V(A2(I),A2(J),0.0D0,0.0D0,ZB)*D2(I)*D2(J)
   20 CONTINUE
C CALCULATE TWO-ELECTRON INTEGRALS
      DO 30 I=1,N
      DO 30 J=1,N
      DO 30 K=1,N
      DO 30 L=1,N
      RAP=A2(I)*R/(A2(I)+A1(J))
      RBP=R-RAP
      RAQ=A2(K)*R/(A2(K)+A1(L))
      RBQ=R-RAQ
      RPQ=RAP-RAQ
      RAP2=RAP*RAP
      RBP2=RBP*RBP
      RAQ2=RAQ*RAQ
      RBQ2=RBQ*RBQ
      RPQ2=RPQ*RPQ
      V1111=V1111+TWOE(A1(I),A1(J),A1(K),A1(L),0.0D0,0.0D0,0.0D0)
     $ *D1(I)*D1(J)*D1(K)*D1(L)
      V2111=V2111+TWOE(A2(I),A1(J),A1(K),A1(L),R2,0.0D0,RAP2)
     $ *D2(I)*D1(J)*D1(K)*D1(L)
      V2121=V2121+TWOE(A2(I),A1(J),A2(K),A1(L),R2,R2,RPQ2)
     $ *D2(I)*D1(J)*D2(K)*D1(L)
      V2211=V2211+TWOE(A2(I),A2(J),A1(K),A1(L),0.0D0,0.0D0,R2)
     $ *D2(I)*D2(J)*D1(K)*D1(L)
      V2221=V2221+TWOE(A2(I),A2(J),A2(K),A1(L),0.0D0,R2,RBQ2)
     $ *D2(I)*D2(J)*D2(K)*D1(L)
      V2222=V2222+TWOE(A2(I),A2(J),A2(K),A2(L),0.0D0,0.0D0,0.0D0)
     $ *D2(I)*D2(J)*D2(K)*D2(L)
   30 CONTINUE
      IF (IOP.EQ.0) GO TO 90
      PRINT 40
   40 FORMAT(3X,'R',10X,'ZETA1',6X,'ZETA2',6X,'S12',8X,'T11'/)
      PRINT 50, R,ZETA1,ZETA2,S12,T11
   50 FORMAT(5F11.6//)
      PRINT 60
   60 FORMAT(3X,'T12',8X,'T22',8X,'V11A',7X,'V12A',7X,'V22A'/)
      PRINT 50, T12,T22,V11A,V12A,V22A
      PRINT 70
   70 FORMAT(3X,4HV11B,7X,4HV12B,7X,4HV22B,7X,'V1111',6X,'V2111'/)
      PRINT 50, V11B,V12B,V22B,V1111,V2111
      PRINT 80
   80 FORMAT(3X,5HV2121,6X,5HV2211,6X,5HV2221,6X,5HV2222/)
      PRINT 50, V2121,V2211,V2221,V2222
   90 RETURN
      END

C*********************************************************************
      FUNCTION F0(ARG)
C
C CALCULATES THE F FUNCTION
C FO ONLY (S-TYPE ORBITALS)
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DATA PI/3.1415926535898D0/
      IF (ARG.LT.1.0D-6) GO TO 10
C F0 IN TERMS OF THE ERROR FUNCTION
      F0=DSQRT(PI/ARG)*DERFOTHER(DSQRT(ARG))/2.0D0
      GO TO 20
C ASYMPTOTIC VALUE FOR SMALL ARGUMENTS
   10 F0=1.0D0-ARG/3.0D0
   20 CONTINUE
      RETURN
      END

C*********************************************************************
      FUNCTION DERFOTHER(ARG)
C
C CALCULATES THE ERROR FUNCTION ACCORDING TO A RATIONAL
C APPROXIMATION FROM M. ARBRAMOWITZ AND I.A. STEGUN,
C HANDBOOK OF MATHEMATICAL FUNCTIONS, DOVER.
C ABSOLUTE ERROR IS LESS THAN 1.5*10**(-7)
C CAN BE REPLACED BY A BUILT-IN FUNCTION ON SOME MACHINES
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(5)
      DATA P/0.3275911D0/
      DATA A/0.254829592D0,-0.284496736D0,1.421413741D0,
     $ -1.453152027D0,1.061405429D0/
      T=1.0D0/(1.0D0+P*ARG)
      TN=T
      POLY=A(1)*TN
      DO 10 I=2,5
      TN=TN*T
      POLY=POLY+A(I)*TN
   10 CONTINUE
      DERFOTHER=1.0D0-POLY*DEXP(-ARG*ARG)
      RETURN
      END

C*********************************************************************
      FUNCTION S(A,B,RAB2)
C
C CALCULATES OVERLAPS FOR UN-NORMALIZED PRIMITIVES
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DATA PI/3.1415926535898D0/
      S=(PI/(A+B))**1.5D0*DEXP(-A*B*RAB2/(A+B))
      RETURN
      END

C*********************************************************************
      FUNCTION T(A,B,RAB2)
C
C CALCULATES KINETIC ENERGY INTEGRALS FOR UN-NORMALIZED PRIMITIVES
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DATA PI/3.1415926535898D0/
      T=A*B/(A+B)*(3.0D0-2.0D0*A*B*RAB2/(A+B))*(PI/(A+B))**1.5D0
     $ *DEXP(-A*B*RAB2/(A+B))
      RETURN
      END

C*********************************************************************
      FUNCTION V(A,B,RAB2,RCP2,ZC)
C
C CALCULATES UN-NORMALIZED NUCLEAR ATTRACTION INTEGRALS
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DATA PI/3.1415926535898D0/
      V=2.0D0*PI/(A+B)*F0((A+B)*RCP2)*DEXP(-A*B*RAB2/(A+B))
      V=-V*ZC
      RETURN
      END

C*********************************************************************
      FUNCTION TWOE(A,B,C,D,RAB2,RCD2,RPQ2)
C
C CALCULATES TWO-ELECTRON INTEGRALS FOR UN-NORMALIZED PRIMITIVES
C A,B,C,D ARE THE EXPONENTS ALPHA, BETA, ETC.
C RAB2 EQUALS SQUARED DISTANCE BETWEEN CENTER A AND CENTER B, ETC.
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DATA PI/3.1415926535898D0/
      TWOE=2.0D0*(PI**2.5D0)/((A+B)*(C+D)*DSQRT(A+B+C+D))
     $ *F0((A+B)*(C+D)*RPQ2/(A+B+C+D))
     $ *DEXP(-A*B*RAB2/(A+B)-C*D*RCD2/(C+D))
      RETURN
      END

C*********************************************************************
      SUBROUTINE COLECT(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C
C THIS TAKES THE BASIC INTEGRALS FROM COMMON AND ASSEMBLES THE
C RELEVENT MATRICES, THAT IS S,H,X,XT, AND TWO-ELECTRON INTEGRALS
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/MATRIX/S(2,2),X(2,2),XT(2,2),H(2,2),F(2,2),G(2,2),C(2,2),
     $ FPRIME(2,2),CPRIME(2,2),P(2,2),OLDP(2,2),TT(2,2,2,2),E(2,2)
      COMMON/INT/S12,T11,T12,T22,V11A,V12A,V22A,V11B,V12B,V22B,
     $ V1111,V2111,V2121,V2211,V2221,V2222
C FORM CORE HAMILTONIAN
      H(1,1)=T11+V11A+V11B
      H(1,2)=T12+V12A+V12B
      H(2,1)=H(1,2)
      H(2,2)=T22+V22A+V22B
C FORM OVERLAP MATRIX
      S(1,1)=1.0D0
      S(1,2)=S12
      S(2,1)=S(1,2)
      S(2,2)=1.0D0
C USE CANONICAL ORTHOGONALIZATION
      X(1,1)=1.0D0/DSQRT(2.0D0*(1.0D0+S12))
      X(2,1)=X(1,1)
      X(1,2)=1.0D0/DSQRT(2.0D0*(1.0D0-S12))
      X(2,2)=-X(1,2)
C TRANSPOSE OF TRANSFORMATION MATRIX
      XT(1,1)=X(1,1)
      XT(1,2)=X(2,1)
      XT(2,1)=X(1,2)
      XT(2,2)=X(2,2)
C MATRIX OF TWO-ELE�CTRON INTEGRALS
      TT(1,1,1,1)=V1111
      TT(2,1,1,1)=V2111
      TT(1,2,1,1)=V2111
      TT(1,1,2,1)=V2111
      TT(1,1,1,2)=V2111
      TT(2,1,2,1)=V2121
      TT(1,2,2,1)=V2121
      TT(2,1,1,2)=V2121
      TT(1,2,1,2)=V2121
      TT(2,2,1,1)=V2211
      TT(1,1,2,2)=V2211
      TT(2,2,2,1)=V2221
      TT(2,2,1,2)=V2221
      TT(2,1,2,2)=V2221
      TT(1,2,2,2)=V2221
      TT(2,2,2,2)=V2222
      IF (IOP.EQ.0) GO TO 40
      CALL MATOUT(S,2,2,2,2,4HS   )
      CALL MATOUT(X,2,2,2,2,4HX   )
      CALL MATOUT(H,2,2,2,2,4HH   )
      PRINT 10
   10 FORMAT(//)
      DO 30 I=1,2
      DO 30 J=1,2
      DO 30 K=1,2
      DO 30 L=1,2
      PRINT 20, I,J,K,L,TT(I,J,K,L)
   20 FORMAT(3X,1H(,4I2,2H ),F10.6)
   30 CONTINUE
   40 RETURN
      END

C*********************************************************************
      SUBROUTINE SCF(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C
C PERFORMS THE SCF ITERATIONS
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/MATRIX/S(2,2),X(2,2),XT(2,2),H(2,2),F(2,2),G(2,2),C(2,2),
     $ FPRIME(2,2),CPRIME(2,2),P(2,2),OLDP(2,2),TT(2,2,2,2),E(2,2)
      DATA PI/3.1415926535898D0/
C CONVERGENCE CRITERION FOR DENSITY MATRIX
      DATA CRIT/1.0D-4/
C MAXIMUM NUMBER OF ITERATIONS
      DATA MAXIT/25/
C ITERATION NUMBER
      ITER=0
C USE CORE-HAMILTONIAN FOR INITIAL GUESS AT F, I.E. (P=0)
      DO 10 I=1,2
      DO 10 J=1,2
   10 P(I,J)=0.0D0
      IF (IOP.LT.2) GO TO 20
      CALL MATOUT(P,2,2,2,2,4HP   )
C START OF ITERATION LOOP
   20 ITER=ITER+1
      IF (IOP.LT.2) GO TO 40
      PRINT 30, ITER
   30 FORMAT(/,4X,28HSTART OF ITERATION NUMBER = ,I2)
   40 CONTINUE
C FORM TWO-ELECTRON PART OF FOCK MATRIX FROM P
      CALL FORMG
      IF (IOP.LT.2) GO TO 50
      CALL MATOUT(G,2,2,2,2,4HG   )
   50 CONTINUE
C ADD CORE HAMILTONIAN TO GET FOCK MATRIX
      DO 60 I=1,2
      DO 60 J=1,2
      F(I,J) = H(I,J)+G(I,J)
   60 CONTINUE
C CALCULATE ELECTRONIC ENERGY
      EN=0.0D0
      DO 70 I=1,2
      DO 70 J=1,2
      EN=EN+0.5D0*P(I,J)*(H(I,J)+F(I,J))
   70 CONTINUE
      IF (IOP.LT.2) GO TO 90
      CALL MATOUT(F,2,2,2,2,4HF   )
      PRINT 80, EN
   80 FORMAT(///,4X,20HELECTRONIC ENERGY = ,D20.12)
   90 CONTINUE
C TRANSFORM FOCK MATRIX USING G FOR TEMPORARY STORAGE
      CALL MULT(F,X,G,2,2)
      CALL MULT(XT,G,FPRIME,2,2)
C DIAGONALIZE TRANSFORMED FOCK MATRIX
      CALL DIAG(FPRIME,CPRIME,E)
C TRANSFORM EIGENVECTORS TO GET MATRIX C
      CALL MULT(X,CPRIME,C,2,2)
C FORM NEW DENSITY MATRIX
      DO 100 I=1,2
      DO 100 J=1,2
C SAVE PRESENT DENSITY MATRIX
C BEFORE CREATING NEW ONE
      OLDP(I,J)=P(I,J)
      P(I,J)=0.0D0
      DO 100 K=1,1
      P(I,J)=P(I,J)+2.0D0*C(I,K)*C(J,K)
  100 CONTINUE
      IF (IOP.LT.2) GO TO 110
      CALL MATOUT(FPRIME,2,2,2,2,"F'  ")
      CALL MATOUT(CPRIME,2,2,2,2,"C'  ")
      CALL MATOUT(E,2,2,2,2,'E   ')
      CALL MATOUT(C,2,2,2,2,'C   ')
      CALL MATOUT(P,2,2,2,2,'P   ')
  110 CONTINUE
C CALCULATE DELTA
      DELTA=0.0D0
      DO 120 I=1,2
      DO 120 J=1,2
      DELTA=DELTA+(P(I,J)-OLDP(I,J))**2
  120 CONTINUE
      DELTA=DSQRT(DELTA/4.0D0)
      IF (IOP.EQ.0) GO TO 140
      PRINT 130, DELTA
  130 FORMAT(/,4X,39HDELTA(CONVERGENCE OF DENSITY MATRIX) =
     $F10.6,/)
  140 CONTINUE
C CHECK FOR CONVERGENCE
      IF (DELTA.LT.CRIT) GO TO 160
C NOT YET CONVERGED
C TEST FOR MAXIMUM NUMBER OF ITERATIONS
C IF MAXIMUM NUMBER NOT YET REACHED
C GO BACK FOR ANOTHER ITERATION
      IF(ITER.LT.MAXIT) GO TO 20
C SOMETHING WRONG HERE
      PRINT 150
  150 FORMAT(4X,21HNO CONVERGENCE IN SCF)
      STOP
  160 CONTINUE
C CALCULATION CONVERGED IF IT GOT HERE
C ADD NUCLEAR REPULSION TO GET TOTAL ENERGY
      ENT=EN+ZA*ZB/R
      IF (IOP.EQ.0) GO TO 180
      PRINT 170, EN, ENT
  170 FORMAT(//,4X,21HCALCULATION CONVERGED,//,
     $4X,20HELECTRONIC ENERGY = ,D20.12,//,
     $4X,20HTOTAL ENERGY =      ,D20.12   )
  180 CONTINUE
      IF (IOP.NE.1) GO TO 190
C PRINT OUT THE FINAL RESULTS IF
C HAVE NOT DONE SO ALREADY
      CALL MATOUT(G,2,2,2,2,4HG   )
      CALL MATOUT(F,2,2,2,2,4HF   )
      CALL MATOUT(E,2,2,2,2,4HE   )
      CALL MATOUT(C,2,2,2,2,4HC   )
      CALL MATOUT(P,2,2,2,2,4HP   )
  190 CONTINUE
C PS MATRIX HAS MULLIKEN POPULATIONS
      CALL MULT(P,S,OLDP,2,2)
      IF(IOP.EQ.0) GO TO 200
      CALL MATOUT(OLDP,2,2,2,2,4HPS   )
  200 CONTINUE
      RETURN
      END

C*********************************************************************
      SUBROUTINE FORMG
C
C CALCULATES THE G MATRIX FROM THE DENSITY MATRIX
C AND TWO-ELECTRON INTEGRALS
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON/MATRIX/S(2,2),X(2,2),XT(2,2),H(2,2),F(2,2),G(2,2),C(2,2),
     $FPRIME(2,2),CPRIME(2,2),P(2,2),OLDP(2,2),TT(2,2,2,2),E(2,2)
      DO 10 I=1,2
      DO 10 J=1,2
      G(I,J)=0.0D0
      DO 10 K=1,2
      DO 10 L=1,2
      G(I,J)=G(I,J)+P(K,L)*(TT(I,J,K,L)-0.5D0*TT(I,L,K,J))
   10 CONTINUE
      RETURN
      END

C*********************************************************************
      SUBROUTINE DIAG(F,C,E)
C
C DIAGONALIZES F TO GIVE EIGENVECTORS IN C AND EIGENVALUES IN E
C THETA IS THE ANGLE DESCRIBING SOLUTION
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION F(2,2),C(2,2),E(2,2)
      DATA PI/3.1415926535898D0/
      IF (DABS(F(1,1)-F(2,2)).GT.1.0D-20) GO TO 10
C HERE IS SYMMETRY DETERMINED SOLUTION (HOMONUCLEAR DIATOMIC)
      THETA=PI/4.0D0
      GO TO 20
   10 CONTINUE
C SOLUTION FOR HETERONUCLEAR DIATOMIC
      THETA=0.5D0*DATAN(2.0D0*F(1,2)/(F(1,1)-F(2,2)))
   20 CONTINUE
      C(1,1)=DCOS(THETA)
      C(2,1)=DSIN(THETA)
      C(1,2)=DSIN(THETA)
      C(2,2)=-DCOS(THETA)
      E(1,1)=F(1,1)*DCOS(THETA)**2+F(2,2)*DSIN(THETA)**2
     $ +F(1,2)*DSIN(2.0D0*THETA)
      E(2,2)=F(2,2)*DCOS(THETA)**2+F(1,1)*DSIN(THETA)**2
     $ -F(1,2)*DSIN(2.0D0*THETA)
      E(2,1)=0.0D0
      E(1,2)=0.0D0
C ORDER EIGENVALUES AND EIGENVECTORS
      IF (E(2,2).GT.E(1,1)) GO TO 30
      TEMP=E(2,2)
      E(2,2)=E(1,1)
      E(1,1)=TEMP
      TEMP=C(1,2)
      C(1,2)=C(1,1)
      C(1,1)=TEMP
      TEMP=C(2,2)
      C(2,2)=C(2,1)
      C(2,1)=TEMP
   30 RETURN
      END

C*********************************************************************
      SUBROUTINE MULT(A,B,C,IM,M)
C
C MULTIPLIES TWO SQUARE MATRICES A AND B TO GET C
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(IM,IM),B(IM,IM),C(IM,IM)
      DO 10 I=1,M
      DO 10 J=1,M
      C(I,J)=0.0D0
      DO 10 K=1,M
   10 C(I,J)=C(I,J)+A(I,K)*B(K,J)
      RETURN
      END

C*********************************************************************
      SUBROUTINE MATOUT(A,IM,IN,M,N,LABEL)
C
C PRINT MATRICES OF SIZE M BY N
C
C*********************************************************************

      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(IM,IN)
      IHIGH=0
   10 LOW=IHIGH+1
      IHIGH=IHIGH+5
      IHIGH=MIN(IHIGH,N)
      PRINT 20, LABEL,(I,I=LOW,IHIGH)
   20 FORMAT(///,3X,5H THE ,A4,6H ARRAY,/,15X,5(10X,I3,6X)//)
      DO 30 I=1,M
   30 PRINT 40, I,(A(I,J),J=LOW,IHIGH)
   40 FORMAT(I10,5X,5(1X,D18.10))
      IF (N-IHIGH) 50,50,10
   50 RETURN
      END