In one of the books there is an ancient Fortran program that has exactly what im looking for but i can't get it to compile. The book says that the original program was compiled using the FORTRAN IV compiler on a PDP-10. The program was subsequently "updated" by another author at some point.
after spending the better part of a day trying to get this thing to work I have finally decided to call on the greater internet to figure out how to compile this thing.
(this is from modern quantum chemistry by atila szabo and neil s. ostlund)
C*********************************************************************
C
C MINIMAL BASIS STO-03 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*********************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IOP=2
N=1
R=1.4632D0
ZETA1=2.0925D0
ZETA2=1.24D0
ZA=2.0D0
ZB=1.0D0
OPEN(UNIT=10, FILE='HeH.out')
CALL HFCALC(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
CLOSE(10)
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
WRITE(10,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 TEH 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 COEFFICEIENTS 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
WRITE(10, 40)
40 FORMAT(3X,'R',10X,'ZETA1',6X,'ZETA2',6X,'S12',8X,'T11'/)
WRITE(10, 50) R,ZETA1,ZETA2,S12,T11
50 FORMAT(5F11.6//)
WRITE(10, 60)
60 FORMAT(3X,'T12',8X,'T22',8X,'V11A',7X,'V12A',7X,'V22A'/)
WRITE(10, 50) T12,T22,V11A,V12A,V22A
WRITE(10, 70)
70 FORMAT(3X,4HV11B,7X,4HV12B,7X,4HV22B,7X,'V1111',6X,'V2111'/)
WRITE(10, 50) V11B,V12B,V22B,V1111,V2111
WRITE(10, 80)
80 FORMAT(3X,5HV2121,6X,5HV2211,6X,5HV2221,6X,5HV2222/)
WRITE(10, 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)*DERF(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 DERF(ARG)
C
C CALCULATES TEH 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
DERF=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-ELECTRON 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 )
WRITE(10, 10)
10 FORMAT(//)
DO 30 I=1,2
DO 30 J=1,2
DO 30 K=1,2
DO 30 L=1,2
WRITE(10, 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
WRITE(10, 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 )
WRITE(10, 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
WRITE(10,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
WRITE(10, 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
WRITE(10, 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)
WRITE(10, 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 WRITE(10, 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