loos/eDFT_FUEG
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
master
eDFT_FUEG/FCI/fci.f

1723 lines
48 KiB
FortranFixed
Raw Permalink Normal View History

First commit
2019-06-16 22:35:10 +02:00
PROGRAM FCI
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
COMMON /CHAMIL/ ISYM,NNCI
COMMON /TAPES / INP,IOUT
DIMENSION MULTT(8,8)
DATA MULTT /1,2,3,4,5,6,7,8, 2,1,4,3,6,5,8,7,
2 3,4,1,2,7,8,5,6, 4,3,2,1,8,7,6,5,
4 5,6,7,8,1,2,3,4, 6,5,8,7,2,1,4,3,
6 7,8,5,6,3,4,1,2, 8,7,6,5,4,3,2,1/
INP = 5
IOUT = 6
WRITE (IOUT,35)
35 FORMAT(' PROGRAM * FCI (Determinant based Full CI)',
> 5X,'Author: P.J. Knowles, 1984')
C... SET UP DEFAULT OR INITIAL VALUES
MEMORY = 100 000
IPRTIM = -1
IPRINT = -1
IPRDIA = -1
IPRHAM = -1
IOPTIO = 0
MAXIT = 25
THR = 1D-5
THRRES = 0.05
NROOT = 1
NACT = 0
NELEC = 0
MS2 = 0
ISYM = 1
DO 2 I=1,8
NOCC(I) = 0
C... D2H MULTIPLICATION TABLE
DO 2 J=1,8
2 MULT(J,I) = MULTT(J,I)
DO 3 I=1,256
3 ITYPEA(I)=1
C
C... READ INPUT
CALL INPDAT
C
C... INITIALIZE MEMORY HANDLING
I = INICOR(MEMORY)
C
C... NUMBERS OF ALPHA, BETA ELECTRONS
NA = (NELEC+MS2)/2
NB = (NELEC-MS2)/2
DO 36 I=1,NACT
36 NOCC(ITYPEA(I)) = NOCC(ITYPEA(I))+1
CALL INITC
C
WRITE (IOUT,25) ISYM
25 FORMAT(/' Space symmetry:',T30,I3)
WRITE (IOUT,21) MAXIT,THR,THRRES,NROOT
21 FORMAT(
4' Maximum iterations:', T30,I3/
5' Convergence threshold:',T30,F12.7/
8' Output threshold:', T30,F12.7/
A' Number of roots:', T30,I3)
NNCI = NCI(ISYM)
IV1 = ICORR(NNCI)
IV2 = ICORR(NNCI)
CALL DAVIDS (Q(IV1),Q(IV2),IOPTIO)
CALL CORLSR(IV1)
C
CALL ENDCOR
END
SUBROUTINE FCMXM (D,Z,E,NK,NIJ,NOPCNT)
CSUBR MATRIX MULTIPLICATION KERNEL FOR FULL CI PROGRAM
CSUBR E = D * Z
CSUBR Z IS INTEGRAL MATRIX AND NOT SPARSE, DIMENSION NIJ=NUMBER OF
CSUBR ORBITAL PAIRS OF GIVEN SYMMETRY
CSUBR D,E LEADING DIMENSION NK=NUMBER OF DETERMINANTS OF A SYMMETRY
CSUBR BLOCK, OR SUB-BATCH OF THIS
CSUBR D IS VERY SPARSE IN 1ST ITERATION OF CI, GRADUALLY RISES TO
CSUBR USUALLY 50-80% POPULATED
CSUBR OPTIMUM FORM OF THIS ROUTINE IS DEPENDENT ON MACHINE ARCHITECTURE
CSUBR USUALLY NK>NIJ, SO BETTER TO VECTORISE ALONG NK, BUT SOMETIMES
CSUBR BETTER TO EXPLOIT SPARSITY OF D
CSUBR ROUTINE SHOULD RETURN NOPCNT -- NUMBER OF FLOATING POINT OPS --
CSUBR IF ADDITIONAL TIMING INFORMATION IS REQUIRED
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NHALF=30)
DIMENSION D(NK,NIJ),Z(NIJ,NIJ),E(NK,NIJ)
c pmwg replaced all code by a dgemm in Jan 2014
call dgemm('n','n',nk,nij,nij,1d0,d,nk,z,nij,0d0,e,nk)
c DO 1 J=1,NK*NIJ
c1 E(J,1) = 0D0
C... EVALUATE SPARSITY OF D
c IDS = 0
c DO 7 J=1,NK*NIJ
c7 IF (D(J,1).NE.0D0) IDS=IDS+1
c IZS = NIJ**2
c IF (IDS*(NIJ+NHALF).LT.IZS*(NK+NHALF)) THEN
C... USE SPARSITY OF D
c NOPCNT = NOPCNT + 2*IDS*NIJ
C CALL MXMA (Z,1,NIJ, D,NK,1, E,NK,1, NIJ,NIJ,NK)
c DO 2 K=1,NK
c DO 2 J=1,NIJ
c IF (D(K,J).EQ.0D0) GOTO 2
c DO 21 I=1,NIJ
c21 E(K,I) = E(K,I) + D(K,J)*Z(J,I)
c2 CONTINUE
c ELSE
C... VECTORISE ALONG D
c NOPCNT = NOPCNT + 2*IZS*NK
C CALL MXMA (D,1,NK, Z,NIJ,1, E,1,NK, NK,NIJ,NIJ)
c DO 3 I=1,NIJ
c DO 3 J=1,NIJ
C IF (Z(J,I).EQ.0D0) GOTO 3
c DO 31 K=1,NK
c31 E(K,I) = E(K,I) + D(K,J)*Z(J,I)
c3 CONTINUE
c END IF
RETURN
END
SUBROUTINE INPDAT
CSUBR DATA INPUT FOR PORTABLE FULL CI PROGRAM
CSUBR NAMELIST INPUT MAY REQUIRE ADJUSTMENT ON SOME MACHINES
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
COMMON /CHAMIL/ ISYM,NNCI
COMMON /TAPES / INP,IOUT
INTEGER ORBSYM(255)
EQUIVALENCE (ORBSYM(1),ITYPEA(1))
NAMELIST /FCI / NORB,NELEC,MS2,ISYM,ORBSYM
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM
> ,INT,MEMORY,CORE
> ,MAXIT,THR,THRRES,NROOT
INT = -1
READ (INP,FCI)
IF (INT.NE.-1) INP=INT
NACT = NORB
RETURN
END
SUBROUTINE INPINT (Z1,Z2)
CSUBR INTEGRAL INPUT FOR PORTABLE FULL CI PROGRAM
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
COMMON /TAPES / INP,IOUT
DIMENSION Z2(*),Z1(*)
cpmg
c rs = 2.0d0
c R = sqrt(dfloat(na+nb))/2 * rs
R = 0.25d0
write(*,*) 'Scaling integrals for R =',R
cpmg
CALL FZERO (Z1,NPAIR(1))
CALL FZERO (Z2,NINT2)
1 READ (INP,*,END=99) Z,I,J,K,L
IF (MAX(I,J,K,L).GT.NACT) GOTO 1
IF (K.NE.0) THEN
C... TWO ELECTRON INTEGRAL
cpmg
z = z/R
cpmg
ISYMIJ = MULT(ITYPEA(I),ITYPEA(J))
IJ = IQ(INTADT+I-1+(J-1)*NACT)
KL = IQ(INTADT+K-1+(L-1)*NACT)
Z2(IJ+(KL-1)*NPAIR(ISYMIJ)+INTOFF(ISYMIJ)) = Z
Z2(KL+(IJ-1)*NPAIR(ISYMIJ)+INTOFF(ISYMIJ)) = Z
ELSE IF (I.NE.0) THEN
C... ONE ELECTRON INTEGRAL
cpmg
z = z/R/R
cpmg
IJ = IQ(INTADT+I-1+(J-1)*NACT)
Z1(IJ) = Z
ELSE
C... CORE REPULSION ENERGY
CORE = Z
END IF
GOTO 1
99 RETURN
END
SUBROUTINE DAVIDS (V1,V2,IOPTIO)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /TAPES / INP,IOUT
C
C... DAVIDSON DIAGONALISER, SINGLE ROOT
C... WARNING, NOT NECESSARILY SAFE WITH DEGENERACIES, ETC.
PARAMETER (IFILD=92,IFIL1=93,IFIL2=94,MAXITR=200)
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
COMMON /CHAMIL/ ISYMH,N
DIMENSION V1(N),V2(N)
DIMENSION ZM(MAXITR,MAXITR),ZV(MAXITR,MAXITR),R(MAXITR)
C
TRIAL = 0
RESULT = 0
MAXITS = MIN(MAXIT,MAXITR,N)
CALL FCDIAG (V2)
IF (IOPTIO.GE.5) CALL OUTVEC (V2,N,'DIAGONAL ELEMENTS')
C
OPEN (IFILD,STATUS='SCRATCH',FORM='UNFORMATTED')
OPEN (IFIL1,STATUS='SCRATCH',FORM='UNFORMATTED')
OPEN (IFIL2,STATUS='SCRATCH',FORM='UNFORMATTED')
REWIND IFILD
REWIND IFIL1
REWIND IFIL2
C
WRITE (IFILD) V2
REWIND IFILD
CALL TRIALV (V1,V2)
C
WRITE (IOUT,'(/'' It Tr CPU'',2X,''Convergence'',6X,
>''Energy''/)')
DO 190 IT=1,MAXITS
WRITE (IFIL1) V1
C... CALL THE CI PROGRAM V2 = H . V1 + V2
DO 30 I=1,N
30 V2(I)=0.0D0
CALL FSIGMA (V1,V2)
WRITE (IFIL2) V2
C
REWIND IFIL1
DO 50 JT=1,IT
READ (IFIL1) V1
ZM(JT,IT) = 0.0D0
DO 40 I=1,N
40 ZM(JT,IT) = ZM(JT,IT) + V1(I)*V2(I)
50 ZM(IT,JT) = ZM(JT,IT)
REWIND IFIL1
REWIND IFIL2
C
C... SMALL MATRIX DIAGONALISER
DO 60 JT=1,IT
DO 60 KT=1,IT
60 ZV(KT,JT) = ZM(KT,JT)
CALL DIAG2 (MAXITR,IT,R,ZV)
C
TEST = 0
DO 61 IROOT=1,MIN(IT,NROOT)
TES = ABS(ZV(IT,IROOT))
IF (TES.LT.TEST) GOTO 61
ITRACK = IROOT
TEST = TES
61 CONTINUE
IF (IT.EQ.1) EREF=R(1)
WRITE (IOUT,70) IT,ITRACK,SECOND(),TEST
>,(R(IROOT)+CORE,IROOT=1,MIN(NROOT,IT))
70 FORMAT(I3,I4,F7.1,F12.8,10F20.9)
IF (TEST.LT.THR) GOTO 210
C
C... MAKE RESIDUAL
DO 80 I=1,N
80 V1(I) = 0.0D0
DO 100 JT=1,IT
READ (IFIL1) V2
DO 90 I=1,N
90 V1(I) = V1(I) + (-R(ITRACK)*ZV(JT,ITRACK))*V2(I)
READ (IFIL2) V2
DO 100 I=1,N
100 V1(I) = V1(I) + ZV(JT,ITRACK) * V2(I)
READ (IFILD) V2
IF (IT.EQ.1) THEN
DO 110 JJ = 1,N
IF (ABS(V2(JJ)-EREF).LT.1D-9) V2(JJ)=1D20
110 CONTINUE
END IF
REWIND IFILD
DO 130 I=1,N
130 V1(I) = V1(I) / (V2(I)-R(1)+1.0D-10)
C
C... ORTHOGONALISE
140 REWIND IFIL1
DO 160 JT=1,IT
READ (IFIL1) V2
ZZ = 0.0D0
DO 150 I=1,N
150 ZZ = ZZ + V1(I)*V2(I)
DO 160 I=1,N
160 V1(I) = V1(I) + (-ZZ) * V2(I)
ZZ = 0.0D0
DO 170 I=1,N
170 ZZ = ZZ + V1(I)**2
IF (ZZ.LT.1D-20) GOTO 210
ZZ = 1.0 / SQRT(ZZ)
DO 180 I=1,N
180 V1(I) = V1(I) * ZZ
C... REORTHOGONALISE IF A LOT HAS BEEN ANNIHILATED
IF (ZZ.GT.1D3) GOTO 140
C
190 CONTINUE
IT = MAXITS
WRITE (IOUT,200)
200 FORMAT(/1X,'*** Convergence not achieved in max iterations')
C
210 DO 280 IROOT=1,NROOT
REWIND IFIL1
DO 220 I=1,N
220 V1(I) = 0.0D0
DO 230 JT=1,IT
READ (IFIL1) V2
DO 230 I=1,N
230 V1(I) = V1(I) + ZV(JT,IROOT) * V2(I)
WRITE (IOUT,'('' State'',I2,5X,''Energy'',F23.12)')
> IROOT,R(IROOT)+CORE
IF (IROOT.EQ.1) WRITE (IOUT,'('' Correlation energy'',F23.12)')
> R(1)-EREF
WRITE (IOUT,'(/'' Final CI vector''/)')
CALL FCVECO (V1,ISYMH,THRRES)
280 CONTINUE
C
CLOSE (IFILD,STATUS='DELETE')
CLOSE (IFIL1,STATUS='DELETE')
CLOSE (IFIL2,STATUS='DELETE')
C
RETURN
END
SUBROUTINE TRIALV (V1,V2)
CSUBR SUBROUTINE TO GENERATE TRIAL CI VECTOR
CSUBR COULD BE REPLACED BY E.G. READING VECTOR IF REQUIRED
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /TAPES / INP,IOUT
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
COMMON /CHAMIL/ ISYMH,N
DIMENSION IMIN(64)
DIMENSION V1(N),V2(N)
TRIAL = 0
IF (NINT(TRIAL).EQ.0) THEN
C... FIND THE DETERMINANT WITH LOWEST ENERGY
VMIN = 1D20
DO 10 I=1,N
IF (V2(I).GT.VMIN) GOTO 10
IMIN(1) = I
VMIN = V2(I)
10 CONTINUE
C... CONSTRUCT A SPIN EIGENFUNCTION CONTAINING THIS DETERMINANT
CALL FCSPAD (V1,IMIN(1))
WRITE (IOUT,20)
20 FORMAT(/1X,'Initial configuration generated:')
NREF = 0
DO 40 I=1,N
IF (V1(I).EQ.0.0D0) GOTO 40
WRITE (IOUT,30) I,V1(I),V2(I)+CORE
NREF = NREF+1
IMIN(NREF) = I
30 FORMAT(I8,2F18.7)
40 CONTINUE
ELSE
C... READ FROM FILE
CONTINUE
END IF
RETURN
END
SUBROUTINE FCVECO (V,ISYM,THRESH)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /TAPES / INP,IOUT
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
CHARACTER*130 TEST
DIMENSION V(*)
COMMON /CLOCAL/ ICGA(32),ICGB(32)
N = NCI(ISYM)
IBASE = ICORR(0)
INTERA = ICORI(NA*NACT)
CALL INITI (IQ(INTERA),NA,NACT)
INTERB = ICORI(NB*NACT)
CALL INITI (IQ(INTERB),NA,NACT)
IC = ICORI(NSTRAA)
ITZ = ICORI(NSTRBB)
CALL FCMIC (IQ(IC),IQ(ITZ),ISYM)
NFROZ=0
DO 60 II=1,NCI(ISYM)
IF (ABS(V(II)).LT.THRESH) GOTO 60
CALL FCSTRG (II,ISYM,ICGA,ICGB)
WRITE (IOUT,'(1X,F20.12,64I3)') V(II)
1,(ICGA(I)+NFROZ,I=1,NA)
2,(ICGB(I)+NFROZ,I=1,NB)
60 CONTINUE
WRITE (IOUT,'(1X,''/EOF'')')
CALL CORLSR(IBASE)
RETURN
END
SUBROUTINE INITC
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
COMMON /TAPES / INP,IOUT
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
DIMENSION ICG(32)
C... LENGTHS OF CI EXPANSIONS
C... NSTRA, NSTRB ARE NUMBERS OF ALPHA,BETA STRINGS OF EACH SYMMETRY
C... NCI ARE NUMBERS OF DETERMINANTS OF EACH SYMMETRY
C... NPAIR ARE NUMBERS OF ORBITAL PAIRS OF EACH SYMMETRY
C.... INTOFF OFFSET FOR SYMMETRY BLOCKS OF TWO ELECTRON INTEGRALS
NSTRAA = 0
NSTRBB = 0
MAXAA = 0
MAXBB = 0
DO 100 ISYM=1,8
ICG(1) = 0
NSTRA(ISYM) = -1
80 NSTRA(ISYM) = NSTRA(ISYM)+1
CALL FCSTRS (NA,ICG,ISYM,ITYPEA)
IF (ICG(1).NE.0) GOTO 80
NSTRAA = NSTRAA + NSTRA(ISYM)
MAXAA = MAX(MAXAA,NSTRA(ISYM))
ICG(1) = 0
NSTRB(ISYM) = -1
90 NSTRB(ISYM) = NSTRB(ISYM)+1
CALL FCSTRS (NB,ICG,ISYM,ITYPEA)
IF (ICG(1).NE.0) GOTO 90
NSTRBB = NSTRBB + NSTRB(ISYM)
MAXBB = MAX(MAXBB,NSTRB(ISYM))
100 CONTINUE
DO 110 ISYM=1,8
NCI(ISYM) = 0
DO 110 ISYMB=1,8
110 NCI(ISYM) = NCI(ISYM) + NSTRA(MULT(ISYMB,ISYM))*NSTRB(ISYMB)
C
C... ORBITAL PAIR ADDRESSING, ETC..
DO 120 I=1,256
120 ICF(I) = (I*(I-1))/2
INTADT = ICORI(NACT*NACT)
NINT2=0
MAXPAR = 0
DO 140 ISYM=1,8
INTOFF(ISYM) = NINT2
NPAIR(ISYM)=0
NPB = 0
DO 130 I=1,NACT
DO 130 J=1,I
IF (ISYM.EQ.MULT(ITYPEA(I),ITYPEA(J))) THEN
NPAIR(ISYM) = NPAIR(ISYM)+1
IQ(INTADT+(I-1)*NACT+(J-1)) = NPAIR(ISYM)
IQ(INTADT+(I-1)+NACT*(J-1)) = NPAIR(ISYM)
END IF
130 CONTINUE
MAXPAR = MAX(NPAIR(ISYM),MAXPAR)
140 NINT2 = NINT2 + NPAIR(ISYM)**2
C... UPPER BOUND ON NUMBER OF SINGLE REPLACEMENTS FROM ANY STRING
MAXREP = MIN(MAXPAR,MAX(NA*(NACT+1-NA),NB*(NACT+1-NB)))
C
WRITE (IOUT,141)
> NACT,(NOCC(I),I=1,8)
>,NA+NB,(NA-NB)*0.5,NPAIR,NSTRA,NSTRB,NCI
141 FORMAT(
1' Active orbitals:', T30,I3,'(',8I3,')'/
2' Active electrons:', T30,I3/
3' Spin quantum number:', T30,F5.1/
4' Orbital pairs:', T30,8I9/
5' Strings:', T30,8I9/T30,8I9/
6' Determinants:', T30,4I12/T30,4I12)
C WRITE (IOUT,'('' ACTIVE ORBITAL SYMMETRIES:'')')
C WRITE (IOUT,150) (ITYPEA(I),I=1,NACT)
C150 FORMAT(1X,40I2)
C
C... LOAD INTEGRALS
INTAA = ICORR(NINT2)
INTHAA = ICORR(NACT)
INTJAA = ICORR(NACT**2)
INTKAA = ICORR(NACT**2)
IBASE = ICORR(0)
IONEA = ICORR(NPAIR(1))
CALL INPINT (Q(IONEA),Q(INTAA))
C
C... INTEGRALS FOR DIAGONAL ELEMENTS OF HAMILTONIAN
DO 10 I=1,NACT
10 Q(INTHAA-1+I) = Q(IONEA-1+IQ(INTADT+(I-1)*(NACT+1)))
DO 11 I=1,NACT
DO 11 J=1,NACT
IJ = IQ(INTADT+(I-1)*NACT+J-1)
II = IQ(INTADT+(I-1)*NACT+I-1)
JJ = IQ(INTADT+(J-1)*NACT+J-1)
ISYMIJ = MULT(ITYPEA(I),ITYPEA(J))
Q(INTJAA+(I-1)*NACT+J-1) = Q(INTAA+(II-1)*NPAIR(1)+JJ-1)
11 Q(INTKAA+(I-1)*NACT+J-1) =
> Q(INTAA+INTOFF(ISYMIJ)+(IJ-1)*(NPAIR(ISYMIJ)+1))
C
C... MODIFY INTEGRALS SUCH THAT HAMILTONIAN BECOMES
C SUM(IJKL) E(IJ) E(KL) (IJ|KL)"
C (NO 1 ELECTRON TERM REMAINS)
C
C... (IJ|KL)' = 0.5 * (IJ|KL)
DO 20 IJKL=1,NINT2
20 Q(INTAA-1+IJKL) = Q(INTAA-1+IJKL)*0.5D0
C
C... H'(IJ) = H(IJ) - SUM(K) (IK|KJ)'
IJ = IONEA-1
DO 40 I=1,NACT
DO 40 J=1,I
IF (ITYPEA(I).NE.ITYPEA(J)) GOTO 40
IJ = IJ+1
DO 30 K=1,NACT
ISYMIK = MULT(ITYPEA(I),ITYPEA(K))
IK =IQ(INTADT+(I-1)*NACT+K-1)-1
JK =IQ(INTADT+(J-1)*NACT+K-1)-1
30 Q(IJ) = Q(IJ) - Q(INTAA+INTOFF(ISYMIK)+IK*NPAIR(ISYMIK)+JK)
40 CONTINUE
C
C... (IJ|KL)" = (IJ|KL)'+(1/2*NELEC)*(DEL(IJ)*H'(KL)+DEL(KL)*H'(IJ))
FAC = 2*(NA+NB)
FAC = 1D0/FAC
DO 70 K=1,NACT
KK = IQ(INTADT+(K-1)*(NACT+1))
DO 50 IJ=1,NPAIR(1)
50 Q(INTAA+(KK-1)*NPAIR(1)+IJ-1) = Q(INTAA+(KK-1)*NPAIR(1)+IJ-1)
> + FAC*Q(IONEA-1+IJ)
DO 60 IJ=1,NPAIR(1)
60 Q(INTAA+(IJ-1)*NPAIR(1)+KK-1) = Q(INTAA+(IJ-1)*NPAIR(1)+KK-1)
> + FAC*Q(IONEA-1+IJ)
70 CONTINUE
CALL CORLSR(IBASE)
WRITE (IOUT,'('' Core energy:'',T30,F20.12)') CORE
RETURN
END
SUBROUTINE FSIGMA (C,S)
CSUBR RETURNS S = H * C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
COMMON /TAPES / INP,IOUT
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
LOGICAL TIMING
DIMENSION IOFFCI(8),C(*),S(*)
COMMON /CHAMIL/ ISYMH,NNCI
TIMING = IPRTIM.GE.0
NOPMXM = 0
NOPMX = 0
NOPOTH = 0
NOPCCA = 0
NOPCCB = 0
NLOCCA = 0
NLOCCB = 0
NOPONE = 0
T0 = SECOND()
TMXM = 0
TCCA = 0
TCCB = 0
TONE = 0
IBASE = ICORR(0)
NNCI = 0
DO 10 ISYMB=1,8
ISYMA = MULT(ISYMH,ISYMB)
IOFFCI(ISYMB) = NNCI
IF (IPRHAM.GT.0) WRITE (IOUT,*)'ISYMB,ISYMA,NNCI,NSTRA,NSTRB '
1,ISYMB,ISYMA,NNCI,NSTRA(ISYMA),NSTRB(ISYMB)
10 NNCI = NNCI + NSTRA(ISYMA)*NSTRB(ISYMB)
IF (IPRHAM.GT.0) CALL OUTVEC (C,NNCI,'c passed to sigma')
IC = ICORI(NSTRAA)
ITZ = ICORI(NSTRBB)
CALL FCMIC (IQ(IC),IQ(ITZ),ISYMH)
DO 170 ISYMA=1,8
DO 170 ISYMB=1,8
IBASES = ICORR(0)
ISYMR = MULT(ISYMH,MULT(ISYMA,ISYMB))
ISYMBE = MULT(ISYMB,ISYMR)
NSBE = NSTRB(ISYMBE)
NSA = NSTRA(ISYMA)
NSB = NSTRB(ISYMB)
NPAIRR = NPAIR(ISYMR)
IF (NSA.EQ.0.OR.NSB.EQ.0.OR.NPAIRR.EQ.0) GOTO 170
MAXREA = MIN(NA*(NACT-NA+1),NPAIRR)
MAXREB = MAX(1,MIN(NB*(NACT-NB+1),NPAIRR))
MAXR2A = MAXREA * 2
MAXR3A = MAXREA * 3
MAXR2B = MAXREB * 2
MAXR3B = MAXREB * 3
IBB = ICORR(0)
IWA = ICORI(MAXR3A*NSA)
NWA = ICORI(NSA)
IWB = ICORI(MAXR3B*NSB)
NWB = ICORI(NSB)
C... CONSTRUCT ALL SINGLE REPLACEMENTS
T1 = SECOND()
INTER = ICORI(NA*NACT)
CALL INITI (IQ(INTER),NA,NACT)
CALL ONELAL (NA,ITYPEA,IQ(INTER),IQ(INTADT),ISYMA,ISYMR
1 ,IQ(IWA),IQ(NWA),IQ(IC),MAXREA)
CALL CORLSI(INTER)
INTER = ICORI(NB*NACT)
CALL INITI (IQ(INTER),NB,NACT)
CALL ONELAL (NB,ITYPEA,IQ(INTER),IQ(INTADT),ISYMB,ISYMR
1 ,IQ(IWB),IQ(NWB),IQ(ITZ),MAXREB)
TONE = TONE+SECOND()-T1
CALL CORLSI(INTER)
C... DETERMINE BLOCKING
NMAT = 2
MWORD = ICORRM()/NMAT
WORDS = MWORD
NWORD = SQRT(WORDS/NPAIRR)
NABLK = (MIN(NWORD,NSA)/2)*2+1
DO 20 NBBLK=(NSB/2)*2+1,1,-2
IF (NABLK*NBBLK*NPAIRR.LE.MWORD) GOTO 30
20 CONTINUE
30 IDA = ICORR(NABLK*NBBLK*NPAIRR)
IEA = ICORR(NABLK*NBBLK*NPAIRR)
IDB = IDA
IEB = IEA
C... LOOP OVER BLOCKS
DO 160 IOFFA=0,NSA-1,NABLK
NAA = MIN(NABLK,NSA-IOFFA)
DO 160 IOFFB=0,NSB-1,NBBLK
NBB = MIN(NBBLK,NSB-IOFFB)
NAABB = NAA*NBB
IF (TIMING) NOPOTH = NOPOTH+NAABB*NPAIRR
CALL FZERO (Q(IDA),NAABB*NPAIRR)
C... ALPHA REPLACEMENTS
IF (TIMING) THEN
T1 = SECOND()
ICNT = 0
END IF
IATABB = IWA+IOFFA*MAXR3A
IDAA = IDA-1-NAABB
DO 61 IA=1,NAA
IATAB = IATABB
IF (TIMING) ICNT = ICNT+IQ(NWA+IOFFA-1+IA)
DO 60 IWW=1,IQ(NWA+IOFFA-1+IA)
ICCC = IOFFB + IQ(IATAB)
IAAA = IDAA+IQ(IATAB+MAXR2A)*NAABB
IF (IQ(IATAB+MAXREA).LT.0) THEN
C... EQUATION (18)
DO 40 IB=1,NBB
40 Q(IAAA+IB) = - C(ICCC+IB)
ELSE
C... EQUATION (18)
DO 50 IB=1,NBB
50 Q(IAAA+IB) = C(ICCC+IB)
END IF
60 IATAB = IATAB+1
IDAA = IDAA + NBB
61 IATABB = IATABB + MAXR3A
IF (TIMING) THEN
TCCA = TCCA + SECOND()-T1
NOPCCA = NOPCCA + ICNT*NBB
NLOCCA = NLOCCA + ICNT
END IF
C... BETA REPLACEMENTS
IF (TIMING) THEN
ICNT = 0
T1 = SECOND()
END IF
IBTABB = IWB+IOFFB*MAXR3B
ICCCC = IOFFCI(ISYMBE)+(IOFFA-1)*NSBE
IDBB = IDB-NBB-NAABB
DO 91 IB=1,NBB
IBTAB = IBTABB
IF (TIMING) ICNT = ICNT+IQ(NWB+IOFFB-1+IB)
DO 90 IWW=1,IQ(NWB-1+IOFFB+IB)
ICCC = ICCCC+IQ(IBTAB)
IBBB = IDBB+IQ(IBTAB+MAXR2B)*NAABB
IF (IQ(IBTAB+MAXREB).LT.0) THEN
C... EQUATION (19)
DO 70 IA=1,NAA
70 Q(IBBB+IA*NBB) = Q(IBBB+IA*NBB) - C(ICCC+IA*NSBE)
ELSE
C... EQUATION (19)
DO 80 IA=1,NAA
80 Q(IBBB+IA*NBB) = Q(IBBB+IA*NBB) + C(ICCC+IA*NSBE)
END IF
90 IBTAB = IBTAB + 1
IDBB = IDBB + 1
91 IBTABB = IBTABB + MAXR3B
IF (TIMING) THEN
TCCB = TCCB + SECOND()-T1
NOPCCB = NOPCCB + ICNT*NAA
NLOCCB = NLOCCB + ICNT
END IF
C
C... MATRIX MULTIPLICATION WITH INTEGRALS
C... EQUATION (20)
IF (TIMING) T1 = SECOND()
IF (IPRHAM.GT.3)
1CALL OUTSQR (Q(IDA),NAABB,NAABB,NPAIRR,'DA')
IF (IPRHAM.GT.3)
1CALL OUTSQR (Q(INTAA+INTOFF(ISYMR)),NPAIRR,NPAIRR,NPAIRR,'INTAA')
CALL FCMXM (Q(IDA),Q(INTAA+INTOFF(ISYMR)),Q(IEA),NAABB,NPAIRR,
>NOPMXM)
NOPMX = NOPMX+2*NAABB*NPAIRR**2
IF (IPRHAM.GT.3)
1CALL OUTSQR (Q(IEA),NAABB,NAABB,NPAIRR,'EA')
IF (TIMING) TMXM = TMXM+SECOND()-T1
C
C... ALPHA REPLACEMENTS
IF (TIMING) THEN
T1 = SECOND()
ICNT = 0
END IF
IATABB = IWA+IOFFA*MAXR3A
IEAA = IEA-1-NAABB
DO 121 IA=1,NAA
IATAB = IATABB
IF (TIMING) ICNT = ICNT+IQ(NWA+IOFFA-1+IA)
DO 120 IWW=1,IQ(NWA+IOFFA-1+IA)
ICCC = IOFFB + IQ(IATAB)
IAAA = IEAA+IQ(IATAB+MAXR2A)*NAABB
IF (IQ(IATAB+MAXREA).LT.0) THEN
C... EQUATION (21)
DO 100 IB=1,NBB
100 S(ICCC+IB) = S(ICCC+IB) - Q(IAAA+IB)
ELSE
C... EQUATION (21)
DO 110 IB=1,NBB
110 S(ICCC+IB) = S(ICCC+IB) + Q(IAAA+IB)
END IF
120 IATAB = IATAB + 1
IEAA = IEAA + NBB
121 IATABB = IATABB + MAXR3A
IF (TIMING) THEN
TCCA = TCCA + SECOND()-T1
NOPCCA = NOPCCA + ICNT*NBB
NLOCCA = NLOCCA + ICNT
END IF
C... BETA REPLACEMENTS
IF (TIMING) THEN
T1 = SECOND()
ICNT = 0
END IF
IBTABB = IWB+IOFFB*MAXR3B
ICCCC = IOFFCI(ISYMBE)+(IOFFA-1)*NSBE
IEBB = IEB-NBB-NAABB
DO 151 IB=1,NBB
IBTAB = IBTABB
IF (TIMING) ICNT = ICNT+IQ(NWB+IOFFB-1+IB)
DO 150 IWW=1,IQ(NWB-1+IOFFB+IB)
ICCC = ICCCC+IQ(IBTAB)
IBBB = IEBB+IQ(IBTAB+MAXR2B)*NAABB
IF (IQ(IBTAB+MAXREB).LT.0) THEN
C... EQUATION (22)
DO 130 IA=1,NAA
130 S(ICCC+IA*NSBE) = S(ICCC+IA*NSBE) - Q(IBBB+IA*NBB)
ELSE
C... EQUATION (22)
DO 140 IA=1,NAA
140 S(ICCC+IA*NSBE) = S(ICCC+IA*NSBE) + Q(IBBB+IA*NBB)
END IF
150 IBTAB = IBTAB + 1
IEBB = IEBB + 1
151 IBTABB = IBTABB + MAXR3B
IF (TIMING) THEN
TCCB = TCCB + SECOND()-T1
NOPCCB = NOPCCB + ICNT*NAA
NLOCCB = NLOCCB + ICNT
END IF
160 CONTINUE
CALL CORLSR(IBASES)
170 CONTINUE
CALL CORLSR(IBASE)
IF (TIMING) THEN
TTOT = SECOND()-T0
TOTH = TTOT-TMXM-TCCA-TCCB-TONE
NOPTOT = NOPOTH+NOPMXM+NOPCCA+NOPCCB+NOPONE
CALL TPRINT ('String c. coeff',NOPONE,TONE,' ',0)
CALL TPRINT ('Coupl coeffs(A)',NOPCCA,TCCA,'Av. v length'
>,NOPCCA/NLOCCA)
CALL TPRINT ('Coupl coeffs(B)',NOPCCB,TCCB,'Av. v length'
>,NOPCCB/NLOCCB)
IAVSP = NINT(1D2-(1D2*NOPMXM)/NOPMX)
CALL TPRINT ('Matrix multiply',NOPMXM,TMXM,'Av. % sparsity',IAVSP)
CALL TPRINT ('Other ',NOPOTH,TOTH,' ',0)
CALL TPRINT ('Total ',NOPTOT,TTOT,' ',0)
END IF
IF (IPRHAM.GT.0)
1CALL OUTVEC (S,NNCI,'S RETURNED BY SIGMA')
RETURN
END
SUBROUTINE ONELAL (N,ITYPE,INTER,INTAD,ISYM,ISYMR,IW,NW,ITZ,MAXRE)
CSUBR CONSTRUCT AND STORE ALL ONE ELECTRON REPLACEMENTS
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C... MAXIMUM NUMBERS OF ORBITALS IN ANY SYMMETRY MUST BE LESS THANTHIS!
PARAMETER (MXORB=64)
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
DIMENSION ITYPE(*),INTER(N,*),INTAD(*)
1 ,IW(*),NW(*),ITZ(*)
C 1 ,IW(MAXRE,3,*),NW(*),ITZ(*)
DIMENSION ICG(32),ICGI(33),JSTOR1(MXORB),JSTOR2(MXORB)
>,JSTOR3(MXORB)
nw(1) = 0
if (n.le.0) return
MAXRE2 = MAXRE*2
MAXRE3 = MAXRE*3
C... SET COUNTERS TO ZERO AND SET UP LOOK-UP FOR BRAS
NINSYM = IBINOM(NACT,N)
ITZI = ICORI(NINSYM)
NI = 0
ICGI(1) = 0
DO 23 II=1,2999999
IISYM = 0
CALL FCSTRS (N,ICGI,IISYM,ITYPE)
IF (ICGI(1).EQ.0) GOTO 24
IF (IISYM.NE.ISYM) GOTO 23
NI=NI+1
IQ(ITZI-1+II) = NI
23 CONTINUE
24 CONTINUE
DO 25 I=1,NI
25 NW(I) = 0
IF (NI.EQ.0) GOTO 999
C... LOOP OVER INTERMEDIATE N-1 ELECTRON STRINGS K
ICG(1) = 0
DO 50 ISK=1,2999999
ISYMK = 0
CALL FCSTRS (N-1,ICG,ISYMK,ITYPE)
IF (ICG(1).EQ.0) GOTO 999
ICG(N) = NACT+1
ISYMIO = MULT(ISYMK,ISYM)
ISYMJO = MULT(ISYMR,ISYMIO)
ISYMJ = MULT(ISYM,ISYMR)
C... CONSTRUCT CREATIONS TO KETS J (SYMMETRY ISYMJ)
C DO 777 IREP=1,10
IREPL = 1
IADR = 1
IPARIT = 1
DO 10 J=1,N-1
10 IADR = IADR + INTER(J+1,ICG(J))
NJ = 0
DO 30 JO=1,NACT
IF(JO.EQ.ICG(IREPL))THEN
IADR = IADR+INTER(IREPL,JO)-INTER(IREPL+1,JO)
IREPL = IREPL+1
IPARIT = -IPARIT
ELSE IF (ITYPE(JO).EQ.ISYMJO) THEN
C... GET CONNECTED N ELECTRON DETERMINANT J
C... STORE INFO
NJ = NJ+1
JSTOR1(NJ) = ITZ(IADR+INTER(IREPL,JO))
JSTOR2(NJ) = IPARIT
JSTOR3(NJ) = JO
END IF
30 CONTINUE
C777 CONTINUE
IF (NJ.EQ.0) GOTO 50
C... CONSTRUCT CREATIONS TO BRAS I (SYMMETRY ISYM)
IREPL = 1
IADR = ITZI
IPARIT = 1
DO 110 J=1,N-1
110 IADR = IADR + INTER(J+1,ICG(J))
DO 130 IO=1,NACT
IF(IO.EQ.ICG(IREPL))THEN
IADR = IADR+INTER(IREPL,IO)-INTER(IREPL+1,IO)
IREPL = IREPL+1
IPARIT = -IPARIT
ELSE IF (ITYPE(IO).EQ.ISYMIO) THEN
C... GET CONNECTED N ELECTRON DETERMINANT I
I = IQ(IADR+INTER(IREPL,IO))
IOFFS = (I-1)*MAXRE3+NW(I)
IOO = (IO-1)*NACT
C... STORE INFO FOR EACH J
IF (IPARIT.GT.0) THEN
CDIR$ SHORTLOOP
CDIR$ IVDEP
C$DIR NO_RECURRENCE
*VOPTION INDEP VEC
DO 45 J=1,NJ
IW(IOFFS +J) = JSTOR1(J)
IW(IOFFS+MAXRE2+J) = JSTOR3(J)+IOO
C IW(IOFFS+MAXRE2+J) = INTAD(JSTOR3(J)+IOO)
45 IW(IOFFS+MAXRE +J) = JSTOR2(J)
ELSE
CDIR$ SHORTLOOP
CDIR$ IVDEP
C$DIR NO_RECURRENCE
*VOPTION INDEP VEC
DO 46 J=1,NJ
IW(IOFFS +J) = JSTOR1(J)
IW(IOFFS+MAXRE2+J) = JSTOR3(J)+IOO
C IW(IOFFS+MAXRE2+J) = INTAD(JSTOR3(J)+IOO)
46 IW(IOFFS+MAXRE +J) = -JSTOR2(J)
END IF
NW(I) = NW(I) + NJ
END IF
130 CONTINUE
50 CONTINUE
999 CALL CORLSI(ITZI)
C... GATHER LOOP FOR PACKED INTEGRAL ADDRESSES
IOFFS = MAXRE2
DO 147 I=1,NI
DO 148 J=1,NW(I)
148 IW(IOFFS+J) = INTAD(IW(IOFFS+J))
147 IOFFS = IOFFS + MAXRE3
RETURN
END
SUBROUTINE FCDIAG (S)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /TAPES / INP,IOUT
COMMON /CHAMIL/ ISYMH,NNCI
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
DIMENSION S(*)
COMMON /CLOCAL/ ICGA(32),ICGB(32)
IF (IPRTIM.GE.0) T0 = SECOND()
NOP = 0
IF (IPRDIA.GT.0) THEN
C WRITE (IOUT,*) 'COMPUTE DIAGONAL ELEMENTS'
CALL OUTVEC (Q(INTHAA),NACT,'HAA')
CALL OUTSQR (Q(INTJAA),NACT,NACT,NACT,'INTJAA')
CALL OUTSQR (Q(INTKAA),NACT,NACT,NACT,'INTKAA')
END IF
C IF (IPRDIA.GT.1) WRITE (IOUT,*) 'MAXAA,MAXBB ',MAXAA,MAXBB
IZB = ICORR(NACT*MAXBB)
IZA = ICORR(NACT*MAXAA)
IZ = ICORR(NACT*MAXBB)
IF = ICORR(MAXBB)
ICCC=0
DO 220 ISYMB=1,8
ISYMA=MULT(ISYMB,ISYMH)
NSA = NSTRA(ISYMA)
NSB = NSTRB(ISYMB)
IF (IPRDIA.GT.1) WRITE (IOUT,*) 'ISYMB,ISYMA,NSA,NSB ',ISYMB
1,ISYMA,NSA,NSB
IF (NSA*NSB.LE.0) GOTO 220
CALL FZERO (Q(IZB),NACT*NSB)
C.. STORE BETA STRING OCCUPANCIES
ICGB(1) = 0
DO 20 IB=1,NSB
CALL FCSTRS (NB,ICGB,ISYMB,ITYPEA)
DO 10 II=1,NB
10 Q(IZB-1+IB+(ICGB(II)-1)*NSB)=1
20 CONTINUE
IF (IPRDIA.GT.2) CALL OUTSQR (Q(IZB),NSB,NSB,NACT,'ZB')
C
C.. STORE ALPHA STRING OCCUPANCIES
ICGA(1) = 0
CALL FZERO (Q(IZA),NACT*NSA)
DO 40 IA=1,NSA
CALL FCSTRS (NA,ICGA,ISYMA,ITYPEA)
DO 40 II=1,NA
40 Q(IZA-1+IA+(ICGA(II)-1)*NSA) = 1
IF (IPRDIA.GT.2) CALL OUTSQR (Q(IZA),NSA,NSA,NACT,'ZA')
C
DO 50 IAB=1,NSA*NSB
50 S(ICCC+IAB) = 0
C
C... SPECIAL CODE FOR RHF ORBITALS
DO 210 IA=1,NSA
C... OCCUPATION NUMBERS OF DETERMINANTS
DO 100 I=1,NACT
DO 100 IB=1,NSB
100 Q(IZ-1+IB+(I-1)*NSB) = Q(IZB-1+IB+(I-1)*NSB)
1 + Q(IZA-1+IA+(I-1)*NSA)
IF (IPRDIA.GT.3) CALL OUTSQR (Q(IZ),NSB,NSB,NACT,'Z')
C
C... CONSTANT EXCHANGE CONTRIBUTION & COULOMB
C... 1/2 SUM(IJ) N(I) N(J) ((II|JJ)-0.5(IJ|JI))
C... EQUATIONS (23), (25)
IJ = -1
DO 130 I=1,NACT
ZZ = Q(INTHAA-1+I)
IF (IPRDIA.GT.3) WRITE (IOUT,*) '1567;I,HAA(I) ',I,ZZ
DO 110 IB=1,NSB
110 S(ICCC+IB) = S(ICCC+IB) + ZZ*Q(IZ-1+IB+(I-1)*NSB)
DO 130 J=1,I
IJ = (I-1)*NACT-1+J
ZZ = Q(INTJAA+IJ)-0.5D0*Q(INTKAA+IJ)
IF (I.EQ.J) ZZ = ZZ*0.5D0
IF (IPRDIA.GT.3) WRITE (IOUT,*) '1567; I,J,ZZ ',I,J,ZZ
DO 120 IB=1,NSB
120 S(ICCC+IB) = S(ICCC+IB)
1 + ZZ*Q(IZ-1+IB+(I-1)*NSB)*Q(IZ-1+IB+(J-1)*NSB)
130 CONTINUE
C
IF (IPRDIA.GT.3) CALL OUTVEC (S(ICCC+1),NSB,'S AFTER CONSTS')
C... VV = 4*S**2, SEE EQUATION (26)
VV = (NA-NB)**2
C.. Z TO HOLD 1 IF SINGLY OCC, 0 OTHERWISE
C... (CAPITAL N = N(2-N), EQUATION (25))
CALL FZERO (Q(IF),NSB)
DO 150 I=1,NACT
DO 140 IB=1,NSB
140 Q(IZ-1+IB+(I-1)*NSB) = Q(IZ-1+IB+(I-1)*NSB)
1 * (2.0D0-Q(IZ-1+IB+(I-1)*NSB))
DO 150 IB=1,NSB
150 Q(IF-1+IB) = Q(IF-1+IB) + Q(IZ-1+IB+(I-1)*NSB)
DO 160 IB=1,NSB
160 Q(IF-1+IB) = MAX(Q(IF-1+IB),1.1D0)
DO 170 IB=1,NSB
170 Q(IF-1+IB) = (VV-Q(IF-1+IB))/( Q(IF-1+IB) * (Q(IF-1+IB)-1.0D0) )
C.. F NOW HOLDS F FACTOR, EQUATION (26)
C
IF (IPRDIA.GT.3) CALL OUTVEC (Q(IF),NSB,'F')
IJ = -1
DO 200 I=1,NACT
IF (I.EQ.1) GOTO 190
DO 180 J=1,I-1
IJ = (I-1)*NACT-1+J
ZZ = -0.5D0*Q(INTKAA+IJ)
DO 180 IB=1,NSB
180 S(ICCC+IB) = S(ICCC+IB)
1 + ZZ*Q(IF-1+IB)*Q(IZ-1+IB+(I-1)*NSB)*Q(IZ-1+IB+(J-1)*NSB)
190 IJ = (I-1)*NACT-1+I
ZZ = -0.25D0*Q(INTKAA+IJ)
IF (IPRDIA.GT.3) WRITE (IOUT,*) '1875; I,ZZ ',I,ZZ
DO 200 IB=1,NSB
200 S(ICCC+IB) = S(ICCC+IB) + ZZ*Q(IZ-1+IB+(I-1)*NSB)
210 ICCC = ICCC+NSB
NOP = NOP+NSA*NSB*(NACT*7+3*NACT*(NACT+1)+5)
220 CONTINUE
CALL CORLSR(IZB)
IF (IPRDIA.GE.0) CALL OUTVEC (S,ICCC,'DIAGONAL ELEMENTS')
IF (IPRTIM.GE.0) CALL TPRINT ('Hamiltonian diag. els.'
>,NOP,SECOND()-T0,' ',0)
RETURN
END
FUNCTION INICOR(MEMREQ)
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INTEGER CORMON,CORLSR,CORLSI,ENDCOR
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
COMMON /TAPES / INP,IOUT
COMMON /CORCTL/ INTREL,LWORD,LTOP,LMAX,LMIN,MREAL,MXBFF,IBASE,MEM
SAVE ICOUNT
DATA ICOUNT/0/
C
C... EXECUTE ENDCOR CODE IF NOT FIRST CALL TO THIS ROUTINE
IF (ICOUNT.NE.0) CALL FMAIN (Q(IBASE+1),MEM)
ICOUNT = ICOUNT+1
C
MEM = MEMREQ
LENN = MEM+2
C
C
CALL GMAINV (Q,IBASE,LENN)
IBASE = IBASE-1
IF (LENN.LT.MEM) THEN
WRITE (IOUT,'('' Request for'',I9,'' words refused'',
> '', '',I9,'' available'')') MEM,LENN
STOP 'Insufficient memory'
END IF
WRITE (IOUT,'('' Variable memory set to '',I10,'' words'')') MEM
LWORD = MEM+IBASE
LTOP = IBASE
LMAX = IBASE
LMIN = IBASE
INICOR = LWORD-IBASE
RETURN
C
ENTRY ICORIM ()
ICORIM = (LWORD - LTOP) * INTREL
RETURN
C
ENTRY ICORRM ()
ICORRM = LWORD - LTOP
RETURN
C
ENTRY ICORI (NWOR)
NWORD = IABS(NWOR)
NW = (NWORD+INTREL-1)/INTREL
MFR = LWORD - LTOP
IF (NW .GT. MFR) GOTO 20
ICORI = LTOP*INTREL+1
LTOP=LTOP+NW
IF (NWOR.LT.0) LMAX = MAX0(LMAX,LTOP)
IF (NWOR.GT.0) LMIN = MAX0(LMIN,LTOP)
RETURN
C
ENTRY ICORR (NWOR)
NWORD = IABS(NWOR)
NW = NWORD
MFR = LWORD - LTOP
IF (NW .GT. MFR) GOTO 20
ICORR = LTOP+1
LTOP=LTOP+NW
IF (NWOR.LT.0) LMAX = MAX0(LTOP,LMAX)
IF (NWOR.GT.0) LMIN = MAX0(LTOP,LMIN)
RETURN
C
20 CONTINUE
WRITE (IOUT,*) 'insufficient memory available - require ',NW,
> ' have ',MFR
IF (NW .NE. NWORD) THEN
WRITE(IOUT,*) 'the request was for ',NWORD,' integer words'
ELSE
WRITE(IOUT,*) 'the request was for real words'
END IF
STOP 'Insufficient memory'
END
SUBROUTINE CORLSI (IAD)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
COMMON /TAPES / INP,IOUT
COMMON /CORCTL/ INTREL,LWORD,LTOP,LMAX,LMIN,MREAL,MXBFF,IBASE,MEM
LTOP = (IAD-2+INTREL)/INTREL
RETURN
C
ENTRY CORLSR (IAD)
LTOP = IAD - 1
RETURN
C
ENTRY CORMON()
WRITE (IOUT,40) LMIN-IBASE,LWORD-LMIN
40 FORMAT(/1X,'Minimal core high water =',I8,
1' real numbers (spare core =',I8,')')
RETURN
C
ENTRY ENDCOR()
CALL FMAIN (Q(IBASE+1),MEM)
WRITE (IOUT,'('' Variable memory released'')')
RETURN
END
FUNCTION FCISTR (IPAR,INTER,ICG,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION ICG(N),INTER(N,*)
DIMENSION ICGG(32)
DO 10 J=1,N
10 ICGG(J) = ICG(J)
IF (N.LT.2) GOTO 40
DO 30 JJ=2,N
JJM = JJ-1
I = JJM
DO 20 II=1,JJM
IF (ICGG(I).LT.ICGG(I+1)) GOTO 20
J1=ICGG(I)
ICGG(I)=ICGG(I+1)
ICGG(I+1)=J1
IPAR=-IPAR
20 I=I-1
30 CONTINUE
40 CONTINUE
FCISTR=1
DO 50 J=1,N
50 FCISTR = FCISTR+INTER(J,ICGG(J))
RETURN
END
SUBROUTINE FCMIC (IC,ITZ,ISYM)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
DIMENSION ICG(32)
DIMENSION IC(*),ITZ(*)
DO 1 I=1,8
NSTRA(I) = 0
1 NSTRB(I) = 0
NSTRAT = IBINOM(NACT,NA)
NSTRBT = IBINOM(NACT,NB)
C... ALPHA STRING INFORMATION
ICG(1) = 0
DO 11 MT=1,NSTRAT
JMT = 0
CALL FCSTRS (NA,ICG,JMT,ITYPEA)
NSTRA(JMT) = NSTRA(JMT)+1
C... STORE SYMMETRY FOR LATER PROCESSING
11 IC(MT) = -JMT
C... BETA STRING INFORMATION
ICG(1) = 0
DO 10 MT=1,NSTRBT
JMT = 0
CALL FCSTRS (NB,ICG,JMT,ITYPEA)
NSTRB(JMT) = NSTRB(JMT)+1
10 ITZ(MT) = NSTRB(JMT)
C
C... LOOP OVER SYMMETRIES BUILDING TOTAL OFFSETS
NNCI = 0
DO 110 ISYMB=1,8
ISYMA = MULT(ISYM,ISYMB)
DO 120 MT=1,NSTRAT
IF (IC(MT).NE.-ISYMA) GOTO 120
IC(MT) = NNCI
NNCI = NNCI + NSTRB(ISYMB)
120 CONTINUE
110 CONTINUE
RETURN
END
SUBROUTINE FCSPAD (C,IREF)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
COMMON /BIG/ Q(1)
INTEGER IQ(1)
EQUIVALENCE (Q(1),IQ(1))
DIMENSION C(*)
COMMON /CHAMIL/ ISYMH,NNCI
DIMENSION ICGA(32),ICGB(32),IDOCA(32)
1 ,IDOCB(32),IALP(32),IBET(32)
DO 10 I=1,NNCI
10 C(I)=0
IBASE = ICORR(0)
INTERA = ICORI(NA*NACT)
INTERB = ICORI(NB*NACT)
IC = ICORI(NSTRAA)
ITZ = ICORI(NSTRBB)
CALL INITI (IQ(INTERA),NA,NACT)
CALL INITI (IQ(INTERB),NB,NACT)
CALL FCMIC (IQ(IC),IQ(ITZ),ISYMH)
CALL FCSTRG (IREF,ISYMH,ICGA,ICGB)
NDOC=0
NALP=0
NBET=0
DO 70 I=1,NACT
ICASE=1
DO 20 J=1,NA
IF (ICGA(J).NE.I) GOTO 20
ICASE=3
JA=J
20 CONTINUE
DO 30 J=1,NB
IF (ICGB(J).NE.I) GOTO 30
ICASE=ICASE+1
JB=J
30 CONTINUE
GOTO (70,40,50,60),ICASE
40 NBET=NBET+1
IBET(NBET)=JB
GOTO 70
50 NALP=NALP+1
IALP(NALP)=JA
GOTO 70
60 NDOC=NDOC+1
IDOCA(NDOC)=JA
IDOCB(NDOC)=JB
70 CONTINUE
FAC = 1.0D0/SQRT(DBLE(2**NBET))
IF (NBET.GT.5) STOP 5
DO 120 I5=1,2
DO 110 I4=1,2
DO 100 I3=1,2
DO 90 I2=1,2
DO 80 I1=1,2
IPAR = 1
MTA = FCISTR (IPAR,IQ(INTERA),ICGA,NA)
MTB = FCISTR (IPAR,IQ(INTERB),ICGB,NB)
C(IQ(IC-1+MTA)+IQ(ITZ-1+MTB)) = DBLE(IPAR)*FAC
IF (NBET.LT.1) GOTO 130
J1=ICGB(IBET(1))
ICGB(IBET(1))=ICGA(IALP(1))
ICGA(IALP(1))=J1
80 CONTINUE
IF (NBET.LT.2) GOTO 130
J1=ICGB(IBET(2))
ICGB(IBET(2))=ICGA(IALP(2))
ICGA(IALP(2))=J1
90 CONTINUE
IF(NBET.LT.3) GOTO 130
J1=ICGB(IBET(3))
ICGB(IBET(3))=ICGA(IALP(3))
ICGA(IALP(3))=J1
100 CONTINUE
IF (NBET.LT.4) GOTO 130
J1=ICGB(IBET(4))
ICGB(IBET(4))=ICGA(IALP(4))
ICGA(IALP(4))=J1
110 CONTINUE
IF (NBET.LT.5) GOTO 130
J1=ICGB(IBET(5))
ICGB(IBET(5))=ICGA(IALP(5))
ICGA(IALP(5))=J1
120 CONTINUE
130 CONTINUE
CALL CORLSR (IBASE)
RETURN
END
SUBROUTINE FCSTRG (II,ISYM,ICGA,ICGB)
CSUBR EXTRACT STRINGS ICGA,ICGB FOR CONFIGURATION II IN SYMMETRY ISYM
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
DIMENSION ICGA(*),ICGB(*)
III = II
DO 10 ISYMB=1,8
ISYMA = MULT(ISYMB,ISYM)
NSB = NSTRB(ISYMB)
NSA = NSTRA(ISYMA)
IF (III.LE.NSB*NSA) GOTO 20
10 III = III - NSB*NSA
20 ISA = (III-1)/NSB+1
ISB = III - (ISA-1)*NSB
ICGA(1) = 0
DO 30 JSA=1,ISA
30 CALL FCSTRS (NA,ICGA,ISYMA,ITYPEA)
ICGB(1) = 0
DO 40 JSB=1,ISB
40 CALL FCSTRS (NB,ICGB,ISYMB,ITYPEA)
RETURN
END
SUBROUTINE FCSTRS (N,ICG,ISYM,ITYPE)
C... GENERATE A NEW ALPHA OR BETA STRING WITH SYMMETRY ISYM
C... ICG SHOULD BE PRESERVED BETWEEN CALLS
C... IF (ISYM.EQ.0) ON ENTRY, THEN JUST NEXT STRING GENERATED
C... ITS SYMMETRY RETURNED IN ISYM
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
DIMENSION ICG(*),ITYPE(*)
C
IF (ICG(1).LE.0) THEN
DO 10 I=1,N
10 ICG(I)=I
IF (N.LE.0) THEN
DO 111 JJ=1,NACT
IF (ITYPE(JJ).EQ.1) GOTO 112
111 CONTINUE
112 ICG(1) = JJ
END IF
GOTO 40
END IF
C
21 DO 20 IU=N,1,-1
ICG(IU) = ICG(IU) + 1
IF (ICG(IU).LE.NACT+IU-N) GOTO 30
20 CONTINUE
GOTO 60
CDIR$ NEXTSCALAR
*VOPTION NOVEC
30 DO 31 JU=IU+1,N
31 ICG(JU) = ICG(IU)+JU-IU
C
40 JMT = ITYPE(ICG(1))
DO 50 J=2,N
50 JMT = MULT(JMT,ITYPE(ICG(J)))
IF (JMT.NE.ISYM.AND.ISYM.NE.0) GOTO 21
IF (ISYM.EQ.0) ISYM = JMT
RETURN
60 ICG(1) = 0
RETURN
END
FUNCTION IBINOM (M,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
TOP = M
NN = MIN0(M-N,N)
IF (NN.LT.0) GOTO 20
BOT = NN
BINOM = 1.0D0
DO 10 I=1,NN
BINOM = BINOM * TOP/BOT
TOP = TOP - 1.0D0
10 BOT = BOT - 1.0D0
IBINOM = BINOM + 0.1D0
IF (DABS(BINOM-DBLE(IBINOM)).GT.0.5D0)
1 STOP 'ARITHMETIC PROBLEM IN COMPUTING BINOMIAL COEFFICIENT'
RETURN
20 IBINOM = 0
RETURN
END
SUBROUTINE INITI (INTER,NITEM,NBOX)
CSUBR SETS UP PARTIAL WEIGHT ARRAY FOR ADDRESSING BINOMIAL DISTRIBUTIONS
DIMENSION INTER(NITEM,*)
N1=NITEM+1
DO 30 K=1,NITEM
DO 10 L=1,NBOX
10 INTER(K,L)=0
DO 20 L=K,NBOX-NITEM+K-1
INTER(K,L+1) = IBINOM(NBOX-L,NITEM-K)+INTER(K,L)
20 CONTINUE
30 CONTINUE
DO 40 K=1,NITEM-1
DO 40 L=K,NBOX-NITEM+K
40 INTER(K,L) = INTER(K,L) - INTER(K+1,L+1)
DO 50 L=NITEM,NBOX
50 INTER(NITEM,L) = L-NITEM
RETURN
END
SUBROUTINE FZERO (A,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DIMENSION A(*)
DO 1 I=1,N
1 A(I) = 0D0
RETURN
END
SUBROUTINE OUTSQR (Q,IDIM,IA,IB,TITLE)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
CHARACTER*(*) TITLE
DIMENSION Q(IDIM,1)
CHARACTER*10 FORMAT
DATA FORMAT /'(/1X,A000)'/
WRITE (FORMAT(7:9),'(I3)') 55+LEN(TITLE)/2
WRITE (6,FORMAT) TITLE
M=1
N=8
10 IF (IB.LT.M) RETURN
N=MIN0(N,IB)
WRITE (6,30) (I,I=M,N)
WRITE (6,40)
DO 20 J=1,IA
WRITE (6,50) J,(Q (J,I) ,I=M,N)
20 CONTINUE
M=M+8
N=N+8
GOTO 10
30 FORMAT(//3X,8I14)
40 FORMAT(/)
50 FORMAT(7X,I3,8F14.7)
END
SUBROUTINE OUTVEC (P,IB,TITLE)
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
COMMON /TAPES / INP,IOUT
CHARACTER*(*) TITLE
DIMENSION P(1)
WRITE (IOUT,'(/1X,(A))') TITLE
M=1
N=10
10 IF (IB.LT.M) RETURN
N=MIN0(N,IB)
WRITE(IOUT,20) M,N,(P (I) ,I=M,N)
M=M+10
N=N+10
GOTO 10
20 FORMAT(I4,'-',I3,10F12.6)
END
SUBROUTINE DIAG2(M,N,D,X)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (MAXDIM=200)
DIMENSION D(M), X(M,M)
DIMENSION E(MAXDIM)
IF(M.GT.MAXDIM) THEN
WRITE(6,10) M,MAXDIM
10 FORMAT(' DIMENSION TOO LARGE IN DIAG2:',2I8)
STOP
END IF
CSTART VAX UNIX-HP UNIX-CONVEX CRAY IBM UNIVAC
EPS=2.5E-13
TOL=2.5D-18
CEND
C
IF(N.EQ.1) GO TO 400
C
C HOUSEHOLDER'S REDUCTION
C SIMULATION OF LOOP DO 150 I=N,2,(-1)
C
DO 150 NI=2,N
I=N+2-NI
L=I-2
H=0.0
G=X(I,I-1)
IF(L) 140,140,20
20 DO 30 K=1,L
30 H=H+X(I,K)**2
S=H+G*G
IF(S.GE.TOL) GO TO 50
40 H=0.0
GO TO 140
50 IF(H) 140,140,60
60 L=L+1
F=G
G=DSQRT(S)
IF(F) 75,75,70
70 G=-G
75 H=S-F*G
X(I,I-1)=F-G
F=0.0
C
DO 110 J=1,L
X(J,I)=X(I,J)/H
S=0.0
DO 80 K=1,J
80 S=S+X(J,K)*X(I,K)
J1=J+1
IF(J1.GT.L) GO TO 100
DO 90 K=J1,L
90 S=S+X(K,J)*X(I,K)
100 E(J)=S/H
110 F=F+S*X(J,I)
C
F=F/(H+H)
C
DO 120 J=1,L
120 E(J)=E(J)-F*X(I,J)
C
DO 130 J=1,L
F=X(I,J)
S=E(J)
DO 130 K=1,J
130 X(J,K)=X(J,K)-F*E(K)-X(I,K)*S
C
140 D(I)=H
150 E(I-1)=G
C
C ACCUMULATION OF TRANSFORMATION MATRICES
C
160 D(1)=X(1,1)
X(1,1)=1.0
DO 220 I=2,N
L=I-1
IF(D(I)) 200,200,170
170 DO 190 J=1,L
S=0.0
DO 180 K=1,L
180 S=S+X(I,K)*X(K,J)
DO 190 K=1,L
190 X(K,J)=X(K,J)-S*X(K,I)
200 D(I)=X(I,I)
X(I,I)=1.0
210 DO 220 J=1,L
X(I,J)=0.0
220 X(J,I)=0.0
C
C DIAGONALIZATION OF THE TRIDIAGONAL MATRIX
C
B=0.0
F=0.0
E(N)=0.0
C
DO 340 L=1,N
H=EPS*(DABS(D(L))+DABS(E(L)))
IF (H.GT.B) B=H
C
C TEST FOR SPLITTING
C
DO 240 J=L,N
IF (DABS(E(J)).LE.B) GOTO 250
240 CONTINUE
C
C TEST FOR CONVERGENCE
C
250 IF(J.EQ.L) GO TO 340
C
C SHIFT FROM UPPER 2*2 MINOR
C
260 P=(D(L+1)-D(L))*0.5/E(L)
R=DSQRT(P*P+1.0)
IF(P) 270,280,280
270 P=P-R
GO TO 290
280 P=P+R
290 H=D(L)-E(L)/P
DO 300 I=L,N
300 D(I)=D(I)-H
F=F+H
C
C QR TRANSFORMATION
C
P=D(J)
C=1.0
S=0.0
C
C SIMULATION OF LOOP DO 330 I=J-1,L,(-1)
C
J1=J-1
DO 330 NI=L,J1
I=L+J1-NI
G=C*E(I)
H=C*P
C
C PROTECTION AGAINST UNDERFLOW OF EXPONENTS
C
IF (DABS(P).LT.DABS(E(I))) GOTO 310
C=E(I)/P
R=DSQRT(C*C+1.0)
E(I+1)=S*P*R
S=C/R
C=1.0/R
GO TO 320
310 C=P/E(I)
R=DSQRT(C*C+1.0)
E(I+1)=S*E(I)*R
S=1.0/R
C=C/R
320 P=C*D(I)-S*G
D(I+1)=H+S*(C*G+S*D(I))
DO 330 K=1,N
H=X(K,I+1)
X(K,I+1)=X(K,I)*S+H*C
330 X(K,I)=X(K,I)*C-H*S
C
E(L)=S*P
D(L)=C*P
IF (DABS(E(L)).GT.B) GO TO 260
C
C CONVERGENCE
C
340 D(L)=D(L)+F
C
C ORDERING OF EIGENVALUES
C
NI=N-1
350 DO 380I=1,NI
K=I
P=D(I)
J1=I+1
DO 360J=J1,N
IF(D(J).GE.P) GOTO 360
K=J
P=D(J)
360 CONTINUE
IF (K.EQ.I) GOTO 380
D(K) =D(I)
D(I)=P
DO 370 J=1,N
P=X(J,I)
X(J,I)=X(J,K)
370 X(J,K)=P
380 CONTINUE
C
C FIXING OF SIGN
C
DO 385 I=1,N
PM=0
DO 386 J=1,N
IF(PM.GT.DABS(X(J,I))) GOTO 386
PM =DABS(X(J,I))
K=J
386 CONTINUE
IF(X(K,I).GE.0) GOTO 385
DO 387 J=1,N
387 X(J,I)=-X(J,I)
385 CONTINUE
390 GO TO 410
C
C SPECIAL TREATMENT OF CASE N = 1
C
400 D(1)=X(1,1)
X(1,1)=1.0
410 RETURN
END
SUBROUTINE TPRINT (TITLE,NOP,T,STITLE,NS)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
COMMON /TAPES / INP,IOUT
COMMON /CFULLR/ CORE,THR,THRRES
COMMON /CFULL / MULT(8,8),ICF(256),NPAIR(8),NSTRA(8),NSTRB(8)
> ,INTOFF(8),ITYPEA(256)
> ,NOCC(8),NCI(8),NELEC,MS2,NA,NB,NACT
> ,INTAA,INTHAA,INTJAA,INTKAA,NINT2
> ,NSTRAA,NSTRBB,MAXAA,MAXBB,MAXREP,MAXPAR,INTADT
> ,IPRINT,IPRDIA,IPRHAM,IPRTIM,MAXIT,NROOT,MEMORY
CHARACTER*(*) TITLE,STITLE
FLOP = 0
IF (T.NE.0D0) FLOP = 1D-6 * NOP/T
IF (STITLE.EQ.' ') THEN
WRITE (IOUT,666) TITLE,NOP,T,FLOP
ELSE
WRITE (IOUT,666) TITLE,NOP,T,FLOP,STITLE,NS
END IF
666 FORMAT(' Timing: ',A,5X,'Operations',I12,5X,'Time',F10.2,
>5X,'Mflops',F10.3,:,5X,A,I6)
RETURN
END
Reference in New Issue Copy Permalink