Personal computing discussed

Moderators: renee, SecretSquirrel, just brew it!

 
Stranger
Graphmaster Gerbil
Topic Author
Posts: 1433
Joined: Thu Mar 13, 2003 8:51 pm
Location: Socialist republic of Ohio
Contact:

Fortran will be the Death of me

Mon Sep 24, 2012 9:39 pm

I'm a chemist who has been picking up programing and recently i started a project implementing a restricted hartree fock solver in java using the COLT matrix libraries. so far everything has been going peachy in the java world. unfortunately I'm majorly lacking in any kind of test case.

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
 
just brew it!
Administrator
Posts: 54500
Joined: Tue Aug 20, 2002 10:51 pm
Location: Somewhere, having a beer

Re: Fortran will be the Death of me

Mon Sep 24, 2012 9:43 pm

What compiler are you trying to use?
Nostalgia isn't what it used to be.
 
just brew it!
Administrator
Posts: 54500
Joined: Tue Aug 20, 2002 10:51 pm
Location: Somewhere, having a beer

Re: Fortran will be the Death of me

Mon Sep 24, 2012 9:51 pm

Been a long time since I've used FORTRAN, but it looks like the whitespace has been mangled. IIRC comments (the lines starting with a 'C') need to begin in column 1, and non-comment statements need to begin in column 7. Statement labels (the numbers at the start of some lines) need to appear to the left of column 6.

Any way you can track down a copy of the file that isn't mangled like this?

Edit: If you want to try and fix the copy you've got:
1. Comments need to be shifted to the left so that the 'C' is in column 1.
2. Lines with numerical statement labels need to be shifted so that the statement itself starts in column 7 and the label appears to the left of column 6.
3. Lines that start with a '$' and a space appear to be continuation lines. Shift these over such that the '$' is in column 6.
4. Anything that ISN'T a comment, a statement with a numerical label, or a continuation needs to be shifted such that the statement starts in column 7.
Nostalgia isn't what it used to be.
 
JustAnEngineer
Gerbil God
Posts: 19673
Joined: Sat Jan 26, 2002 7:00 pm
Location: The Heart of Dixie

Re: Fortran will be the Death of me

Mon Sep 24, 2012 10:04 pm

Remember that FORTRAN IV and FORTRAN 77 were originally written on punchcards, so each line of code is column sensitive. A "C" in the first column indicates a comment. Numbers in columns 1-5 are statement numbers for jumps. A character in column 6 indicates that the current line is a continuation of the previous line.

Oops, JBI! beat me to it.
 
Stranger
Graphmaster Gerbil
Topic Author
Posts: 1433
Joined: Thu Mar 13, 2003 8:51 pm
Location: Socialist republic of Ohio
Contact:

Re: Fortran will be the Death of me

Mon Sep 24, 2012 10:35 pm

so is the first character supposed to be a space? by column to do mean something like a character?(i'm guessing this means there is a specific number of columns on a line) It unfortunately came to me mangled. I had tried to fix it with a tab. I found a random ancient fortran compiler g77 that im guessing someone compiled for windows in the early ninties. maybe i should switch to linux for this.
 
Ouroboros
Gerbil In Training
Posts: 8
Joined: Sat Jan 22, 2011 2:48 am

Re: Fortran will be the Death of me

Mon Sep 24, 2012 10:42 pm

You can try using GFortran on Linux. Since it accepts F90, you will likely need to convert the line continuations. See this:

http://people.sc.fsu.edu/~jburkardt/f_s ... ixcon.html
 
just brew it!
Administrator
Posts: 54500
Joined: Tue Aug 20, 2002 10:51 pm
Location: Somewhere, having a beer

Re: Fortran will be the Death of me

Mon Sep 24, 2012 10:45 pm

Yes, by column I mean the position of the character on the line. First character position is 1, second is 2, etc.

The first character on each line should be a 'C' for the lines that are comments, otherwise it should be a space. All lines need to be positioned as I described in my previous post; use spaces to make things line up as indicated, do NOT use tabs.

As JAE notes, FORTRAN was originally designed for 80-column punch cards. Anything after the 72nd column is ignored. Wikipedia has a description of the FORTRAN line format: http://en.wikipedia.org/wiki/Fortran#Fixed_layout

Your old copy of g77 is probably OK for this; FORTRAN IV has been around since 1962! If you want a newer compiler you can switch to Linux, or install Cygwin if you'd prefer to stick with Windows.
Nostalgia isn't what it used to be.
 
PopcornMachine
Gerbil
Posts: 16
Joined: Tue Dec 02, 2008 11:25 pm

Re: Fortran will be the Death of me

Mon Sep 24, 2012 10:48 pm

I used to code Fortran a "few" years back.

Unless that's a very ancient Fortran, I'd say that it's actually written in very early BASIC.

All those gotos and do loops. Very ugly stuff.

One could rewrite it in a more modern language, but not a trivial task.
 
just brew it!
Administrator
Posts: 54500
Joined: Tue Aug 20, 2002 10:51 pm
Location: Somewhere, having a beer

Re: Fortran will be the Death of me

Mon Sep 24, 2012 10:52 pm

PopcornMachine wrote:
I used to code Fortran a "few" years back.

Unless that's a very ancient Fortran, I'd say that it's actually written in very early BASIC.

As noted in the OP, it is FORTRAN IV. FORTRAN IV has been around since 1962. My guess is that you were using FORTRAN 90, which does not look much like FORTRAN IV!
Nostalgia isn't what it used to be.
 
PopcornMachine
Gerbil
Posts: 16
Joined: Tue Dec 02, 2008 11:25 pm

Re: Fortran will be the Death of me

Mon Sep 24, 2012 10:58 pm

just brew it! wrote:
PopcornMachine wrote:
I used to code Fortran a "few" years back.

Unless that's a very ancient Fortran, I'd say that it's actually written in very early BASIC.

As noted in the OP, it is FORTRAN IV. FORTRAN IV has been around since 1962. My guess is that you were using FORTRAN 90, which does not look much like FORTRAN IV!


Just try to help. Sorry I hit a nerve. Looks more basic to me.

In either case, I don't see it getting complied now.

And by a few years back I meant in the 80s.
 
just brew it!
Administrator
Posts: 54500
Joined: Tue Aug 20, 2002 10:51 pm
Location: Somewhere, having a beer

Re: Fortran will be the Death of me

Mon Sep 24, 2012 11:22 pm

No nerve, sorry if it seemed like I had taken offense, that was not my intent!

If it was the 80s you were probably using FORTRAN 77.

Pretty sure the current GNU FORTRAN 77 compiler will still compile FORTRAN IV code (IIRC FORTRAN 77 is still backwards compatible with FORTRAN IV).
Nostalgia isn't what it used to be.
 
Stranger
Graphmaster Gerbil
Topic Author
Posts: 1433
Joined: Thu Mar 13, 2003 8:51 pm
Location: Socialist republic of Ohio
Contact:

Re: Fortran will be the Death of me

Tue Sep 25, 2012 1:46 am

alright after spending some personal time with the text editor I've got...

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 Laboriously 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


this now gives me......

c:\Fort99\G77\bin>
c:\Fort99\G77\bin>
c:\Fort99\G77\bin>g77.exe RHF.f
RHF.f:187: warning:
         F0=DSQRT(PI/ARG)*DERF(DSQRT(ARG))/2.0D0
                          1
RHF.f:196: (continued):
         FUNCTION DERF(ARG)
                  2
Same name `derf' used for global at (2) and intrinsic at (1) [info -f g77 M INTGLOB]
RHF.f: In subroutine `colect':
RHF.f:330:
         CALL MATOUT(S,2,2,2,2,4HS)
              1                    2
Invalid token at (2) in expression or subexpression at (1)
RHF.f:331:
         CALL MATOUT(X,2,2,2,2,4HX)
              1                    2
Invalid token at (2) in expression or subexpression at (1)
RHF.f:332:
         CALL MATOUT(H,2,2,2,2,4HH)
              1                    2
Invalid token at (2) in expression or subexpression at (1)
RHF.f: In subroutine `scf':
RHF.f:367:
         CALL MATOUT(P,2,2,2,2,4HP)
              1                    2
Invalid token at (2) in expression or subexpression at (1)
RHF.f:377:
         CALL MATOUT(G,2,2,2,2,4HG)
              1                    2
Invalid token at (2) in expression or subexpression at (1)
RHF.f:391:
         CALL MATOUT(F,2,2,2,2,4HF)
              1                    2
Invalid token at (2) in expression or subexpression at (1)
RHF.f:429: warning:
        $ F10.6,/)
          ^
Missing comma in FORMAT statement at (^)
RHF.f:450: warning:
        $ 4X,20HTOTAL ENERGY = ,D20.12)
                                    ^
Missing comma in FORMAT statement at (^)
RHF.f:450:
        $ 4X,20HTOTAL ENERGY = ,D20.12)
                                      ^
Invalid form for FORMAT statement at (^)
RHF.f:455:
         CALL MATOUT(G,2,2,2,2,4HG )
              1                     2
Invalid token at (2) in expression or subexpression at (1)
RHF.f:456:
         CALL MATOUT(F,2,2,2,2,4HF )
              1                     2
Invalid token at (2) in expression or subexpression at (1)
RHF.f:457:
         CALL MATOUT(E,2,2,2,2,4HE )
              1                     2
Invalid token at (2) in expression or subexpression at (1)
RHF.f:458:
         CALL MATOUT(C,2,2,2,2,4HC )
              1                     2
Invalid token at (2) in expression or subexpression at (1)
RHF.f:459:
         CALL MATOUT(P,2,2,2,2,4HP )
              1                     2
Invalid token at (2) in expression or subexpression at (1)
RHF.f:464:
         CALL MATOUT(OLDP,2,2,2,2,4HPS )
              1                         2
Invalid token at (2) in expression or subexpression at (1)
RHF.f: In subroutine `matout':
RHF.f:560:
         10 LOW=IHIGH+1
            ^
Invalid form for assignment statement at (^)
RHF.f:568:
         IF (N-IHIGH) 50,50,10
                            ^
Undefined label, first referenced at (^)


this error code in particular has me stumped...

         CALL MATOUT(F,2,2,2,2,4HF)
              1                    2
Invalid token at (2) in expression or subexpression at (1)
 
just brew it!
Administrator
Posts: 54500
Joined: Tue Aug 20, 2002 10:51 pm
Location: Somewhere, having a beer

Re: Fortran will be the Death of me

Tue Sep 25, 2012 2:01 am

Looks like the compiler provides its own implementation of the DERF() function, this is causing a conflict. Should probably just delete the implementation of this function from your program. (Note the comment in the code: "CAN BE REPLACED BY A BUILT-IN FUNCTION ON SOME MACHINES")

Some of the calls to MATOUT() are screwed up. Looks like the last argument is a string and should be quoted?

Looks like some of the FORMAT statements are screwed up as well. My FORTRAN's rather rusty, but it looks like maybe some strings are missing their enclosing quotes here as well, just like the last argument to the MATOUT() calls.

You appear to have incorrectly indented line 560 (this is also indirectly the cause of the error being reported on line 568).

Edit: Antique programming languages are fun, eh? :lol:
Nostalgia isn't what it used to be.
 
Stranger
Graphmaster Gerbil
Topic Author
Posts: 1433
Joined: Thu Mar 13, 2003 8:51 pm
Location: Socialist republic of Ohio
Contact:

Re: Fortran will be the Death of me

Tue Sep 25, 2012 2:29 am

I think this whole thing is going to drive me crazy : )


CALL MATOUT(S,2,2,2,2,4HS)

needed to be....

CALL MATOUT(S,2,2,2,2,4HS )

why this works I have no earthly idea...

So here's what I'm down to now...

c:\Fort99\G77\bin>g77.exe -Wall RHFV2.f
RHFV2.f:187: warning:
         F0=DSQRT(PI/ARG)*DERF(DSQRT(ARG))/2.0D0
                          1
RHFV2.f:196: (continued):
         FUNCTION DERF(ARG)
                  2
Same name `derf' used for global at (2) and intrinsic at (1) [info -f g77 M INTGLOB]
RHFV2.f: In subroutine `scf':
RHFV2.f:429: warning:
        $ F10.6,/)
          ^
Missing comma in FORMAT statement at (^)
RHFV2.f:355: warning: unused variable `pi'
C:\Users\Andrew\AppData\Local\Temp\cc2lcaaa.o(.text+0x1560):RHFV2.f: undefined reference to `erf'

c:\Fort99\G77\bin>
 
Stranger
Graphmaster Gerbil
Topic Author
Posts: 1433
Joined: Thu Mar 13, 2003 8:51 pm
Location: Socialist republic of Ohio
Contact:

Re: Fortran will be the Death of me

Tue Sep 25, 2012 2:32 am

renaming the error function DERF to DERFOTHER cuts the list down to....

c:\Fort99\G77\bin>g77.exe -Wall RHFV2.f
RHFV2.f: In subroutine `scf':
RHFV2.f:429: warning:
        $ F10.6,/)
          ^
Missing comma in FORMAT statement at (^)
RHFV2.f:355: warning: unused variable `pi'

c:\Fort99\G77\bin>


I think this means its at the point where it will produce a program.... i think.
 
Stranger
Graphmaster Gerbil
Topic Author
Posts: 1433
Joined: Thu Mar 13, 2003 8:51 pm
Location: Socialist republic of Ohio
Contact:

Re: Fortran will be the Death of me

Tue Sep 25, 2012 2:41 am

holy **** I'm a FORTRAN programer now. I feel like I went spelunking into the long lost depths of computing pasts.

    THE G    ARRAY
                           1                  2
         1        0.0000000000E+00   0.0000000000E+00
         2        0.0000000000E+00   0.0000000000E+00



    THE F    ARRAY
                           1                  2
         1       -0.2379140570E+01  -0.1277051444E+01
         2       -0.1277051444E+01  -0.1691423286E+01



    THE F' C ARRAY
                           1                  2
         1       -0.2325185985E+01  -0.3797836810E+00
         2       -0.3797836810E+00  -0.1317620358E+01



    THE C' E ARRAY
                           1                  2
         1        0.9482923002E+00   0.3173983512E+00
         2        0.3173983512E+00  -0.9482923002E+00



    THE E C  ARRAY
                           1                  2
         1       -0.2452301553E+01   0.0000000000E+00
         2        0.0000000000E+00  -0.1190504790E+01



    THE C P  ARRAY
                           1                  2
         1        0.8576680677E+00  -0.6958973545E+00
         2        0.2659508699E+00   0.1071978441E+01



    THE P E  ARRAY
                           1                  2
         1        0.1471189029E+01   0.4561951374E+00
         2        0.4561951374E+00   0.1414597304E+00



    THE G    ARRAY
                           1                  2
         1        0.1182093547E+01   0.3349608710E+00
         2        0.3349608710E+00   0.9512187968E+00



    THE F    ARRAY
                           1                  2
         1       -0.1197047024E+01  -0.9420905731E+00
         2       -0.9420905731E+00  -0.7402044894E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1341281306E+01  -0.2522858496E+00
         2       -0.2522858496E+00  -0.4611170162E-01



    THE C' E ARRAY
                           1                  2
         1        0.9827988313E+00   0.1846793359E+00
         2        0.1846793359E+00  -0.9827988313E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1388688752E+01   0.0000000000E+00
         2        0.0000000000E+00   0.1295745021E-02



    THE C P  ARRAY
                           1                  2
         1        0.7543989898E+00  -0.8066906511E+00
         2        0.4101062778E+00   0.1025514741E+01



    THE P E  ARRAY
                           1                  2
         1        0.1138235672E+01   0.6187675233E+00
         2        0.6187675233E+00   0.3363743182E+00



    THE G    ARRAY
                           1                  2
         1        0.1151879951E+01   0.2834392447E+00
         2        0.2834392447E+00   0.8960293942E+00



    THE F    ARRAY
                           1                  2
         1       -0.1227260619E+01  -0.9936121994E+00
         2       -0.9936121994E+00  -0.7953938921E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1407423890E+01  -0.2384932576E+00
         2       -0.2384932576E+00  -0.3078446362E-01



    THE C' E ARRAY
                           1                  2
         1        0.9861259982E+00   0.1659985414E+00
         2        0.1659985414E+00  -0.9861259982E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1447570416E+01   0.0000000000E+00
         2        0.0000000000E+00   0.9362062250E-02



    THE C P  ARRAY
                           1                  2
         1        0.7389570966E+00  -0.8208593374E+00
         2        0.4294904866E+00   0.1017548802E+01



    THE P E  ARRAY
                           1                  2
         1        0.1092115181E+01   0.6347500861E+00
         2        0.6347500861E+00   0.3689241562E+00



    THE G    ARRAY
                           1                  2
         1        0.1147660358E+01   0.2773268383E+00
         2        0.2773268383E+00   0.8885012439E+00



    THE F    ARRAY
                           1                  2
         1       -0.1231480213E+01  -0.9997246059E+00
         2       -0.9997246059E+00  -0.8029220424E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1415837997E+01  -0.2366661464E+00
         2       -0.2366661464E+00  -0.3036994843E-01



    THE C' E ARRAY
                           1                  2
         1        0.9864832833E+00   0.1638619289E+00
         2        0.1638619289E+00  -0.9864832833E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1455149937E+01   0.0000000000E+00
         2        0.0000000000E+00   0.8941991171E-02



    THE C P  ARRAY
                           1                  2
         1        0.7371771533E+00  -0.8224581978E+00
         2        0.4316937724E+00   0.1016616019E+01



    THE P E  ARRAY
                           1                  2
         1        0.1086860311E+01   0.6364695725E+00
         2        0.6364695725E+00   0.3727190263E+00



    THE G    ARRAY
                           1                  2
         1        0.1147179051E+01   0.2766463140E+00
         2        0.2766463140E+00   0.8876453143E+00



    THE F    ARRAY
                           1                  2
         1       -0.1231961519E+01  -0.1000405130E+01
         2       -0.1000405130E+01  -0.8037779720E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1416785065E+01  -0.2364592651E+00
         2       -0.2364592651E+00  -0.3034925810E-01



    THE C' E ARRAY
                           1                  2
         1        0.9865226529E+00   0.1636247392E+00
         2        0.1636247392E+00  -0.9865226529E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1456004221E+01   0.0000000000E+00
         2        0.0000000000E+00   0.8869897651E-02



    THE C P  ARRAY
                           1                  2
         1        0.7369793843E+00  -0.8226354171E+00
         2        0.4319381899E+00   0.1016512195E+01



    THE P E  ARRAY
                           1                  2
         1        0.1086277226E+01   0.6366590825E+00
         2        0.6366590825E+00   0.3731411998E+00



    THE G    ARRAY
                           1                  2
         1        0.1147125639E+01   0.2765710040E+00
         2        0.2765710040E+00   0.8875503626E+00



    THE F    ARRAY
                           1                  2
         1       -0.1232014932E+01  -0.1000480440E+01
         2       -0.1000480440E+01  -0.8038729237E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1416890006E+01  -0.2364363258E+00
         2       -0.2364363258E+00  -0.3034729831E-01



    THE C' E ARRAY
                           1                  2
         1        0.9865270049E+00   0.1635984984E+00
         2        0.1635984984E+00  -0.9865270049E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1456098895E+01   0.0000000000E+00
         2        0.0000000000E+00   0.8861590727E-02



    THE C P  ARRAY
                           1                  2
         1        0.7369575026E+00  -0.8226550199E+00
         2        0.4319652282E+00   0.1016500705E+01



    THE P E  ARRAY
                           1                  2
         1        0.1086212721E+01   0.6366800315E+00
         2        0.6366800315E+00   0.3731879167E+00



    THE PS   ARRAY
                           1                  2
         1        0.1356512402E+01   0.1097826820E+01
         2        0.7951153049E+00   0.6434875977E+00

c:\Fort99\G77\bin>
 
Stranger
Graphmaster Gerbil
Topic Author
Posts: 1433
Joined: Thu Mar 13, 2003 8:51 pm
Location: Socialist republic of Ohio
Contact:

Re: Fortran will be the Death of me

Tue Sep 25, 2012 2:51 am

To ensure that this program isn't lost to the world again I'm going to post the output

Source Code(compiled with g77 on windows 7)
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 Cleaned up and debugged again by Andrew Long (2012)
C Compiled with g77 on windows 7
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
      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 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
      PRINT 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//)
      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 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
      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-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   )
      PRINT 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)
      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
    


H-He+ STO-1G
    STO-1G FOR ATOMIC NUMBERS  2.00 AND  1.00


   R          ZETA1      ZETA2      S12        T11

   1.463200   2.092500   1.240000   0.424546   1.779555


   T12        T22        V11A       V12A       V22A

   0.219882   0.624919  -3.476243  -1.125477  -1.286345


   V11B       V12B       V22B       V1111      V2111

  -0.682453  -0.371456  -1.029998   1.229037   0.447815


   V2121      V2211      V2221      V2222

   0.182076   0.612241   0.291359   0.728318



    THE S    ARRAY
                           1                  2
         1        0.1000000000E+01   0.4245455607E+00
         2        0.4245455607E+00   0.1000000000E+01



    THE X    ARRAY
                           1                  2
         1        0.5924433518E+00   0.9321365336E+00
         2        0.5924433518E+00  -0.9321365336E+00



    THE H    ARRAY
                           1                  2
         1       -0.2379140570E+01  -0.1277051444E+01
         2       -0.1277051444E+01  -0.1691423286E+01






    THE P    ARRAY
                           1                  2
         1        0.0000000000E+00   0.0000000000E+00
         2        0.0000000000E+00   0.0000000000E+00



    THE G    ARRAY
                           1                  2
         1        0.0000000000E+00   0.0000000000E+00
         2        0.0000000000E+00   0.0000000000E+00



    THE F    ARRAY
                           1                  2
         1       -0.2379140570E+01  -0.1277051444E+01
         2       -0.1277051444E+01  -0.1691423286E+01



    THE F' C ARRAY
                           1                  2
         1       -0.2325185985E+01  -0.3797836810E+00
         2       -0.3797836810E+00  -0.1317620358E+01



    THE C' E ARRAY
                           1                  2
         1        0.9482923002E+00   0.3173983512E+00
         2        0.3173983512E+00  -0.9482923002E+00



    THE E C  ARRAY
                           1                  2
         1       -0.2452301553E+01   0.0000000000E+00
         2        0.0000000000E+00  -0.1190504790E+01



    THE C P  ARRAY
                           1                  2
         1        0.8576680677E+00  -0.6958973545E+00
         2        0.2659508699E+00   0.1071978441E+01



    THE P E  ARRAY
                           1                  2
         1        0.1471189029E+01   0.4561951374E+00
         2        0.4561951374E+00   0.1414597304E+00



    THE G    ARRAY
                           1                  2
         1        0.1182093547E+01   0.3349608710E+00
         2        0.3349608710E+00   0.9512187968E+00



    THE F    ARRAY
                           1                  2
         1       -0.1197047024E+01  -0.9420905731E+00
         2       -0.9420905731E+00  -0.7402044894E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1341281306E+01  -0.2522858496E+00
         2       -0.2522858496E+00  -0.4611170162E-01



    THE C' E ARRAY
                           1                  2
         1        0.9827988313E+00   0.1846793359E+00
         2        0.1846793359E+00  -0.9827988313E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1388688752E+01   0.0000000000E+00
         2        0.0000000000E+00   0.1295745021E-02



    THE C P  ARRAY
                           1                  2
         1        0.7543989898E+00  -0.8066906511E+00
         2        0.4101062778E+00   0.1025514741E+01



    THE P E  ARRAY
                           1                  2
         1        0.1138235672E+01   0.6187675233E+00
         2        0.6187675233E+00   0.3363743182E+00



    THE G    ARRAY
                           1                  2
         1        0.1151879951E+01   0.2834392447E+00
         2        0.2834392447E+00   0.8960293942E+00



    THE F    ARRAY
                           1                  2
         1       -0.1227260619E+01  -0.9936121994E+00
         2       -0.9936121994E+00  -0.7953938921E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1407423890E+01  -0.2384932576E+00
         2       -0.2384932576E+00  -0.3078446362E-01



    THE C' E ARRAY
                           1                  2
         1        0.9861259982E+00   0.1659985414E+00
         2        0.1659985414E+00  -0.9861259982E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1447570416E+01   0.0000000000E+00
         2        0.0000000000E+00   0.9362062250E-02



    THE C P  ARRAY
                           1                  2
         1        0.7389570966E+00  -0.8208593374E+00
         2        0.4294904866E+00   0.1017548802E+01



    THE P E  ARRAY
                           1                  2
         1        0.1092115181E+01   0.6347500861E+00
         2        0.6347500861E+00   0.3689241562E+00



    THE G    ARRAY
                           1                  2
         1        0.1147660358E+01   0.2773268383E+00
         2        0.2773268383E+00   0.8885012439E+00



    THE F    ARRAY
                           1                  2
         1       -0.1231480213E+01  -0.9997246059E+00
         2       -0.9997246059E+00  -0.8029220424E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1415837997E+01  -0.2366661464E+00
         2       -0.2366661464E+00  -0.3036994843E-01



    THE C' E ARRAY
                           1                  2
         1        0.9864832833E+00   0.1638619289E+00
         2        0.1638619289E+00  -0.9864832833E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1455149937E+01   0.0000000000E+00
         2        0.0000000000E+00   0.8941991171E-02



    THE C P  ARRAY
                           1                  2
         1        0.7371771533E+00  -0.8224581978E+00
         2        0.4316937724E+00   0.1016616019E+01



    THE P E  ARRAY
                           1                  2
         1        0.1086860311E+01   0.6364695725E+00
         2        0.6364695725E+00   0.3727190263E+00



    THE G    ARRAY
                           1                  2
         1        0.1147179051E+01   0.2766463140E+00
         2        0.2766463140E+00   0.8876453143E+00



    THE F    ARRAY
                           1                  2
         1       -0.1231961519E+01  -0.1000405130E+01
         2       -0.1000405130E+01  -0.8037779720E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1416785065E+01  -0.2364592651E+00
         2       -0.2364592651E+00  -0.3034925810E-01



    THE C' E ARRAY
                           1                  2
         1        0.9865226529E+00   0.1636247392E+00
         2        0.1636247392E+00  -0.9865226529E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1456004221E+01   0.0000000000E+00
         2        0.0000000000E+00   0.8869897651E-02



    THE C P  ARRAY
                           1                  2
         1        0.7369793843E+00  -0.8226354171E+00
         2        0.4319381899E+00   0.1016512195E+01



    THE P E  ARRAY
                           1                  2
         1        0.1086277226E+01   0.6366590825E+00
         2        0.6366590825E+00   0.3731411998E+00



    THE G    ARRAY
                           1                  2
         1        0.1147125639E+01   0.2765710040E+00
         2        0.2765710040E+00   0.8875503626E+00



    THE F    ARRAY
                           1                  2
         1       -0.1232014932E+01  -0.1000480440E+01
         2       -0.1000480440E+01  -0.8038729237E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1416890006E+01  -0.2364363258E+00
         2       -0.2364363258E+00  -0.3034729831E-01



    THE C' E ARRAY
                           1                  2
         1        0.9865270049E+00   0.1635984984E+00
         2        0.1635984984E+00  -0.9865270049E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1456098895E+01   0.0000000000E+00
         2        0.0000000000E+00   0.8861590727E-02



    THE C P  ARRAY
                           1                  2
         1        0.7369575026E+00  -0.8226550199E+00
         2        0.4319652282E+00   0.1016500705E+01



    THE P E  ARRAY
                           1                  2
         1        0.1086212721E+01   0.6366800315E+00
         2        0.6366800315E+00   0.3731879167E+00



    THE PS   ARRAY
                           1                  2
         1        0.1356512402E+01   0.1097826820E+01
         2        0.7951153049E+00   0.6434875977E+00
PS C:\>


STO-2G
   STO-2G FOR ATOMIC NUMBERS  2.00 AND  1.00


   R          ZETA1      ZETA2      S12        T11

   T12        T22        V11A       V12A       V22A

   0.158710   0.734883  -4.013825  -1.079394  -1.262417


   V11B       V12B       V22B       V1111      V2111

  -0.678042  -0.409128  -1.189281   1.301613   0.432134


   V2121      V2211      V2221      V2222

   0.174271   0.605340   0.308926   0.771326



    THE S    ARRAY
                           1                  2
         1        0.1000000000E+01   0.4441615663E+00
         2        0.4441615663E+00   0.1000000000E+01



    THE X    ARRAY
                           1                  2
         1        0.5884060246E+00   0.9484418647E+00
         2        0.5884060246E+00  -0.9484418647E+00



    THE H    ARRAY
                           1                  2
         1       -0.2599171798E+01  -0.1329811047E+01
         2       -0.1329811047E+01  -0.1716815310E+01






    THE P    ARRAY
                           1                  2
         1        0.0000000000E+00   0.0000000000E+00
         2        0.0000000000E+00   0.0000000000E+00



    THE G    ARRAY
                           1                  2
         1        0.0000000000E+00   0.0000000000E+00
         2        0.0000000000E+00   0.0000000000E+00



    THE F    ARRAY
                           1                  2
         1       -0.2599171798E+01  -0.1329811047E+01
         2       -0.1329811047E+01  -0.1716815310E+01



    THE F' C ARRAY
                           1                  2
         1       -0.2415106926E+01  -0.4924157210E+00
         2       -0.4924157210E+00  -0.1489969850E+01



    THE C' E ARRAY
                           1                  2
         1        0.9177887745E+00   0.3970689680E+00
         2        0.3970689680E+00  -0.9177887745E+00



    THE E C  ARRAY
                           1                  2
         1       -0.2628143964E+01   0.0000000000E+00
         2        0.0000000000E+00  -0.1276932811E+01



    THE C P  ARRAY
                           1                  2
         1        0.9166292767E+00  -0.6368315237E+00
         2        0.1634356117E+00   0.1104107070E+01



    THE P E  ARRAY
                           1                  2
         1        0.1680418462E+01   0.2996197332E+00
         2        0.2996197332E+00   0.5342239838E-01



    THE G    ARRAY
                           1                  2
         1        0.1250786900E+01   0.3589713450E+00
         2        0.3589713450E+00   0.9839638378E+00



    THE F    ARRAY
                           1                  2
         1       -0.1348384898E+01  -0.9708397019E+00
         2       -0.9708397019E+00  -0.7328514724E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1392820536E+01  -0.3435100663E+00
         2       -0.3435100663E+00  -0.1255373491E+00



    THE C' E ARRAY
                           1                  2
         1        0.9693102812E+00   0.2458405555E+00
         2        0.2458405555E+00  -0.9693102812E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1479943006E+01   0.0000000000E+00
         2        0.0000000000E+00  -0.3841487951E-01



    THE C P  ARRAY
                           1                  2
         1        0.8035134840E+00  -0.7746803867E+00
         2        0.3371825343E+00   0.1063988515E+01



    THE P E  ARRAY
                           1                  2
         1        0.1291267838E+01   0.5418614257E+00
         2        0.5418614257E+00   0.2273841229E+00



    THE G    ARRAY
                           1                  2
         1        0.1192353794E+01   0.2917635025E+00
         2        0.2917635025E+00   0.9242296635E+00



    THE F    ARRAY
                           1                  2
         1       -0.1406818004E+01  -0.1038047544E+01
         2       -0.1038047544E+01  -0.7925856467E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1480270227E+01  -0.3427839803E+00
         2       -0.3427839803E+00  -0.1109212268E+00



    THE C' E ARRAY
                           1                  2
         1        0.9731890091E+00   0.2300068531E+00
         2        0.2300068531E+00  -0.9731890091E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1561284977E+01   0.0000000000E+00
         2        0.0000000000E+00  -0.2990647646E-01



    THE C P  ARRAY
                           1                  2
         1        0.7907784047E+00  -0.7876757805E+00
         2        0.3544821474E+00   0.1058350617E+01



    THE P E  ARRAY
                           1                  2
         1        0.1250660971E+01   0.5606336540E+00
         2        0.5606336540E+00   0.2513151856E+00



    THE G    ARRAY
                           1                  2
         1        0.1186439886E+01   0.2859115437E+00
         2        0.2859115437E+00   0.9182155638E+00



    THE F    ARRAY
                           1                  2
         1       -0.1412731912E+01  -0.1043899503E+01
         2       -0.1043899503E+01  -0.7985997465E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1488452111E+01  -0.3427280667E+00
         2       -0.3427280667E+00  -0.1111228056E+00



    THE C' E ARRAY
                           1                  2
         1        0.9734628514E+00   0.2288450938E+00
         2        0.2288450938E+00  -0.9734628514E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1569021839E+01   0.0000000000E+00
         2        0.0000000000E+00  -0.3055307817E-01



    THE C P  ARRAY
                           1                  2
         1        0.7898376740E+00  -0.7886190901E+00
         2        0.3557451389E+00   0.1057926754E+01



    THE P E  ARRAY
                           1                  2
         1        0.1247687102E+01   0.5619618261E+00
         2        0.5619618261E+00   0.2531092078E+00



    THE G    ARRAY
                           1                  2
         1        0.1186008092E+01   0.2854912932E+00
         2        0.2854912932E+00   0.9177766872E+00



    THE F    ARRAY
                           1                  2
         1       -0.1413163706E+01  -0.1044319754E+01
         2       -0.1044319754E+01  -0.7990386230E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1489044556E+01  -0.3427241141E+00
         2       -0.3427241141E+00  -0.1111499442E+00



    THE C' E ARRAY
                           1                  2
         1        0.9734821006E+00   0.2287631960E+00
         2        0.2287631960E+00  -0.9734821006E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1569582928E+01   0.0000000000E+00
         2        0.0000000000E+00  -0.3061157206E-01



    THE C P  ARRAY
                           1                  2
         1        0.7897713249E+00  -0.7886855360E+00
         2        0.3558341406E+00   0.1057896821E+01



    THE P E  ARRAY
                           1                  2
         1        0.1247477491E+01   0.5620552014E+00
         2        0.5620552014E+00   0.2532358713E+00



    THE G    ARRAY
                           1                  2
         1        0.1185977664E+01   0.2854617149E+00
         2        0.2854617149E+00   0.9177457614E+00



    THE F    ARRAY
                           1                  2
         1       -0.1413194134E+01  -0.1044349332E+01
         2       -0.1044349332E+01  -0.7990695488E+00



    THE F' C ARRAY
                           1                  2
         1       -0.1489086279E+01  -0.3427238362E+00
         2       -0.3427238362E+00  -0.1111519207E+00



    THE C' E ARRAY
                           1                  2
         1        0.9734834533E+00   0.2287574397E+00
         2        0.2287574397E+00  -0.9734834533E+00



    THE E C  ARRAY
                           1                  2
         1       -0.1569622448E+01   0.0000000000E+00
         2        0.0000000000E+00  -0.3061575235E-01



    THE C P  ARRAY
                           1                  2
         1        0.7897666614E+00  -0.7886902060E+00
         2        0.3558403961E+00   0.1057894717E+01



    THE P E  ARRAY
                           1                  2
         1        0.1247462759E+01   0.5620617632E+00
         2        0.5620617632E+00   0.2532447749E+00



    THE PS   ARRAY
                           1                  2
         1        0.1497108992E+01   0.1116136776E+01
         2        0.6745433591E+00   0.5028910080E+00
Last edited by Stranger on Tue Sep 25, 2012 3:24 am, edited 1 time in total.
 
just brew it!
Administrator
Posts: 54500
Joined: Tue Aug 20, 2002 10:51 pm
Location: Somewhere, having a beer

Re: Fortran will be the Death of me

Tue Sep 25, 2012 3:03 am

Stranger wrote:
I feel like I went spelunking into the long lost depths of computing pasts.

You did! And lived to tell the tale... :lol:

(And it looks like we got you to the finish line before anyone even replied to your post over on the Ubuntu forums... :wink:)
Nostalgia isn't what it used to be.
 
trackerben
Minister of Gerbil Affairs
Posts: 2188
Joined: Mon Jun 15, 2009 12:28 am
Location: 'Tween oceans...

Re: Fortran will be the Death of me

Tue Sep 25, 2012 3:49 am

TechReport = TechSupport

It's no wonder I come here a lot. Great articles, helpful folks, literate posts (well, mostly)
 
Glorious
Gerbilus Supremus
Posts: 12343
Joined: Tue Aug 27, 2002 6:35 pm

Re: Fortran will be the Death of me

Tue Sep 25, 2012 6:12 am

Sorry I didn't see this sooner, I program in Fortran 77 every day (DEC compiler though). Glad to see it was sorted out.

Man, do some of you guys ever sleep?
 
JustAnEngineer
Gerbil God
Posts: 19673
Joined: Sat Jan 26, 2002 7:00 pm
Location: The Heart of Dixie

Re: Fortran will be the Death of me

Tue Sep 25, 2012 6:21 am

I learned FORTRAN on a VAX 11/780 some 28 years ago. We were the first group at our university that didn't have to submit our code through JCL for batch processing on the IBM mainframe. :) Real-time compiling with immediate error messages was way better than waiting 10 minutes to 4 hours for your results to come back from the mainframe telling you that you'd left out a space. :lol:

Your smart phone today has hundreds of times the processing power of the system that we routinely had 100 simultaneous users logged into back then.
 
rootbear
Gerbil In Training
Posts: 5
Joined: Thu Sep 26, 2002 8:17 pm
Location: Bowie, MD (near Washington, DC)

Re: Fortran will be the Death of me

Tue Sep 25, 2012 8:27 am

Stranger wrote:
I think this whole thing is going to drive me crazy : )

CALL MATOUT(S,2,2,2,2,4HS)

needed to be....

CALL MATOUT(S,2,2,2,2,4HS )

why this works I have no earthly idea...


Just so you'll know, the 4HSbbb (where b is a blank) construct is a Hollerith constant, which is how strings were represented in old FORTRAN. Basically, it's an integer count (the length of the string), then an H, then the string itself. For example:

13HHello, World!

Surely the ugliest string syntax ever invented, but it made sense on the very, very simple machines on which FORTRAN was developed. Reading over this code gave me flash backs to when I was a co-op student at the Naval Research Lab and helped maintain their FORTRAN code (or Fortran, these days). Fortran 77 was the first version that had modern control structures like IF-THEN-ELSE, but the code I worked on at NRL was before that and looked a lot like the code posted here. I had a hard time explaining to the scientists I worked with why this new language I was learning at college, Pascal, was so much nicer because it used 'structured programming'!

I work at NASA now and there is still a ton of Fortran around, but in much more modern dialects. Aside from the quirky syntax, it's really not that much different from C these days.
 
just brew it!
Administrator
Posts: 54500
Joined: Tue Aug 20, 2002 10:51 pm
Location: Somewhere, having a beer

Re: Fortran will be the Death of me

Tue Sep 25, 2012 9:09 am

Jeez, I had completely forgotten about Hollerith constants. I'd forgotten about arithmetic IF as well, but remembered it when I saw it in the code.

Last time I actually used FORTRAN was probably around 1990...
Nostalgia isn't what it used to be.
 
bthylafh
Maximum Gerbil
Posts: 4320
Joined: Mon Dec 29, 2003 11:55 pm
Location: Southwest Missouri, USA

Re: Fortran will be the Death of me

Tue Sep 25, 2012 9:14 am

just brew it! wrote:
Jeez, I had completely forgotten about Hollerith constants. I'd forgotten about arithmetic IF as well, but remembered it when I saw it in the code.


How much alcohol did you have to consume to achieve this? :D
Hakkaa päälle!
i7-8700K|Asus Z-370 Pro|32GB DDR4|Asus Radeon RX-580|Samsung 960 EVO 1TB|1988 Model M||Logitech MX 518 & F310|Samsung C24FG70|Dell 2209WA|ATH-M50x
 
just brew it!
Administrator
Posts: 54500
Joined: Tue Aug 20, 2002 10:51 pm
Location: Somewhere, having a beer

Re: Fortran will be the Death of me

Tue Sep 25, 2012 12:39 pm

bthylafh wrote:
just brew it! wrote:
Jeez, I had completely forgotten about Hollerith constants. I'd forgotten about arithmetic IF as well, but remembered it when I saw it in the code.

How much alcohol did you have to consume to achieve this? :D

Only a little... :wink:
Nostalgia isn't what it used to be.
 
notfred
Maximum Gerbil
Posts: 4610
Joined: Tue Aug 10, 2004 10:10 am
Location: Ottawa, Canada

Re: Fortran will be the Death of me

Tue Sep 25, 2012 1:17 pm

My first job was programming in FORTRAN 77 in 1993, this thread is bringing back a bunch of stuff I'd forgotten as well.
 
pattakosn
Gerbil In Training
Posts: 1
Joined: Tue Sep 25, 2012 1:35 pm

Re: Fortran will be the Death of me

Tue Sep 25, 2012 1:45 pm

Do yourself a favor and convert it to free form fortran. this can be easily achieved by substituting all "C" in C style comments with "!" . Free style means that you no longer need to follow the 1,2-6,7,8-72 column rules and it will make your life easier.

FYI, a common source of problems for novices compiling fortran codes is which version of fortran the compiler thinks the source is written in. gfortran and ifort use the file's ending to determine it. You can either control it with compiler switces or by simply renaming the source file.

http://gcc.gnu.org/onlinedocs/gfortran/GNU-Fortran-and-GCC.html

If I may give my 2 cents of advice, always write using the IMPLICIT NONE statement and never write/modify fixed form source code.

There is no need to upset yourself, fortran is very reasonable :D, just google your problem.
 
Stranger
Graphmaster Gerbil
Topic Author
Posts: 1433
Joined: Thu Mar 13, 2003 8:51 pm
Location: Socialist republic of Ohio
Contact:

Re: Fortran will be the Death of me

Tue Sep 25, 2012 3:30 pm

Thanks for the help every one : ) I would have been going in circles for hours without it : )
 
Kalium
Gerbil In Training
Posts: 2
Joined: Fri Aug 16, 2013 9:37 am
Location: Italy
Contact:

Re: Fortran will be the Death of me

Fri Aug 16, 2013 9:48 am

I've found some errors that prevent the printing of some information (i.e. write(10 statement instead of print). The corrected code is posted below:
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 Compiled with g77 on windows 7
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
      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 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
      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
         


P.S.: The program converges on the same numbers reported in this thread. However, in the original book of Szabo the reported output is somewhat different.
Mine compiler was GNU Fortran (GCC) 4.8.0 20120319 (experimental) [trunk revision 185521]
 
Waco
Maximum Gerbil
Posts: 4850
Joined: Tue Jan 20, 2009 4:14 pm
Location: Los Alamos, NM

Re: Fortran will be the Death of me

Fri Aug 16, 2013 10:20 am

Debugging Fortran isn't fun.

Debugging parallel Fortran is worse. :cry:
Victory requires no explanation. Defeat allows none.

Who is online

Users browsing this forum: No registered users and 1 guest
GZIP: On