博文
MINIMAL BASIS STO-3G CALCULATION ON HEH+
||
PROGRAM STO3G
C
C MINIMAL BASIS STO-3G CALCULATION ON HEH+
C The main work was done by my senior sister apprentice(Honghui Shang) ,and I've just
C modified the SCF subroutine so that it will be convenient to calculate
C other system with more than two orbitals .
IMPLICIT INTEGER(I-N)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IOP=2
N=3
R=1.4632
ZETA1=2.0925
ZETA2=1.24
ZA=2.0
ZB=1.0
CALL HFCALC(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
END
C************************************************************************
SUBROUTINE HFCALC(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C
C DOES A HARTREE-FOCK CALCULATION FOR A TWO-ELECTRON DIATOMIC
C USING THE 1S MINIMAL STO-NG BASIS SET
C MINIMAL BASIS SET HAS BASIS FUNCTIONS 1 & 2 ON NUCLEI A & B
C
C IOP=0 NO PRINTING AT ALL(ie. OPTIMIZE EXPONENTS)
C IOP=1 PRINT ONLY CONVERGED RESULTS
C IOP=2 PRINT EVERY ITERATION
C N 1,2 OR 3
C R BOND LENGTH
C ZETA1 SLATER ORBITAL (EXPONENT FUNCTION 1)
C ZETA2 SLATER ORBITAL (EXPONENT FUNCTION 2)
C ZA ATOMIC # (ATOM A)
C ZB ATOMIC # (ATOM B)
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IMPLICIT INTEGER(I-N)
IF (IOP.EQ.0) GO TO 20
WRITE(*,10) N,ZA,ZB
10 FORMAT(2X,'STO-',I1,'G FOR ATOMIC NUMBERS ',F5.2,' AND ',
$ F5.2,//)
20 CONTINUE
C CALC. ALL ONE & TWO ELECTRON INTEGRALS
CALL INTGRL(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C PLACE ALL INTEGRALS IN ARRAYS
CALL COLECT(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C DO THE SCF CALC.
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 BASIC INTEGRALS NEEDED FOR SCF CALC.
C
C************************************************************************
IMPLICIT INTEGER(I-N)
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.1415926535898/
C THESE ARE THE CONTRACTION COEFICIENTS AND EXPONENTS FOR A
C NORMALIZED 1S SLATER ORBITAL WITH EXPONENT 1.0 IN TERMS OF
C NORMALIZED 1S PRIMITIVE GAUSSIANS
DATA COEF,EXPON/1.0,2*0.0,0.678914,0.430129,0.0,
$ 0.444635,0.535328,0.154329,0.270950,2*0.0,0.151623,
$ 0.851819,0.0,0.109818,0.405771,2.22766/
R2=R*R
C SCALE THE EXPONENTS (A) OF PRIMITIVE GAUSSIANS
C INCLUDE NORMALIZATION IN CONTRACTION COEFFICIENTS (D)
DO I=1,N
A1(I)=EXPON(I,N)*(ZETA1**2)
D1(I)=COEF(I,N)*((2.0*A1(I)/PI)**0.75)
A2(I)=EXPON(I,N)*(ZETA2**2)
D2(I)=COEF(I,N)*((2.0*A2(I)/PI)**0.75)
ENDDO
C D & A ARE NOW THE CONTRACTION COEFFICIENTS AND EXPONENTS
C IN TERMS OF UNNORMALIZED PRIMITIVE GAUSSIANS
S12=0.0
T11=0.0
T12=0.0
T22=0.0
V11A=0.0
V12A=0.0
V22A=0.0
V11B=0.0
V12B=0.0
V22B=0.0
V1111=0.0
V2111=0.0
V2121=0.0
V2211=0.0
V2221=0.0
V2222=0.0
C CALC ONE-ELECTRON INTEGRALS
C CENTER A IS FIRST ATOM, CENTER B IS SECOND ATOM
C ORIGIN IS ON CENTER A
C V12A=OFF-DIAGONAL NUCLEAR ATTRACTION TO CENTER, ETC.
DO I=1,N
DO 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)
ENDDO
ENDDO
C CALC TWO-ELECTRON INTEGRALS
DO I=1,N
DO J=1,N
DO K=1,N
DO 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)
ENDDO
ENDDO
ENDDO
ENDDO
IF(IOP.EQ.0) GO TO 90
WRITE(*,40)
40 FORMAT(3X,'R',10X,'ZETA1',6X,'ZETA2',6X,'S12',8X,'T11'/)
WRITE(*,50) R,ZETA1,ZETA2,S12,T11
50 FORMAT(5F11.6//)
WRITE(*,60)
60 FORMAT(3X,'T12',8X,'T22',8X,'V11A',7X,'V12A',7X,'V22A'/)
WRITE(*,50) T12,T22,V11A,V12A,V22A
WRITE(*,70)
70 FORMAT(3X,'V11B',7X,'V12B',7X,'V221',7X,'V1111',6X,'V2111'/)
WRITE(*,50) V11B,V12B,V22B,V1111,V2111
WRITE(*,80)
80 FORMAT(3X,'V2121',6X,'V2211',6X,'V2221',6X,'V2222'/)
WRITE(*,50) V2121,V2211,V2221,V2222
90 RETURN
END
C************************************************************************
FUNCTION FO(ARG)
C
C CALC. THE F FUNCTION
C FO ONLY (S-TYPE ORBITALS)
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IMPLICIT INTEGER(I-N)
DATA PI/3.1415926535898/
IF(ARG.LT.1.0D-6) GO TO 10
C FO IN TERMS OF THE ERROR FUNCTION,Here DERF is an intrinsic function
FO=SQRT(PI/ARG)*DERF(SQRT(ARG))/2.0
GO TO 20
C ASYMPTOTIC VALUE FOR SMALL ARGUMENTS
10 FO=1.0-ARG/3.0
20 CONTINUE
RETURN
END
C************************************************************************
C************************************************************************
FUNCTION S(A,B,RAB2)
C
C CALC. OVERLAPS FOR UN-NORMALIZED PRIMITIVES
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IMPLICIT INTEGER(I-N)
DATA PI/3.1415926535898/
S=(PI/(A+B))**1.5*EXP(-A*B*RAB2/(A+B))
RETURN
END
C************************************************************************
FUNCTION T(A,B,RAB2)
C
C CALC. KINETIC ENERGY INTEGRALS FOR UN-NORMALIZED PRIMITIVE
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DATA PI/3.1415926535898/
T=A*B/(A+B)*(3.0-2.0*A*B*RAB2/(A+B))*(PI/(A+B))**1.5
$ *EXP(-A*B*RAB2/(A+B))
RETURN
END
C************************************************************************
FUNCTION V(A,B,RAB2,RCP2,ZC)
C
C CALC. UN-NORMALIZED NUCLEAR ATTRACTION INTEGRALS
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DATA PI/3.1415926535898/
V=2.0*PI/(A+B)*FO((A+B)*RCP2)*EXP(-A*B*RAB2/(A+B))
V=-V*ZC
RETURN
END
C************************************************************************
FUNCTION TWOE(A,B,C,D,RAB2,RCD2,RPQ2)
C
C CALC. TWO-ELECTRON INTEGRALS FOR UN-NORMALIZED PRIMITIVES
C A,B,C,D ARE THE EXPONENTS ALPHA,BETA, ETC.
C RAB2 IS SQUARED DISTANCE BETWEEN CENTER A & CENTER B, ETC.
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DATA PI/3.1415926535898/
TWOE=2.0*(PI**2.5)/((A+B)*(C+D)*SQRT(A+B+C+D))
$ *FO((A+B)*(C+D)*RPQ2/(A+B+C+D))
$ *EXP(-A*B*RAB2/(A+B)-C*D*RCD2/(C+D))
RETURN
END
C**********************************************************************
SUBROUTINE COLECT(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C
C TAKES BASIC INTEGRALS FROM COMMON AND ASSEMBLES THE RELEVANT
C MATRICES--S,H,X,XT & TWO-ELECTRON INTEGRALS.
C
C**********************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IMPLICIT INTEGER(I-N)
COMMON/MATRIX/S(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 FROM 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.0
S(1,2)=S12
S(2,1)=S(1,2)
S(2,2)=1.0
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,'S ')
CALL MATOUT(H,2,2,2,2,'H ')
WRITE(*,10)
10 FORMAT(//)
DO I=1,2
DO J=1,2
DO K=1,2
DO L=1,2
WRITE(6,20) I,J,K,L,TT(I,J,K,L)
20 FORMAT(' ',3X,'(',4I2,' )',F10.6)
ENDDO
ENDDO
ENDDO
ENDDO
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)
IMPLICIT INTEGER(I-N)
real(kind=8)::L(2,2),L_inv(2,2),L_inv_T(2,2),F_in(2,2),Eig(2)
integer,parameter::dim=2
COMMON/MATRIX/S(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.1415926535898/
C CONVERGENCE CRITERION FOR DENSITY MATRIX
DATA CRIT/1.0D-6/
C MAX. NO. ITERATIONS
DATA MAXIT/100/
C ITERATION NUMBER
ITER=0
C USE CORE-HAMILTONIAN FOR INITIAL GUESS AT F I.E. (P=0)
DO I=1,2
DO J=1,2
P(I,J)=0.0
ENDDO
ENDDO
IF(IOP.LT.2) GO TO 20
CALL MATOUT(P,2,2,2,2,'P ')
C START OF ITERATION LOOP
20 ITER=ITER+1
IF(IOP.LT.2) GO TO 40
WRITE(*,30) ITER
30 FORMAT(/,4X,'START 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,'G ')
50 CONTINUE
C ADD CORE HAMILTONIAN TO GET FOCK MATRIX
DO I=1,2
DO J=1,2
F(I,J)=H(I,J)+G(I,J)
ENDDO
ENDDO
C CALC. ELECTRONIC ENERGY
EN=0.0
DO I=1,2
DO J=1,2
EN=EN+0.5*P(I,J)*(H(I,J)+F(I,J))
ENDDO
ENDDO
IF(IOP.LT.2) GO TO 90
CALL MATOUT(F,2,2,2,2,'F ')
WRITE(*,80) EN
80 FORMAT(///,4X,'ELECTRONIC ENERGY = ',D20.12)
90 CONTINUE
C Decompose the S matrix
call Chol(S,L,L_inv,dim)
C Calc the FPRIME and C PRIME Matrix
L_inv_T=TRANSPOSE(L_inv)
FPRIME=MATMUL(MATMUL(L_inv,F),L_inv_T)
C Calc the eigenvalue and eigenvector for Equation Fprime*Cprime=Cprime*E
C Fprime=L_inv*F*L_inv' and Cprime=L'*C
F_in=FPRIME
CPRIME=MATMUL(TRANSPOSE(L),C)
call DSYEVJ(F_in,CPRIME,E)
C calc the C matrix
C=matmul(L_inv_T,CPRIME)
C FORM NEW DENSITY MATRIX
DO I=1,2
DO J=1,2
C SAVE PRESENT DENSITY MATRIX BEFORE CREATING NEW ONE
OLDP(I,J)=P(I,J)
P(I,J)=0.0
DO K=1,1
P(I,J)=P(I,J)+2.0*C(I,K)*C(J,K)
ENDDO
ENDDO
ENDDO
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 CALC. DELTA
DELTA=0.0
DO I=1,2
DO J=1,2
DELTA=DELTA+(P(I,J)-OLDP(I,J))**2
ENDDO
ENDDO
DELTA=SQRT(DELTA/4.0)
IF(IOP.EQ.0) GO TO 140
WRITE(*,130) DELTA
130 FORMAT(/,4X,'DELTA(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 MAX. NUMBER NOT YET REACHED THEN GO BACK FOR
C ANOTHER ITERATION
IF(ITER.LT.MAXIT) GO TO 20
C SOMETHING WRONG HERE
WRITE(*,150)
150 FORMAT(4X,'NO 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(*,170)EN,ENT
170 FORMAT(//,4X,'CALCULATION CONVERGED',//,
$ 4X, 'ELECTRONIC ENERGY = ',D20.12,//,
$ 4X, 'TOTAL ENERGY = ' ,D20.12)
180 CONTINUE
IF(IOP.NE.1) GO TO 190
C PRINT OUT THE FINAL RESULTS IF NOT ALREADY DONE SO
CALL MATOUT(G,2,2,2,2,'G ')
CALL MATOUT(F,2,2,2,2,'F ')
CALL MATOUT(E,2,2,2,2,'E ')
CALL MATOUT(C,2,2,2,2,'C ')
CALL MATOUT(P,2,2,2,2,'P ')
190 CONTINUE
C PS MATRIX HAS MULIKEN POPULATIONS
OLDP=MATMUL(P,S)
IF(IOP.EQ.0) GO TO 200
CALL MATOUT(OLDP,2,2,2,2,'PS ')
200 CONTINUE
RETURN
END
C***********************************************************************
SUBROUTINE FORMG
C
C CALC THE G MATRIX FROM THE DENSITY MATRIX & TWO-ELECTRON
C INTEGRALS
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IMPLICIT INTEGER(I-N)
COMMON/MATRIX/S(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.0
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.5*TT(I,L,K,J))
10 CONTINUE
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)
IMPLICIT INTEGER(I-N)
Character*2 Label
DIMENSION A(IM,IN)
IHIGH=0
10 LOW=IHIGH+1
IHIGH=IHIGH+5
IHIGH=MIN(IHIGH,N)
WRITE(6,*)' '
WRITE(6,20) LABEL
20 FORMAT(' ',3X,'THE ',A,' ARRAY')
WRITE(6,*)' ', (I,I=LOW,IHIGH)
DO 30 I=1,M
30 WRITE(*,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
SUBROUTINE DSYEVJ(A, Q, W)
* ************************************************
* Calculates the eigenvalues and normalized eigenvectors of a symmetric
* matrix A using the Jacobi algorithm.
* The upper triangular part of A is destroyed during the calculation,
* the diagonal elements are read but not destroyed, and the lower
* triangular elements are not referenced at all.
* ----------------------------------------------------------------------------
* Parameters:
* A: The symmetric input matrix
* Q: Storage buffer for eigenvectors
* W: Storage buffer for eigenvalues
* ----------------------------------------------------------------------------
* .. Arguments ..
DOUBLE PRECISION A(2,2)
DOUBLE PRECISION Q(2,2)
DOUBLE PRECISION W(2,2)
* .. Parameters ..
INTEGER N
PARAMETER ( N = 2 )
* .. Local Variables ..
DOUBLE PRECISION SD, SO
DOUBLE PRECISION S, C, T
DOUBLE PRECISION G, H, Z, THETA
DOUBLE PRECISION THRESH
INTEGER I, X, Y, R
* Initialize Q to the identitity matrix
* --- This loop can be omitted if only the eigenvalues are desired ---
DO 10 X = 1, N
Q(X,X) = 1.0D0
DO 11, Y = 1, X-1
Q(X, Y) = 0.0D0
Q(Y, X) = 0.0D0
11 CONTINUE
10 CONTINUE
* Initialize W to diag(A)
DO X = 1, N
Do Y=1,N
if(X==Y)then
W(X,Y) = A(X, X)
else
W(X,Y)=0
end if
end do
end do
* Calculate SQR(tr(A))
SD = 0.0D0
DO 30 X = 1, N
SD = SD + ABS(W(X,X))
30 CONTINUE
SD = SD**2
* Main iteration loop
DO 40 I = 1, 50
* Test for convergence
SO = 0.0D0
DO 50 X = 1, N
DO 51 Y = X+1, N
SO = SO + ABS(A(X, Y))
51 CONTINUE
50 CONTINUE
IF (SO .EQ. 0.0D0) THEN
RETURN
END
PROGRAM STO3G
C
C MINIMAL BASIS STO-3G CALCULATION ON HEH+
C The main work was done by my senior sister apprentice ,and I've just
C modified the SCF subroutine so that it will be convenient to calculate
C other system with more than two orbitals .
IMPLICIT INTEGER(I-N)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IOP=2
N=3
R=1.4632
ZETA1=2.0925
ZETA2=1.24
ZA=2.0
ZB=1.0
CALL HFCALC(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
END
C************************************************************************
SUBROUTINE HFCALC(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C
C DOES A HARTREE-FOCK CALCULATION FOR A TWO-ELECTRON DIATOMIC
C USING THE 1S MINIMAL STO-NG BASIS SET
C MINIMAL BASIS SET HAS BASIS FUNCTIONS 1 & 2 ON NUCLEI A & B
C
C IOP=0 NO PRINTING AT ALL(ie. OPTIMIZE EXPONENTS)
C IOP=1 PRINT ONLY CONVERGED RESULTS
C IOP=2 PRINT EVERY ITERATION
C N 1,2 OR 3
C R BOND LENGTH
C ZETA1 SLATER ORBITAL (EXPONENT FUNCTION 1)
C ZETA2 SLATER ORBITAL (EXPONENT FUNCTION 2)
C ZA ATOMIC # (ATOM A)
C ZB ATOMIC # (ATOM B)
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IMPLICIT INTEGER(I-N)
IF (IOP.EQ.0) GO TO 20
WRITE(*,10) N,ZA,ZB
10 FORMAT(2X,'STO-',I1,'G FOR ATOMIC NUMBERS ',F5.2,' AND ',
$ F5.2,//)
20 CONTINUE
C CALC. ALL ONE & TWO ELECTRON INTEGRALS
CALL INTGRL(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C PLACE ALL INTEGRALS IN ARRAYS
CALL COLECT(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C DO THE SCF CALC.
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 BASIC INTEGRALS NEEDED FOR SCF CALC.
C
C************************************************************************
IMPLICIT INTEGER(I-N)
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.1415926535898/
C THESE ARE THE CONTRACTION COEFICIENTS AND EXPONENTS FOR A
C NORMALIZED 1S SLATER ORBITAL WITH EXPONENT 1.0 IN TERMS OF
C NORMALIZED 1S PRIMITIVE GAUSSIANS
DATA COEF,EXPON/1.0,2*0.0,0.678914,0.430129,0.0,
$ 0.444635,0.535328,0.154329,0.270950,2*0.0,0.151623,
$ 0.851819,0.0,0.109818,0.405771,2.22766/
R2=R*R
C SCALE THE EXPONENTS (A) OF PRIMITIVE GAUSSIANS
C INCLUDE NORMALIZATION IN CONTRACTION COEFFICIENTS (D)
DO I=1,N
A1(I)=EXPON(I,N)*(ZETA1**2)
D1(I)=COEF(I,N)*((2.0*A1(I)/PI)**0.75)
A2(I)=EXPON(I,N)*(ZETA2**2)
D2(I)=COEF(I,N)*((2.0*A2(I)/PI)**0.75)
ENDDO
C D & A ARE NOW THE CONTRACTION COEFFICIENTS AND EXPONENTS
C IN TERMS OF UNNORMALIZED PRIMITIVE GAUSSIANS
S12=0.0
T11=0.0
T12=0.0
T22=0.0
V11A=0.0
V12A=0.0
V22A=0.0
V11B=0.0
V12B=0.0
V22B=0.0
V1111=0.0
V2111=0.0
V2121=0.0
V2211=0.0
V2221=0.0
V2222=0.0
C CALC ONE-ELECTRON INTEGRALS
C CENTER A IS FIRST ATOM, CENTER B IS SECOND ATOM
C ORIGIN IS ON CENTER A
C V12A=OFF-DIAGONAL NUCLEAR ATTRACTION TO CENTER, ETC.
DO I=1,N
DO 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)
ENDDO
ENDDO
C CALC TWO-ELECTRON INTEGRALS
DO I=1,N
DO J=1,N
DO K=1,N
DO 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)
ENDDO
ENDDO
ENDDO
ENDDO
IF(IOP.EQ.0) GO TO 90
WRITE(*,40)
40 FORMAT(3X,'R',10X,'ZETA1',6X,'ZETA2',6X,'S12',8X,'T11'/)
WRITE(*,50) R,ZETA1,ZETA2,S12,T11
50 FORMAT(5F11.6//)
WRITE(*,60)
60 FORMAT(3X,'T12',8X,'T22',8X,'V11A',7X,'V12A',7X,'V22A'/)
WRITE(*,50) T12,T22,V11A,V12A,V22A
WRITE(*,70)
70 FORMAT(3X,'V11B',7X,'V12B',7X,'V221',7X,'V1111',6X,'V2111'/)
WRITE(*,50) V11B,V12B,V22B,V1111,V2111
WRITE(*,80)
80 FORMAT(3X,'V2121',6X,'V2211',6X,'V2221',6X,'V2222'/)
WRITE(*,50) V2121,V2211,V2221,V2222
90 RETURN
END
C************************************************************************
FUNCTION FO(ARG)
C
C CALC. THE F FUNCTION
C FO ONLY (S-TYPE ORBITALS)
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IMPLICIT INTEGER(I-N)
DATA PI/3.1415926535898/
IF(ARG.LT.1.0D-6) GO TO 10
C FO IN TERMS OF THE ERROR FUNCTION,Here DERF is an intrinsic function
FO=SQRT(PI/ARG)*DERF(SQRT(ARG))/2.0
GO TO 20
C ASYMPTOTIC VALUE FOR SMALL ARGUMENTS
10 FO=1.0-ARG/3.0
20 CONTINUE
RETURN
END
C************************************************************************
C************************************************************************
FUNCTION S(A,B,RAB2)
C
C CALC. OVERLAPS FOR UN-NORMALIZED PRIMITIVES
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IMPLICIT INTEGER(I-N)
DATA PI/3.1415926535898/
S=(PI/(A+B))**1.5*EXP(-A*B*RAB2/(A+B))
RETURN
END
C************************************************************************
FUNCTION T(A,B,RAB2)
C
C CALC. KINETIC ENERGY INTEGRALS FOR UN-NORMALIZED PRIMITIVE
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DATA PI/3.1415926535898/
T=A*B/(A+B)*(3.0-2.0*A*B*RAB2/(A+B))*(PI/(A+B))**1.5
$ *EXP(-A*B*RAB2/(A+B))
RETURN
END
C************************************************************************
FUNCTION V(A,B,RAB2,RCP2,ZC)
C
C CALC. UN-NORMALIZED NUCLEAR ATTRACTION INTEGRALS
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DATA PI/3.1415926535898/
V=2.0*PI/(A+B)*FO((A+B)*RCP2)*EXP(-A*B*RAB2/(A+B))
V=-V*ZC
RETURN
END
C************************************************************************
FUNCTION TWOE(A,B,C,D,RAB2,RCD2,RPQ2)
C
C CALC. TWO-ELECTRON INTEGRALS FOR UN-NORMALIZED PRIMITIVES
C A,B,C,D ARE THE EXPONENTS ALPHA,BETA, ETC.
C RAB2 IS SQUARED DISTANCE BETWEEN CENTER A & CENTER B, ETC.
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
DATA PI/3.1415926535898/
TWOE=2.0*(PI**2.5)/((A+B)*(C+D)*SQRT(A+B+C+D))
$ *FO((A+B)*(C+D)*RPQ2/(A+B+C+D))
$ *EXP(-A*B*RAB2/(A+B)-C*D*RCD2/(C+D))
RETURN
END
C**********************************************************************
SUBROUTINE COLECT(IOP,N,R,ZETA1,ZETA2,ZA,ZB)
C
C TAKES BASIC INTEGRALS FROM COMMON AND ASSEMBLES THE RELEVANT
C MATRICES--S,H,X,XT & TWO-ELECTRON INTEGRALS.
C
C**********************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IMPLICIT INTEGER(I-N)
COMMON/MATRIX/S(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 FROM 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.0
S(1,2)=S12
S(2,1)=S(1,2)
S(2,2)=1.0
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,'S ')
CALL MATOUT(H,2,2,2,2,'H ')
WRITE(*,10)
10 FORMAT(//)
DO I=1,2
DO J=1,2
DO K=1,2
DO L=1,2
WRITE(6,20) I,J,K,L,TT(I,J,K,L)
20 FORMAT(' ',3X,'(',4I2,' )',F10.6)
ENDDO
ENDDO
ENDDO
ENDDO
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)
IMPLICIT INTEGER(I-N)
real(kind=8)::L(2,2),L_inv(2,2),L_inv_T(2,2),F_in(2,2),Eig(2)
integer,parameter::dim=2
COMMON/MATRIX/S(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.1415926535898/
C CONVERGENCE CRITERION FOR DENSITY MATRIX
DATA CRIT/1.0D-6/
C MAX. NO. ITERATIONS
DATA MAXIT/100/
C ITERATION NUMBER
ITER=0
C USE CORE-HAMILTONIAN FOR INITIAL GUESS AT F I.E. (P=0)
DO I=1,2
DO J=1,2
P(I,J)=0.0
ENDDO
ENDDO
IF(IOP.LT.2) GO TO 20
CALL MATOUT(P,2,2,2,2,'P ')
C START OF ITERATION LOOP
20 ITER=ITER+1
IF(IOP.LT.2) GO TO 40
WRITE(*,30) ITER
30 FORMAT(/,4X,'START 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,'G ')
50 CONTINUE
C ADD CORE HAMILTONIAN TO GET FOCK MATRIX
DO I=1,2
DO J=1,2
F(I,J)=H(I,J)+G(I,J)
ENDDO
ENDDO
C CALC. ELECTRONIC ENERGY
EN=0.0
DO I=1,2
DO J=1,2
EN=EN+0.5*P(I,J)*(H(I,J)+F(I,J))
ENDDO
ENDDO
IF(IOP.LT.2) GO TO 90
CALL MATOUT(F,2,2,2,2,'F ')
WRITE(*,80) EN
80 FORMAT(///,4X,'ELECTRONIC ENERGY = ',D20.12)
90 CONTINUE
C Decompose the S matrix
call Chol(S,L,L_inv,dim)
C Calc the FPRIME and C PRIME Matrix
L_inv_T=TRANSPOSE(L_inv)
FPRIME=MATMUL(MATMUL(L_inv,F),L_inv_T)
C Calc the eigenvalue and eigenvector for Equation Fprime*Cprime=Cprime*E
C Fprime=L_inv*F*L_inv' and Cprime=L'*C
F_in=FPRIME
CPRIME=MATMUL(TRANSPOSE(L),C)
call DSYEVJ(F_in,CPRIME,E)
C calc the C matrix
C=matmul(L_inv_T,CPRIME)
C FORM NEW DENSITY MATRIX
DO I=1,2
DO J=1,2
C SAVE PRESENT DENSITY MATRIX BEFORE CREATING NEW ONE
OLDP(I,J)=P(I,J)
P(I,J)=0.0
DO K=1,1
P(I,J)=P(I,J)+2.0*C(I,K)*C(J,K)
ENDDO
ENDDO
ENDDO
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 CALC. DELTA
DELTA=0.0
DO I=1,2
DO J=1,2
DELTA=DELTA+(P(I,J)-OLDP(I,J))**2
ENDDO
ENDDO
DELTA=SQRT(DELTA/4.0)
IF(IOP.EQ.0) GO TO 140
WRITE(*,130) DELTA
130 FORMAT(/,4X,'DELTA(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 MAX. NUMBER NOT YET REACHED THEN GO BACK FOR
C ANOTHER ITERATION
IF(ITER.LT.MAXIT) GO TO 20
C SOMETHING WRONG HERE
WRITE(*,150)
150 FORMAT(4X,'NO 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(*,170)EN,ENT
170 FORMAT(//,4X,'CALCULATION CONVERGED',//,
$ 4X, 'ELECTRONIC ENERGY = ',D20.12,//,
$ 4X, 'TOTAL ENERGY = ' ,D20.12)
180 CONTINUE
IF(IOP.NE.1) GO TO 190
C PRINT OUT THE FINAL RESULTS IF NOT ALREADY DONE SO
CALL MATOUT(G,2,2,2,2,'G ')
CALL MATOUT(F,2,2,2,2,'F ')
CALL MATOUT(E,2,2,2,2,'E ')
CALL MATOUT(C,2,2,2,2,'C ')
CALL MATOUT(P,2,2,2,2,'P ')
190 CONTINUE
C PS MATRIX HAS MULIKEN POPULATIONS
OLDP=MATMUL(P,S)
IF(IOP.EQ.0) GO TO 200
CALL MATOUT(OLDP,2,2,2,2,'PS ')
200 CONTINUE
RETURN
END
C***********************************************************************
SUBROUTINE FORMG
C
C CALC THE G MATRIX FROM THE DENSITY MATRIX & TWO-ELECTRON
C INTEGRALS
C
C************************************************************************
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
IMPLICIT INTEGER(I-N)
COMMON/MATRIX/S(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.0
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.5*TT(I,L,K,J))
10 CONTINUE
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)
IMPLICIT INTEGER(I-N)
Character*2 Label
DIMENSION A(IM,IN)
IHIGH=0
10 LOW=IHIGH+1
IHIGH=IHIGH+5
IHIGH=MIN(IHIGH,N)
WRITE(6,*)' '
WRITE(6,20) LABEL
20 FORMAT(' ',3X,'THE ',A,' ARRAY')
WRITE(6,*)' ', (I,I=LOW,IHIGH)
DO 30 I=1,M
30 WRITE(*,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
SUBROUTINE DSYEVJ(A, Q, W)
* ************************************************
* Calculates the eigenvalues and normalized eigenvectors of a symmetric
* matrix A using the Jacobi algorithm.
* The upper triangular part of A is destroyed during the calculation,
* the diagonal elements are read but not destroyed, and the lower
* triangular elements are not referenced at all.
* ----------------------------------------------------------------------------
* Parameters:
* A: The symmetric input matrix
* Q: Storage buffer for eigenvectors
* W: Storage buffer for eigenvalues
* ----------------------------------------------------------------------------
* .. Arguments ..
DOUBLE PRECISION A(2,2)
DOUBLE PRECISION Q(2,2)
DOUBLE PRECISION W(2,2)
* .. Parameters ..
INTEGER N
PARAMETER ( N = 2 )
* .. Local Variables ..
DOUBLE PRECISION SD, SO
DOUBLE PRECISION S, C, T
DOUBLE PRECISION G, H, Z, THETA
DOUBLE PRECISION THRESH
INTEGER I, X, Y, R
* Initialize Q to the identitity matrix
* --- This loop can be omitted if only the eigenvalues are desired ---
DO 10 X = 1, N
Q(X,X) = 1.0D0
DO 11, Y = 1, X-1
Q(X, Y) = 0.0D0
Q(Y, X) = 0.0D0
11 CONTINUE
10 CONTINUE
* Initialize W to diag(A)
DO X = 1, N
Do Y=1,N
if(X==Y)then
W(X,Y) = A(X, X)
else
W(X,Y)=0
end if
end do
end do
* Calculate SQR(tr(A))
SD = 0.0D0
DO 30 X = 1, N
SD = SD + ABS(W(X,X))
30 CONTINUE
SD = SD**2
* Main iteration loop
DO 40 I = 1, 50
* Test for convergence
SO = 0.0D0
DO 50 X = 1, N
DO 51 Y = X+1, N
SO = SO + ABS(A(X, Y))
51 CONTINUE
50 CONTINUE
IF (SO .EQ. 0.0D0) THEN
RETURN
END
https://blog.sciencenet.cn/blog-276702-829704.html
上一篇：Cholesky decomposition for positive definite matrix
下一篇：原子电子结构计算程序atom-0.1.0
