C-*- Fortran -*-
PROGRAM EPSNES
C
C molecule, from 1D ring
C
C -----> all integrals (ov|ov) in-core <------
C
C perturbation theory, MP2, MP3, Epstein-Nesbet 2nd and 3rd order
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C ORTHO - programs for ab-initio calculations in localized orbitals
C Copyright (C) 2008 Peter Reinhardt (Univ. Paris VI, France)
C
C This program is free software: you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation, either version 3 of the License, or
C (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public License
C along with this program. If not, see .
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C..FILE 'flow.h'
C..FILE 'counters.h'
C..FILE 'nbull.h'
C..FILE 'common_syst.h'
C..FILE 'common_intu.h'
C..FILE 'common_vect.h'
C..FILE 'common_denom.h'
C..FILE 'common_wbuf.h'
C
C..DELCF WE HAVE A ROUTINE FLU
C DELCD Debug
C DELCT test conditions
C
C if we want to do MP3, we have to store integrals (oo|oo), (oo|vv), and
C (vv|vv) as well
C
C
INCLUDE 'param.h'
PARAMETER (NBULL=43000000)
COMMON /TWOI/ H0(NBULL),IH0(4,NBULL),NUMINT
C.. INCLUDE 'nbull.h'
COMMON /FLOW/ THRESH,ICUT,NFROZ,LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX
$ ,L3RDF,L4THF,LINFF,LCUT,LTHRE,LMART,L34EPS,LMP3
LOGICAL LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX,L3RDF,L4THF,LINFF,LCUT
$ ,LTHRE,LMART,L34EPS,LMP3
C.. INCLUDE 'flow.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCCS(NBASM),IOCC(NBASM)
C.. INCLUDE 'common_vect.h'
C
CHARACTER*8 PNAME
LOGICAL LHIGH,LDUMP,L3
PNAME='EPSNES '
INCLUDE 'compiler_stamp'
WRITE(6,*)
WRITE(6,*) ' E P S N E S - single molecule'
WRITE(6,*) ' correlation by MP2 and Epstein-Nesbet perturbation'
WRITE(6,*) ' MP3 as well possible '
WRITE(6,*)
C
X=CPTIME(3)
CALL DATING(PNAME,1)
L3=.FALSE.
LHIGH=.FALSE.
LDUMP=.FALSE.
C
IUNITR=26
OPEN(UNIT=IUNITR,FILE='SYSTEM.ORTHO',STATUS='OLD',
- FORM='FORMATTED',ERR=901)
READ(IUNITR,*) NATOM
WRITE(6,*)
NSHL=0
NBAS=1
LMAX=0
DO IAT=1,NATOM
READ(IUNITR,*) IDUM,NSHLAT
DO ISH=1,NSHLAT
READ(IUNITR,*) ITYPE,NPRIM
LMAX=MAX(LMAX,ITYPE+1)
NBAS=NBAS+ITYPE+ITYPE+1
DO III=1,NPRIM
READ(IUNITR,*) XDUM
END DO
END DO
END DO
NBAS=NBAS-1
WRITE(6,*)
WRITE(6,*) ' READ INFORMATION ON SYSTEM'
WRITE(6,*) ' NATOMS ',NATOM
WRITE(6,*) ' NBAS ',NBAS
WRITE(6,*)
WRITE(6,*)
WRITE(6,*) ' COMPILED MAXIMUM DIMENSIONS'
WRITE(6,*) ' NBAS ',NBASM
WRITE(6,*)
IF (NBAS.GT.NBASM) THEN
WRITE(6,*) 'NBAS, MAXIMUM: ',NBAS,NBASM
STOP 'INCREASE NBASM IN param.h AND RECOMPILE'
END IF
CALL RDINP
C
C read vector
C
IUNITO=31
WRITE(6,*)
WRITE(6,*) ' READING VECTOR FROM UNIT',IUNITO
OPEN(UNIT=IUNITO,FILE='VECTOR',FORM='FORMATTED',STATUS='OLD',ERR
$ =902)
READ(IUNITO,*) ((CI(I,J),I=1,NBAS),J=1,NBAS)
READ(IUNITO,*) (IOCC(I),I=1,NBAS)
READ(IUNITO,*) (IOCCS(I),I=1,NBAS)
CLOSE(IUNITO)
C
C CALL PFUNC(CI,NBAS,NCELL)
C
NOCC=0
DO I=1,NBAS
IF(IOCCS(I).EQ.IOCC(I)) NOCC=NOCC+1
END DO
NVIRT=NBAS-NOCC
WRITE(6,*) ' NOCC ',NOCC
WRITE(6,*)
NVIRT=NBAS-NOCC
C
IF (LUNFOR) THEN
WRITE(6,*) ' READING UNFORMATTED INTEGRALS, BUFFER LENGTH '
$ ,IWBULL
ELSE
WRITE(6,*) ' READING INTEGRALS FORMATTED'
END IF
WRITE(6,*) ' THRESHOLD FOR BIELECTRONIC INTEGRALS: ',THRESH
IF (LMP2PR) WRITE(6,*) ' ALL MP2 INTEGRALS WILL BE PRINTED '
IF (LCUT) WRITE(6,*) ' CUT-OFF RADIUS SET TO ',ICUT
C
IF (L3RDF) WRITE(6,*) ' 3rd order localization diagrams '
IF (L4THF) WRITE(6,*) ' 4th order paired F '
IF (L34EPS) WRITE(6,*) ' 3rd and 4th order F in Epstein-Nesbet '
IF (LINFF) WRITE(6,*)
$ ' infinite summation over pairs (geometrical series) '
IF (LPRCEX) WRITE(6,*) ' OUTPUT: J, K '
IF (LPRIFO) WRITE(6,*) ' Fock matrix will be printed '
WRITE(6,*)
C
C read Fock matrix
C
IUNITF=23
OPEN(UNIT=IUNITF,FILE='FOCK',FORM='FORMATTED',STATUS='OLD',ERR
$ =904)
READ(IUNITF,*) IDUM1
WRITE(6,*) ' READING FOCK MATRIX FROM UNIT',IUNITF
WRITE(6,*)
READ(IUNITF,*) ((F(I,J),I=1,NBAS),J=1,NBAS)
READ(IUNITF,*) EDUM,EN,E1,E2
CLOSE(IUNITF)
WRITE(6,*) ' E0 E1 E2 ',EN,E1,E2
C
C F in MOs
C
CALL TRANSF(F,CI)
C
WRITE(6,*)
WRITE(6,*) ' TRANSFORMED FOCK MATRIX:'
WRITE(6,*)
WRITE(6,*)' -----------------------------------------------'
WRITE(6,*)
IF (LPRIFO) CALL PFUNC(F,NBAS)
DO I=1,NBAS
ORBEN(I)=F(I,I)
END DO
WRITE(6,*) ' THE ORBITAL ENERGIES:'
WRITE(6,*) ' occupied:'
WRITE(6,'(4(I4,F20.12))') (I,ORBEN(I),I=1,NOCC)
WRITE(6,*) ' virtual:'
WRITE(6,'(4(I4,F20.12))') (I,ORBEN(I),I=1+NOCC,NBAS)
WRITE(6,*)
C
IUNIT4=53
CALL READ53(IUNIT4)
C
C calculate two-electron energy
C
CALL TOTCAL(E1,E2,EN)
C
C MP2, EN
C
CALL PERTET
CALL TIMING('PERT')
C
C MP3
C
IF (LMP3) THEN
CALL MP3
c$$$ CALL MP3BIS
CALL TIMING('MP3 ')
END IF
C
WRITE(6,*)
X=CPTIME(4)
CALL DATING(PNAME,2)
STOP
901 CONTINUE
WRITE(6,*) ' NO FILE FOUND ... '
STOP
902 CONTINUE
WRITE(6,*) ' NO FILE FOUND ... '
STOP
904 CONTINUE
WRITE(6,*) ' NO FILE FOUND ... '
STOP
END
C
SUBROUTINE PERTET
INCLUDE 'param.h'
PARAMETER (NBULL=43000000)
COMMON /TWOI/ H0(NBULL),IH0(4,NBULL),NUMINT
C.. INCLUDE 'nbull.h'
COMMON /FLOW/ THRESH,ICUT,NFROZ,LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX
$ ,L3RDF,L4THF,LINFF,LCUT,LTHRE,LMART,L34EPS,LMP3
LOGICAL LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX,L3RDF,L4THF,LINFF,LCUT
$ ,LTHRE,LMART,L34EPS,LMP3
C.. INCLUDE 'flow.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCCS(NBASM),IOCC(NBASM)
C.. INCLUDE 'common_vect.h'
PARAMETER (N2MX=NBASM*NBASM,N3MX=NBASM*N2MX,N4MX=N3MX*NBASM,NDMX
$ =1+N4MX/64+N3MX/32+N2MX/16)
C COMMON /DENOM/ DENSUM(NDMX)
C.. INCLUDE 'common_denom.h'
LOGICAL LIJ,LKL,L8
C
IF (LINFF.AND.L3RDF) L8=.TRUE.
CD WRITE(6,*) ' PERTET: LCUT, ICH, ICL ',LCUT,ICUH,ICUL
CD WRITE(6,*) ' PERTET: NOCC, NVIRT, NBAS ',NOCC,NVIRT,NBAS
CD WRITE(6,*) ' L3RDF ',L3RDF
CD WRITE(6,*) ' L4THF ',L4THF
CD WRITE(6,*) ' LINFF ',LINFF
CD WRITE(6,*) ' L33EPS ',L34EPS
IF (LINFF) THEN
NOV=NOCC+1
DO I=1,NOCC
EI=ORBEN(I)
DO J=1,I
EIJ=EI+ORBEN(J)
EJ=ORBEN(J)
DO K=NOV,NBAS
EIJK=EIJ-ORBEN(K)
EK=ORBEN(K)
DO L=NOV,K
EL=ORBEN(L)
EIJKL=EIJK-ORBEN(L)
T1=0.D0
T2=0.D0
T3=0.D0
T4=0.D0
DO M=1,NOCC
IF (M.NE.I) THEN
FF=F(I,M)
FF=FF*FF/(EIJKL-EI+ORBEN(M))
T1=T1+FF
END IF
IF (M.NE.J) THEN
FF=F(J,M)
FF=FF*FF/(EIJKL-EJ+ORBEN(M))
T2=T2+FF
END IF
END DO
DO M=NOV,NBAS
IF (M.NE.K) THEN
FF=F(K,M)
FF=FF*FF/(EIJKL+EK-ORBEN(M))
T3=T3+FF
END IF
IF (M.NE.L) THEN
FF=F(L,M)
FF=FF*FF/(EIJKL+EL-ORBEN(M))
T4=T4+FF
END IF
END DO
C
CALL FABINX(I,J,K,L,INDX)
C DENSUM(INDX)=EIJKL-T1-T2-T3-T4
C WRITE(6,9883) I,J,K,L,EIJKL,T1,T2,T3,T4,DENSUM(I,J,K,L)
9883 FORMAT('DENSUM:',4I3,F9.6,F11.6,3F9.6,F11.6)
END DO
END DO
END DO
END DO
END IF
C
C loop over determinants
C
C
C MP2 contributions
C
EM1=0.D0
EM2=0.D0
EM3=0.D0
EM4=0.D0
C
C 3rd order
C
EM13=0.D0
EM23=0.D0
EM33=0.D0
EM43=0.D0
C
C 4th order
C
EM14=0.D0
EM24=0.D0
EM34=0.D0
EM44=0.D0
C
C Martin Albrecht's 3rd order + summation
C
EM15=0.D0
EM25=0.D0
EM35=0.D0
EM45=0.D0
C
C Martin Albrecht's 4th order + summation
C
EM16=0.D0
EM26=0.D0
EM36=0.D0
EM46=0.D0
C
C my infinite summation in pairs
C
EM17=0.D0
EM27=0.D0
EM37=0.D0
EM47=0.D0
C
C my summation onto 3rd order
C
EM18=0.D0
EM28=0.D0
EM38=0.D0
EM48=0.D0
C
C Epstein-Nesbet contributions
C
EE1=0.D0
EE2=0.D0
EE3=0.D0
EE4=0.D0
C
EE1S=0.D0
EE2S=0.D0
EE3S=0.D0
EE4S=0.D0
C
EE13=0.D0
EE23=0.D0
EE33=0.D0
EE43=0.D0
C
EE13S=0.D0
EE23S=0.D0
EE33S=0.D0
EE43S=0.D0
C
EE14=0.D0
EE24=0.D0
EE34=0.D0
EE44=0.D0
C
EE14S=0.D0
EE24S=0.D0
EE34S=0.D0
EE44S=0.D0
C
C
XTT7=0.D0
XTT17=0.D0
C
NOV=NOCC+1
DO I=1,NOCC
EI=ORBEN(I)
DO J=I,NOCC
INDJ=J
EJ=ORBEN(J)
IF (I.EQ.INDJ) THEN
LIJ=.TRUE.
ELSE
LIJ=.FALSE.
END IF
DO K=NOV,NBAS
EK=ORBEN(K)
INDK=K
DO L=K,NBAS
INDL=L
EL=ORBEN(L)
IF (INDK.EQ.INDL) THEN
LKL=.TRUE.
ELSE
LKL=.FALSE.
END IF
C WRITE(6,*) ' I, INDJ, INDK, INDL ',I,INDJ,INDK,INDL
C
C MP2
C
IF (LIJ) THEN
IF (LKL) THEN
EKI=EI+EI-EK-EK
HHH=HEXC(I,K)
XNOM=HHH*HHH
C ii -> kk
EM1=EM1+XNOM/EKI
IF (LINFF) THEN
CALL FABINX(I,J,K,L,INDX)
C EM17=EM17+XNOM/DENSUM(INDX)
END IF
C Epstein - Nesbet
ENKI=EKI-HCOU(I,I)-HCOU(K,K)+
- 4.D0*HCOU(I,K)-2.D0*HEXC(I,K)
EE1=EE1+XNOM/ENKI
C --- start L3, L4
IF (L3RDF.OR.L4THF.OR.L34EPS) THEN
XTT=HHH/EKI
EM131=0.D0
EM141=0.D0
EM151=0.D0
EM161=0.D0
EM181=0.D0
EE131=0.D0
EE131S=0.D0
EE141=0.D0
EE141S=0.D0
DO II=1,NOCC
INDII=II
IF (INDII.NE.I) THEN
H3=HFIND(I,INDK,INDII,INDK)
EII=ORBEN(II)
EDEN1=EKI-EI+EII
FF=F(I,II)
EM131=EM131-H3*FF/EDEN1
EM141=EM141+2.D0*FF*FF/EDEN1
EM151=EM151-H3*FF/(EDEN1-FF*FF/EKI)
EM161=EM161+2.D0*FF*FF/(EDEN1*(1.D0-FF*FF/EDEN1
$ /EKI))
IF (L8) THEN
CALL FABINX(II,I,K,K,INDX)
C EM181=EM181-H3*FF/DENSUM(INDX)
END IF
C
C the right Epstein-Nesbet denominators:
C I,II,K,K -> EDC3N
IF (L34EPS) THEN
CALL EDC3N(I,II,K,EDEN1,E,ES)
EE131=EE131-H3*FF/E
EE141=EE141+2.D0*FF*FF/E
EE131S=EE131S-H3*FF/ES
EE141S=EE141S+2.D0*FF*FF/ES
END IF
END IF
END DO
DO KK=NOV,NBAS
INDKK=KK
IF (INDKK.NE.INDK) THEN
H3=HFIND(I,INDK,I,INDKK)
EKK=ORBEN(KK)
EDEN2=EKI+EK-EKK
FF=F(K,KK)
EM131=EM131+H3*FF/EDEN2
EM141=EM141+2.D0*FF*FF/EDEN2
EM151=EM151+H3*FF/(EDEN2-FF*FF/EKI)
EM161=EM161+2.D0*FF*FF/(EDEN2*(1.D0-FF*FF/EDEN2
$ /EKI))
IF (L8) THEN
CALL FABINX(I,I,KK,K,INDX)
C EM181=EM181+H3*FF/DENSUM(INDX)
END IF
C Eps-Nes
IF (L34EPS) THEN
CALL EDC2N(I,K,KK,EDEN2,E,ES)
EE131=EE131+H3*FF/E
EE141=EE141+2.D0*FF*FF/E
EE131S=EE131S+H3*FF/ES
EE141S=EE141S+2.D0*FF*FF/ES
END IF
END IF
END DO
EM13=EM13+XTT*EM131
EM14=EM14+XNOM*EM141/(EKI*EKI)
EM15=EM15+XTT*EM151
EM16=EM16+XNOM*EM161/(EKI*EKI)
IF (L8) EM18=EM18+XTT*EM181
EE13=EE13+XTT*EE131
EE14=EE14+XNOM*EE141/(ENKI*ENKI)
EE13S=EE13S+XTT*EE131S
EE14S=EE14S+XNOM*EE141S/(ENKI*ENKI)
END IF
C --- end L3, L4
ELSE
EKLI=EI+EI-EK-EL
HIJKL=HFIND(I,INDK,I,INDL)
XTT=HIJKL/EKLI
C ii -> kl
EM2=EM2+HIJKL*XTT
C Epstein-Nesbet
C
CALL EDC2N(I,K,L,EKLI,ENKLI,ENKLIS)
EE2=EE2+HIJKL*HIJKL/ENKLI
EE2S=EE2S+HIJKL*HIJKL/ENKLIS
IF (LINFF) THEN
CALL FABINX(I,I,K,L,INDX)
C XTT7=HIJKL/DENSUM(INDX)
EM27=EM27+HIJKL*XTT7
END IF
C --- start L3, L4
IF (L3RDF.OR.L4THF.OR.L34EPS) THEN
EM231=0.D0
EM241=0.D0
EM251=0.D0
EM261=0.D0
EM281=0.D0
EE231=0.D0
EE231S=0.D0
EE241=0.D0
EE241S=0.D0
DO II=1,NOCC
INDII=II
IF (INDII.NE.I) THEN
EII=ORBEN(II)
H3=HFIND(I,INDK,INDII,INDL)
H4=HFIND(I,INDL,INDII,INDK)
EDEN1=EKLI-EI+EII
FF=F(I,II)
EM231=EM231-(H3+H4)*FF/EDEN1
EM241=EM241+2.D0*FF*FF/EDEN1
EM251=EM251-(H3+H4)*FF/(EDEN1-FF*FF/EKLI)
EM261=EM261+2.D0*FF*FF/(EDEN1*(1.D0-FF*FF/EDEN1
$ /EKLI))
IF (L8) THEN
CALL FABINX(I,II,K,L,INDX)
C EM281=EM281-(H3+H4)*FF/DENSUM(INDX)
END IF
C Eps-Nes denom.
C case 4: I,II,K,L
IF (L34EPS) THEN
CALL EDC4N(I,II,K,L,H3,H4,EDEN1,E,ES)
EE231=EE231-(H3+H4)*FF/E
EE241=EE241+2.D0*FF*FF/E
EE231S=EE231S-(H3+H4)*FF/ES
EE241S=EE241S+2.D0*FF*FF/ES
END IF
END IF
END DO
DO KK=NOV,NBAS
INDKK=KK
IF (INDKK.NE.INDK) THEN
EKK=ORBEN(KK)
EDEN2=EKLI+EK-EKK
H3=HFIND(I,INDKK,I,INDL)
FF=F(KK,K)
EM231=EM231+H3*FF/EDEN2
EM241=EM241+FF*FF/EDEN2
EM251=EM251+H3*FF/(EDEN2-FF*FF/EKLI)
EM261=EM261+FF*FF/(EDEN2*(1.D0-FF*FF/EKLI/EDEN2))
IF (L8) THEN
CALL FABINX(I,I,KK,L,INDX)
C EM281=EM281+H3*FF/DENSUM(INDX)
END IF
IF (L34EPS) THEN
C Eps-Nes
IF (INDKK.NE.INDLL) THEN
CALL EDC2N(I,KK,L,EDEN2,E,ES)
ELSE
E=EDEN2-HCOU(I,I)-HCOU(L,L)+
- 4.D0*HCOU(I,L)-2.D0*HEXC(I,L)
ES=E
END IF
EE231=EE231+H3*FF/E
EE241=EE241+FF*FF/E
EE231S=EE231S+H3*FF/ES
EE241S=EE241S+FF*FF/ES
END IF
END IF
IF (INDKK.NE.INDL) THEN
EKK=ORBEN(KK)
EDEN3=EKLI+EL-EKK
FF=F(KK,L)
H3=HFIND(I,INDK,I,INDKK)
EM231=EM231+H3*FF/EDEN3
EM241=EM241+FF*FF/EDEN3
EM251=EM251+H3*FF/(EDEN3-FF*FF/EKLI)
EM261=EM261+FF*FF/(EDEN3*(1.D0-FF*FF/EKLI/EDEN3))
IF (L8) THEN
CALL FABINX(I,I,K,KK,INDX)
C EM281=EM281+H3*FF/DENSUM(INDX)
END IF
IF (L34EPS) THEN
C Eps-Nes
IF (INDKK.NE.INDK) THEN
CALL EDC2N(I,K,KK,EDEN3,E,ES)
ELSE
E=EDEN3-HCOU(I,I)-HCOU(KK,KK)+
- 4.D0*HCOU(I,KK)-2.D0*HEXC(I,KK)
ES=E
END IF
EE231=EE231+H3*FF/E
EE241=EE241+FF*FF/E
EE231S=EE231S+H3*FF/ES
EE241S=EE241S+FF*FF/ES
END IF
END IF
END DO
EM23=EM23+XTT*EM231
EM24=EM24+EM241*XTT*XTT
EM25=EM25+XTT*EM251
EM26=EM26+EM261*XTT*XTT
IF (L8) EM28=EM28+XTT7*EM281
EE23=EE23+HIJKL/ENKLI*EE231
EE24=EE24+EM241*(HIJKL/ENKLI)*(HIJKL/ENKLI)
EE23S=EE23S+HIJKL/ENKLIS*EE231S
EE24S=EE24S+EE241S*(HIJKL/ENKLIS)*(HIJKL/ENKLIS)
END IF
C --- end L3, L4
END IF
ELSE
IF (LKL) THEN
EKKIJ=EJ+EI-EK-EK
HIJKL=HFIND(I,INDK,INDJ,INDK)
XTT=HIJKL/EKKIJ
C ij -> kk
C
c$$$ IF (ABS(XTT*HIJKL).GT.1.D-10) WRITE(6,9335) I,INDK
c$$$ $ ,INDJ,INDK,EKKIJ
c$$$ 9335 FORMAT(' EM3: I,INDK,INDJ,INDK,TERM ',4I5,2F20.12)
C ij -> kk
EM3=EM3+XTT*HIJKL
EM37=EM37+XTT7*HIJKL
C Epstein-Nesbet
C
CALL EDC3N(I,J,K,EKKIJ,ENKIJ,ENKIJS)
EE3=EE3+HIJKL*HIJKL/ENKIJ
EE3S=EE3S+HIJKL*HIJKL/ENKIJS
IF (LINFF) THEN
CALL FABINX(I,J,K,K,INDX)
C XTT7=HIJKL/DENSUM(INDX)
END IF
C --- start L3, L4
IF (L3RDF.OR.L4THF.OR.L34EPS) THEN
EM331=0.D0
EM341=0.D0
EM351=0.D0
EM361=0.D0
EM381=0.D0
EE331=0.D0
EE331S=0.D0
EE341=0.D0
EE341S=0.D0
DO II=1,NOCC
INDII=II
IF (INDII.NE.I) THEN
EII=ORBEN(II)
EDEN1=EKKIJ-EI+EII
INDIII=INDII
INDJJ=INDJ
INDKK=INDK
C CALL TURN4(INDJ,INDII,INDK,IVJ,INDJJ,INDIII,INDKK)
H3=HFIND(J,INDKK,INDIII,INDKK)
FF=F(I,II)
EM331=EM331-H3*FF/EDEN1
EM341=EM341+FF*FF/EDEN1
EM351=EM351-H3*FF/(EDEN1-FF*FF/EKKIJ)
EM361=EM361+FF*FF/(EDEN1*(1.D0-FF*FF/EKKIJ/EDEN1
$ ))
IF (L8) THEN
CALL FABINX(II,J,K,K,INDX)
C EM381=EM381-H3*FF/DENSUM(INDX)
END IF
C Eps-Nes
IF (INDII.NE.INDJ) THEN
CALL EDC3N(II,J,K,EDEN1,E,ES)
ELSE
E=EDEN1-HCOU(J,J)-HCOU(K,K)+
- 4.D0*HCOU(J,K)-2.D0*HEXC(J,K)
ES=E
END IF
EE331=EE331-H3*FF/E
EE341=EE341+FF*FF/E
EE331S=EE331S-H3*FF/ES
EE341S=EE341S+FF*FF/ES
END IF
IF (INDII.NE.INDJ) THEN
H3=HFIND(I,INDK,INDII,INDK)
FF=F(II,J)
EDEN2=EKKIJ-EJ+EII
EM331=EM331-H3*FF/EDEN2
EM341=EM341+FF*FF/EDEN2
EM351=EM351-H3*FF/(EDEN2-FF*FF/EKKIJ)
EM361=EM361+FF*FF/(EDEN2*(1.D0-FF*FF/EKKIJ/EDEN2))
IF (L8) THEN
CALL FABINX(I,II,K,K,INDX)
C EM381=EM381-H3*FF/DENSUM(INDX)
END IF
C Eps-Nes
IF (INDII.NE.I) THEN
CALL EDC3N(I,II,K,EDEN2,E,ES)
ELSE
E=EDEN2-HCOU(I,I)-HCOU(K,K)+
- 4.D0*HCOU(I,K)-2.D0*HEXC(I,K)
ES=E
END IF
EE331=EE331-H3*FF/E
EE341=EE341+FF*FF/E
EE331S=EE331S-H3*FF/ES
EE341S=EE341S+FF*FF/ES
END IF
END DO
DO KK=NOV,NBAS
INDKK=KK
IF (INDKK.NE.INDK) THEN
EKK=ORBEN(KK)
EDEN2=EKKIJ+EK-EKK
FF=F(KK,K)
H3=HFIND(I,INDKK,INDJ,INDK)
H4=HFIND(I,INDK,INDJ,INDKK)
EM331=EM331+(H3+H4)*FF/EDEN2
EM341=EM341+2.D0*FF*FF/EDEN2
IF (L8) THEN
CALL FABINX(I,J,KK,K,INDX)
C EM381=EM381+(H3+H4)*FF/DENSUM(INDX)
END IF
IF (L34EPS) THEN
C Eps-Nes
CALL EDC4N(I,J,K,KK,H3,H4,EDEN2,E
$ ,ES)
EE331=EE331+(H3+H4)*FF/E
EE341=EE341+2.D0*FF*FF/E
EE331S=EE331S+(H3+H4)*FF/ES
EE341S=EE341S+2.D0*FF*FF/ES
END IF
END IF
END DO
EM33=EM33+XTT*EM331
EM34=EM34+EM341*XTT*XTT
EM35=EM35+XTT*EM351
EM36=EM36+EM361*XTT*XTT
IF (L8) EM38=EM38+XTT7*EM381
IF (L34EPS) THEN
EE33=EE33+HIJKL*HIJKL/ENKIJ*EE331
EE34=EE34+EE341*(HIJKL/ENKIJ)*(HIJKL/ENKIJ)
EE33S=EE33S+HIJKL/ENKIJS*EE331S
EE34S=EE34S+EE341S*(HIJKL/ENKIJS)*(HIJKL/ENKIJS)
END IF
END IF
C --- end L3, L4
ELSE
EKLIJ=EI+EJ-EK-EL
H1=HFIND(I,INDK,INDJ,INDL)
H2=HFIND(I,INDL,INDJ,INDK)
H11=H1*H1
H22=H2*H2
H12=H1*H2
XNOM=(H11+H22-H12)
C ij -> kl
EM4=EM4+XNOM/EKLIJ
IF (LINFF) THEN
CALL FABINX(I,J,K,L,INDX)
C EM47=EM47+XNOM/DENSUM(INDX)
END IF
C Epstein-Nesbet
C
CALL EDC4N(I,J,K,L,H1,H2,EKLIJ,YKLIJ
$ ,YKLIJS)
EE4=EE4+XNOM/YKLIJ
EE4S=EE4S+XNOM/YKLIJS
C
C --- start L3, L4
IF (L3RDF.OR.L4THF.OR.L34EPS) THEN
XTT1=H1/EKLIJ
XTT2=H2/EKLIJ
IF (LINFF) THEN
C XTT17=H1/DENSUM(INDX)
C XTT27=H2/DENSUM(INDX)
END IF
EM431=0.D0
EM432=0.D0
EM441=0.D0
EM451=0.D0
EM452=0.D0
EM461=0.D0
EM481=0.D0
EM482=0.D0
EE431=0.D0
EE432=0.D0
EE431S=0.D0
EE432S=0.D0
EE441=0.D0
EE441S=0.D0
C
DO II=1,NOCC
INDII=II
IF (INDII.NE.I) THEN
EII=ORBEN(II)
FF=F(I,II)
EDEN1=EKLIJ-EJ+EII
FDEN1=FF/EDEN1
INDIII=INDII
INDKK=INDK
INDLL=INDL
C CALL TURN4(INDII,INDK,INDL,IVJ,INDIII,INDKK,INDLL)
H3=HFIND(J,INDLL,INDIII,INDKK)
H5=HFIND(J,INDKK,INDIII,INDLL)
EM431=EM431 - (H5+H5-H3)*FDEN1
EM432=EM432 - (H3+H3-H5)*FDEN1
EM441=EM441+FF*FDEN1
EM451=EM451 - (H5+H5-H3)*FF/(EDEN1-FF*FF/EKLIJ)
EM452=EM452 - (H3+H3-H5)*FF/(EDEN1-FF*FF/EKLIJ)
EM461=EM461+FF*FF/(EDEN1*(1.D0-FF*FF/EDEN1/EKLIJ
$ ))
C
IF (L8) THEN
CALL FABINX(II,J,K,L,INDX)
C FDEN7=FF/DENSUM(INDX)
EM481=EM481 - (H5+H5-H3)*FDEN7
EM482=EM482 - (H3+H3-H5)*FDEN7
END IF
C Eps-Nes
IF (INDII.NE.INDJ) THEN
CALL EDC4N(II,J,K,L,H3,H5
$ ,EDEN1,E,ES)
ELSE
CALL EDC3N(J,K,L,EDEN1,E,ES)
END IF
EE431=EE431 - (H5+H5-H3)*FF/E
EE432=EE432 - (H3+H3-H5)*FF/E
EE441=EE441+FF*FF/E
EE431S=EE431S - (H5+H5-H3)*FF/ES
EE432S=EE432S - (H3+H3-H5)*FF/ES
EE441S=EE441S+FF*FF/ES
END IF
IF (INDII.NE.INDJ) THEN
FF=F(II,J)
EDEN2=EKLIJ-EI+EII
FDEN2=FF/EDEN2
H4=HFIND(I,INDK,INDII,INDL)
H6=HFIND(I,INDL,INDII,INDK)
EM432=EM432-(H4+H4-H6)*FDEN2
EM431=EM431-(H6+H6-H4)*FDEN2
EM441=EM441+FF*FDEN2
EM452=EM452-(H4+H4-H6)*FF/(EDEN2-FF*FF/EKLIJ)
EM451=EM451-(H6+H6-H4)*FF/(EDEN2-FF*FF/EKLIJ)
EM461=EM461+FF*FF/(EDEN2*(1.D0-FF*FF/EKLIJ/EDEN2
$ ))
IF (L8) THEN
CALL FABINX(I,II,K,L,INDX)
C FDEN7=FF/DENSUM(INDX)
EM482=EM482-(H4+H4-H6)*FDEN7
EM481=EM481-(H6+H6-H4)*FDEN7
END IF
C Eps-Nes
IF (INDII.NE.I) THEN
CALL EDC4N(I,II,K,L,H4,H6,EDEN2,E
$ ,ES)
ELSE
CALL EDC3N(I,K,L,EDEN2,E,ES)
END IF
EE432=EE432-(H4+H4-H6)*FF/E
EE431=EE431-(H6+H6-H4)*FF/E
EE441=EE441+FF*FF/E
EE432S=EE432S-(H4+H4-H6)*FF/ES
EE431S=EE431S-(H6+H6-H4)*FF/ES
EE441S=EE441S+FF*FF/ES
END IF
END DO
DO KK=NOV,NBAS
INDKK=KK
IF (INDKK.NE.INDK) THEN
EKK=ORBEN(KK)
FF=F(KK,K)
EDEN1=EKLIJ+EK-EKK
FDEN1=FF/EDEN1
H3=HFIND(I,INDKK,INDJ,INDL)
H5=HFIND(I,INDL,INDJ,INDKK)
EM431=EM431+(H5+H5-H3)*FDEN1
EM432=EM432+(H3+H3-H5)*FDEN1
EM441=EM441+FF*FDEN1
EM451=EM451+(H5+H5-H3)*FF/(EDEN1-FF*FF/EKLIJ)
EM452=EM452+(H3+H3-H5)*FF/(EDEN1-FF*FF/EKLIJ)
EM461=EM461+FF*FF/(EDEN1*(1.D0-FF*FF/EDEN1/EKLIJ
$ ))
IF (L8) THEN
CALL FABINX(I,J,KK,L,INDX)
C FDEN7=FF/DENSUM(INDX)
EM481=EM481+(H5+H5-H3)*FDEN7
EM482=EM482+(H3+H3-H5)*FDEN7
END IF
C Eps-Nes
IF (INDKK.NE.INDL) THEN
CALL EDC4N(I,J,KK,L,H3,H5,EDEN1,E
$ ,ES)
ELSE
CALL EDC3N(I,J,L,EDEN1,E,ES)
END IF
C WRITE(6,*) ' 3A+ ',I,INDJ,INDK,INDL,INDII,E,ES
EE431=EE431+(H5+H5-H3)*FF/E
EE432=EE432+(H3+H3-H5)*FF/E
EE441=EE441+FF*FF/E
EE431S=EE431S+(H5+H5-H3)*FF/ES
EE432S=EE432S+(H3+H3-H5)*FF/ES
EE441S=EE441S+FF*FF/ES
END IF
IF (INDKK.NE.INDL) THEN
FF=F(KK,L)
EDEN2=EKLIJ+EL-EKK
FDEN2=FF/EDEN2
H4=HFIND(I,INDK,INDJ,INDKK)
H6=HFIND(I,INDKK,INDJ,INDK)
EM432=EM432+(H4+H4-H6)*FDEN2
EM431=EM431+(H6+H6-H4)*FDEN2
EM441=EM441+FF*FDEN2
EM452=EM452+(H4+H4-H6)*FF/(EDEN2-FF*FF/EKLIJ)
EM451=EM451+(H6+H6-H4)*FF/(EDEN2-FF*FF/EKLIJ)
EM461=EM461+FF*FF/(EDEN2*(1.D0-FF*FF/EDEN2/EKLIJ
$ ))
IF (L8) THEN
CALL FABINX(I,J,K,KK,IDNX)
C FDEN7=FF/DENSUM(INDX)
EM482=EM482+(H4+H4-H6)*FDEN7
EM481=EM481+(H6+H6-H4)*FDEN7
END IF
C Eps-Nes
IF (INDKK.NE.INDK) THEN
CALL EDC4N(I,J,K,KK,H4,H6,EDEN2,E
$ ,ES)
ELSE
CALL EDC3N(I,J,K,EDEN2,E,ES)
END IF
C WRITE(6,*) ' 4A+ ',I,INDJ,INDK,INDL,INDII,E,ES
EE432=EE432+(H4+H4-H6)*FF/E
EE431=EE431+(H6+H6-H4)*FF/E
EE441=EE441+FF*FF/E
EE432S=EE432S+(H4+H4-H6)*FF/ES
EE431S=EE431S+(H6+H6-H4)*FF/ES
EE441S=EE441S+FF*FF/ES
END IF
END DO
EM43=EM43+XTT1*EM432+XTT2*EM431
EM44=EM44+EM441*XNOM/(EKLIJ*EKLIJ)
EM45=EM45+XTT1*EM452+XTT2*EM451
EM46=EM46+EM461*XNOM/(EKLIJ*EKLIJ)
IF (L8) EM48=EM48+XTT17*EM482+XTT27*EM481
EE43=EE43+1.D0/YKLIJ*(H1*EE432+H2*EE431)
EE44=EE44+EE441*XNOM/(YKLIJ*YKLIJ)
EE43S=EE43S+1.D0/YKLIJS*(H1*EE432S+H2*EE431S)
EE44S=EE44S+EE441S*XNOM/(YKLIJS*YKLIJS)
END IF
C --- end L3
END IF
END IF
END DO
END DO
END DO
END DO
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
C we have: no index --> 2nd order
C 3 --> thrd order F
C 4 --> 4th order F
C 5 --> 3th order F, Martin's diagrams -> infty
C 6 --> 4th order F, Martin's diagrams -> infty
C 7 --> geom. summation of pairs
C 8 --> 3rd order F + geom. summ. of pairs
EE1S=EE1
C
EM13 =EM13 *2.D0
IF (L8) EM18 =EM18 *2.D0
EE13 =EE13 *2.D0
EE13S=EE13S*2.D0
EM2 =EM2 *2.D0
EM23 =EM23 *2.D0
EM24 =EM24 *2.D0
EM27 =EM27 *2.D0
IF (L8) EM28 =EM28 *2.D0
EE2 =EE2 *2.D0
EE23 =EE23 *2.D0
EE24 =EE24 *2.D0
EE2S =EE2S *2.D0
EE23S=EE23S*2.D0
EE24S=EE24S*2.D0
C
EM3 =EM3 *2.D0
EM33 =EM33 *2.D0
EM34 =EM34 *2.D0
EM37 =EM37 *2.D0
IF (L8) EM38 =EM38 *2.D0
EE3 =EE3 *2.D0
EE33 =EE33 *2.D0
EE34 =EE34 *2.D0
EE3S =EE3S *2.D0
EE33S=EE33S*2.D0
EE34S=EE34S*2.D0
C
EM4 =EM4 *4.D0
EM43 =EM43 *2.D0
EM44 =EM44 *4.D0
EM47 =EM47 *4.D0
IF (L8) EM48 =EM48 *2.D0
EE4 =EE4 *4.D0
EE43 =EE43 *2.D0
EE44 =EE44 *4.D0
EE4S =EE4S *4.D0
EE43S=EE43S*2.D0
EE44S=EE44S*4.D0
C
EM15 =EM15 *2.D0
EM25 =EM25 *2.D0
EM35 =EM35 *2.D0
EM45 =EM45 *2.D0
EM26 =EM26 *2.D0
EM36 =EM36 *2.D0
EM46 =EM46 *4.D0
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
EMP2 =EM1 +EM2 +EM3 +EM4
EMP23=EM13+EM23+EM33+EM43
EMP24=EM14+EM24+EM34+EM44
EMP25=EM15+EM25+EM35+EM45-EMP23
EMP26=EM16+EM26+EM36+EM46-EMP24
EMP27=EM17+EM27+EM37+EM47-EMP2-EMP24
IF (L8) EMP28=EM18+EM28+EM38+EM48-EMP23
C
EEN2 =EE1 +EE2 +EE3 +EE4
EEN23=EE13+EE23+EE33+EE43
EEN24=EE14+EE24+EE34+EE44
C
EEN2S =EE1S +EE2S +EE3S +EE4S
EEN23S=EE13S+EE23S+EE33S+EE43S
EEN24S=EE14S+EE24S+EE34S+EE44S
C
WRITE(6,9913) '(MP2 LOC.) ',
- EMP2,EM1,EM2,EM3,EM4,EMP2
IF (L3RDF) WRITE(6,9913) '(3rd Fock) '
$ ,EMP23,EM13,EM23,EM33,EM43,EMP23
IF (L4THF) WRITE(6,9913) '(4th Fock) '
$ ,EMP24,EM14,EM24,EM34,EM44,EMP24
IF (LMART) THEN
WRITE(6,9913) '(3rd Fock, Martin) ',EMP25,EM15
$ -EM13,EM25-EM23,EM35-EM33,EM45-EM43,EMP25
WRITE(6,9913) '(4th Fock, Martin) ',EMP26,EM16
$ -EM14,EM26-EM24,EM36-EM34,EM46-EM44,EMP26
END IF
IF (LINFF) WRITE(6,9913) '(approx. geo. series in pairs on MP2)'
$ ,EMP27,EM17-EM1-EM14,EM27-EM2-EM24,EM37-EM3-EM34
$ ,EM47-EM4-EM44,EMP27
IF (L8) THEN
WRITE(6,9913) '(geom. series on 3rd Fock) ',
- EMP28,EM18,EM28,EM38,EM48,EMP28
END IF
WRITE(6,9913) '(EN2 Loc.)',EEN2,EE1,EE2,EE3,EE4,EEN2
IF (L34EPS) THEN
WRITE(6,9913) '(thrd order F on EN2)',EEN23,EE13,EE23
$ ,EE33,EE43,EEN23
WRITE(6,9913) '(4th order F on EN2)',EEN24,EE14,EE24
$ ,EE34,EE44,EEN24
END IF
WRITE(6,9913) '(S-A EN2 loc.)',
- EEN2S,EE1S,EE2S,EE3S,EE4S,EEN2S
IF (L34EPS) THEN
WRITE(6,9913) '(S-A EN2 + 3rd order F)',EEN23S,EE13S
$ ,EE23S,EE33S,EE43S,EEN23S
WRITE(6,9913) '(S-A EN2 + 4th order pairs)',EEN24S
$ ,EE14S,EE24S,EE34S,EE44S,EEN24S
END IF
9913 FORMAT(/,' CORRELATION ENERGY ',A15,':',F20.12,//,
- ' CONTRIBUTIONS II -> KK ',F20.12,/,
- ' CONTRIBUTIONS II -> KL ',F20.12,/,
- ' CONTRIBUTIONS IJ -> KK ',F20.12,/,
- ' CONTRIBUTIONS IJ -> KL ',F20.12,/,
- 28X,'----------------------',/,13X,'SUM',10X,F20.12,/)
C
C CONCENTRATED OUTPUT
C
WRITE(6,7715)
WRITE(6,7716) 'MP2 ',EMP2, EMP2
EEEE=EMP2+EMP23
IF (L3RDF) WRITE(6,7716) 'MP2/L3 ',EMP23,EEEE
EEEE=EEEE+EMP24
IF (L4THF)WRITE(6,7716) 'MP2/P4 ',EMP24,EEEE
IF (LMART) THEN
EEEE=EMP2+EMP25
WRITE(6,7716) 'MP2/M3 ',EMP25-EMP23,EEEE
EEEE=EEEE+EMP26
WRITE(6,7716) 'MP2/M4 ',EMP26-EMP24,EEEE
EEEE=EEEE-EMP25-EMP26+EMP24+EMP23
END IF
EEEE=EMP27+EEEE
IF (LINFF) WRITE(6,7716) 'MP2/G ',EMP27,EEEE
IF (L8) THEN
EEEE=EMP28+EEEE
WRITE(6,7716) 'MP2/G3 ',EMP28,EEEE
END IF
WRITE(6,*)
WRITE(6,7716) 'EN2 ',EEN2 , EEN2
EEEE=EEN2+EEN23
IF (L34EPS) THEN
WRITE(6,7716) 'EN2/F3 ',EEN23, EEEE
EEEE=EEEE+EEN24
WRITE(6,7716) 'EN2/P4 ',EEN24, EEEE
END IF
WRITE(6,*)
C
EEEE=EEN2S
WRITE(6,7716) 'EN2S ',EEN2S, EEEE
EEEE=EEEE+EEN23S
IF (L34EPS) THEN
WRITE(6,7716) 'EN2S/F3',EEN23S, EEEE
EEEE=EEEE+EEN24S
WRITE(6,7716) 'EN2S/P4',EEN24S, EEEE
END IF
WRITE(6,*)
7715 FORMAT(' METHOD',9X,'total',9X,'cumul. total',/,74('-'))
7716 FORMAT(A8,4F16.10)
C
RETURN
END
C
SUBROUTINE TRANSF(A,CI)
INCLUDE 'param.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
DIMENSION AJKLM(NBASM,NBASM)
DIMENSION A(NBASM,NBASM),CI(NBASM,NBASM)
C
DO J=1,NBAS
DO I=1,NBAS
AJKLM(I,J)=0.D0
END DO
END DO
C
DO J=1,NBAS
DO K=1,NBAS
BB=0.D0
DO L=1,NBAS
BB=BB+CI(L,J)*A(K,L)
END DO
AJKLM(K,J)=AJKLM(K,J)+BB
END DO
END DO
C
DO J=1,NBAS
DO I=1,NBAS
A(I,J)=0.D0
END DO
END DO
C
DO I=1,NBAS
C only J=I..NBAS, copying is cheaper
DO J=I,NBAS
SSS=0.D0
DO K=1,NBAS
SSS=SSS+AJKLM(K,J)*CI(K,I)
END DO
A(I,J)=A(I,J)+SSS
END DO
END DO
C
C copy A(I,J)=A(J,I)
C
DO I=1,NBAS
DO J=I+1,NBAS
A(J,I)=A(I,J)
END DO
END DO
C
RETURN
END
C
SUBROUTINE PFUNC(CI,NBP)
INCLUDE 'param.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
DIMENSION CI(NBASM,NBASM)
DO IVEC=1,NBP
WRITE(6,*) ' function index ',IVEC
WRITE(6,'(5X,4E15.8)') (CI(J,IVEC),J=1,NBAS)
END DO
RETURN
END
C
FUNCTION HFIND(I1,J1,K1,L1)
INCLUDE 'param.h'
PARAMETER (NBULL=43000000)
COMMON /TWOI/ H0(NBULL),IH0(4,NBULL),NUMINT
C.. INCLUDE 'nbull.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
LOGICAL LFOU,LUP
C
CT IF (I1.EQ.0) STOP 'I1.EQ.0'
CT IF (J1.EQ.0) STOP 'J1.EQ.0'
CT IF (K1.EQ.0) STOP 'K1.EQ.0'
CT IF (L1.EQ.0) STOP 'L1.EQ.0'
C I1,K1 besetzt, J1,L1 leer
C ordne Paare (I1,J1) und (K1,L1)
I=I1
J=J1
K=K1
L=L1
IF (I.GT.J) THEN
IDUM=I
I=J
J=IDUM
END IF
IF (K.GT.L) THEN
IDUM=K
K=L
L=IDUM
END IF
IF (I.GT.K) THEN
IDUM=I
I=K
K=IDUM
IDUM=J
J=L
L=IDUM
ELSE IF (I.EQ.K) THEN
IF (J.GT.L) THEN
IDUM=J
J=L
L=IDUM
END IF
END IF
C
C Integrals are lexically ordered according to IH0
C
IND=(NUMINT+1)/2
INC=(IND+1)/2+1
LFOU=.FALSE.
100 CONTINUE
C WRITE(6,*) '395: ',NUMINT,IND,INC,(IH0(JJJ,IND),JJJ=1,4),I,J,K,L
IF (IH0(1,IND).LT.I) THEN
LUP=.TRUE.
ELSE IF (IH0(1,IND).GT.I) THEN
LUP=.FALSE.
ELSE IF (IH0(1,IND).EQ.I) THEN
IF (IH0(2,IND).LT.J) THEN
LUP=.TRUE.
ELSE IF (IH0(2,IND).GT.J) THEN
LUP=.FALSE.
ELSE IF (IH0(2,IND).EQ.J) THEN
IF (IH0(3,IND).LT.K) THEN
LUP=.TRUE.
ELSE IF (IH0(3,IND).GT.K) THEN
LUP=.FALSE.
ELSE IF (IH0(3,IND).EQ.K) THEN
IF (IH0(4,IND).LT.L) THEN
LUP=.TRUE.
ELSE IF (IH0(4,IND).GT.L) THEN
LUP=.FALSE.
ELSE IF (IH0(4,IND).EQ.L) THEN
HFIND=H0(IND)
RETURN
END IF
END IF
END IF
END IF
C
IF (LUP) THEN
IND=MIN(NUMINT,IND+INC)
INC=(INC+1)/2
ELSE
IND=MAX(1,IND-INC)
INC=(INC+1)/2
END IF
IF (IND.LT.1) IND=1
IF (INC.EQ.1) GO TO 200
GO TO 100
C
200 CONTINUE
C perhaps IND
INDX=IND
IF (IH0(1,INDX).EQ.I) THEN
IF (IH0(2,INDX).EQ.J) THEN
IF (IH0(3,INDX).EQ.K) THEN
IF (IH0(4,INDX).EQ.L) THEN
HFIND=H0(INDX)
RETURN
END IF
END IF
END IF
END IF
C perhaps IND-1
IF (IND.GT.1) THEN
INDX=IND-1
IF (IH0(1,INDX).EQ.I) THEN
IF (IH0(2,INDX).EQ.J) THEN
IF (IH0(3,INDX).EQ.K) THEN
IF (IH0(4,INDX).EQ.L) THEN
HFIND=H0(INDX)
RETURN
END IF
END IF
END IF
END IF
END IF
C or IND+1
IF (IND.LT.NUMINT) THEN
INDX=IND+1
IF (IH0(1,INDX).EQ.I) THEN
IF (IH0(2,INDX).EQ.J) THEN
IF (IH0(3,INDX).EQ.K) THEN
IF (IH0(4,INDX).EQ.L) THEN
HFIND=H0(INDX)
RETURN
END IF
END IF
END IF
END IF
END IF
C
C no, not found
C
HFIND=0.D0
C WRITE(6,*) 'WE HAVE NO INTEGRAL ...',I,J,K,L
RETURN
C STOP
END
C
FUNCTION IFIND(I1,J1,K1,L1)
INCLUDE 'param.h'
PARAMETER (NBULL=43000000)
COMMON /TWOI/ H0(NBULL),IH0(4,NBULL),NUMINT
C.. INCLUDE 'nbull.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
LOGICAL LUP
C
CT IF (I1.EQ.0) STOP 'I1.EQ.0'
CT IF (J1.EQ.0) STOP 'J1.EQ.0'
CT IF (K1.EQ.0) STOP 'K1.EQ.0'
CT IF (L1.EQ.0) STOP 'L1.EQ.0'
C I1,K1 besetzt, J1,L1 leer
C ordne Paare (I1,J1) und (K1,L1)
I=I1
J=J1
K=K1
L=L1
IF (I.GT.K) THEN
IDUM=I
I=K
K=IDUM
IDUM=J
J=L
L=IDUM
ELSE IF (I.EQ.K) THEN
IF (J.GT.L) THEN
IDUM=J
J=L
L=IDUM
END IF
END IF
C
C Integrals are lexically ordered according to IH0
C
IND=(NUMINT+1)/2
INC=(IND+1)/2
100 CONTINUE
C WRITE(6,*) '395: ',NUMINT,IND,INC,(IH0(JJJ,IND),JJJ=1,4),I,J,K,L
IF (IH0(1,IND).LT.I) THEN
LUP=.TRUE.
ELSE IF (IH0(1,IND).GT.I) THEN
LUP=.FALSE.
ELSE IF (IH0(1,IND).EQ.I) THEN
IF (IH0(2,IND).LT.J) THEN
LUP=.TRUE.
ELSE IF (IH0(2,IND).GT.J) THEN
LUP=.FALSE.
ELSE IF (IH0(2,IND).EQ.J) THEN
IF (IH0(3,IND).LT.K) THEN
LUP=.TRUE.
ELSE IF (IH0(3,IND).GT.K) THEN
LUP=.FALSE.
ELSE IF (IH0(3,IND).EQ.K) THEN
IF (IH0(4,IND).LT.L) THEN
LUP=.TRUE.
ELSE IF (IH0(4,IND).GT.L) THEN
LUP=.FALSE.
ELSE IF (IH0(4,IND).EQ.L) THEN
IFIND=IND
RETURN
END IF
END IF
END IF
END IF
C
IF (LUP) THEN
IND=IND+INC
INC=(INC+1)/2
ELSE
IND=IND-INC
INC=(INC+1)/2
END IF
IF (IND.LT.1) IND=1
IF (INC.EQ.1) GO TO 200
GO TO 100
C
200 CONTINUE
C perhaps IND
INDX=IND
IF (IH0(1,INDX).EQ.I) THEN
IF (IH0(2,INDX).EQ.J) THEN
IF (IH0(3,INDX).EQ.K) THEN
IF (IH0(4,INDX).EQ.L) THEN
IFIND=INDX
RETURN
END IF
END IF
END IF
END IF
C perhaps IND-1
IF (IND.GT.1) THEN
INDX=IND-1
IF (IH0(1,INDX).EQ.I) THEN
IF (IH0(2,INDX).EQ.J) THEN
IF (IH0(3,INDX).EQ.K) THEN
IF (IH0(4,INDX).EQ.L) THEN
IFIND=INDX
RETURN
END IF
END IF
END IF
END IF
END IF
C or IND+1
IF (IND.LT.NUMINT) THEN
INDX=IND+1
IF (IH0(1,INDX).EQ.I) THEN
IF (IH0(2,INDX).EQ.J) THEN
IF (IH0(3,INDX).EQ.K) THEN
IF (IH0(4,INDX).EQ.L) THEN
IFIND=INDX
RETURN
END IF
END IF
END IF
END IF
END IF
C
C no, not found
C
IFIND=NUMINT+1
C WRITE(6,*) 'WE HAVE NO INDEX ...',I,J,K,L
RETURN
C STOP
END
C
SUBROUTINE INTSRT(IH0,H0,NUMINT)
INCLUDE 'param.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
DIMENSION H0(NUMINT),IH0(4,NUMINT)
ISRT=0
C
C check lexical order of integrals
C
I=1
11 CONTINUE
IF (IH0(1,I).GT.IH0(1,I+1)) THEN
C WRITE(6,'(2I7,3I4,I8,3I4)') I,(IH0(J,I),J=1,4),(IH0(K,I+1),K=1,4)
ISRT=1
ELSE IF (IH0(1,I).EQ.IH0(1,I+1)) THEN
IF (IH0(2,I).GT.IH0(2,I+1)) THEN
ISRT=1
ELSE IF (IH0(2,I).EQ.IH0(2,I+1)) THEN
IF (IH0(3,I).GT.IH0(3,I+1)) THEN
ISRT=1
ELSE IF (IH0(3,I).EQ.IH0(3,I+1)) THEN
IF (IH0(4,I).GT.IH0(4,I+1)) THEN
ISRT=1
ELSE IF (IH0(4,I).EQ.IH0(4,I+1)) THEN
WRITE(6,'(I7,3I4,I8,3I4)') (IH0(J,I),J=1,4),(IH0(K,I+1),K=1,4)
STOP 'IDENTICAL INDICES FOR DIFFERENT INTEGRALS ENCOUNTERED'
END IF
END IF
END IF
END IF
IF (ISRT.EQ.0) THEN
I=I+1
IF (I.LT.NUMINT-1) GO TO 11
END IF
C
CALL TIMING('CHCK')
IF (ISRT.EQ.1) THEN
C
WRITE(6,*) ' INTEGRALS HAVE TO BE SORTED '
WRITE(6,*) ' WE HAVE ',NUMINT,' INTEGRALS in core'
CALL LEXSRT(IH0,H0,NUMINT)
CALL TIMING('SORT')
ELSE
WRITE(6,*) ' MP2 INTEGRALS ARE WELL-ORDERED'
END IF
C
RETURN
END
C
SUBROUTINE READ53(IUNIT4)
INCLUDE 'param.h'
PARAMETER (NBULL=43000000)
COMMON /TWOI/ H0(NBULL),IH0(4,NBULL),NUMINT
C.. INCLUDE 'nbull.h'
COMMON /FLOW/ THRESH,ICUT,NFROZ,LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX
$ ,L3RDF,L4THF,LINFF,LCUT,LTHRE,LMART,L34EPS,LMP3
LOGICAL LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX,L3RDF,L4THF,LINFF,LCUT
$ ,LTHRE,LMART,L34EPS,LMP3
C.. INCLUDE 'flow.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /WBUF/ HVAL(IWBULL),IIND(IWBULL),JIND(IWBULL),
- KIND(IWBULL),LIND(IWBULL),IPOS,IUNIT,IBLE
C.. INCLUDE 'common_wbuf.h'
COMMON /CNTERS/ NMP2,NMP3,NOOOO,NOVOV,NOOVV,NVVVV,NCOU,NEXC,IREJ
COMMON /BIBUF/HOOOO(IWBULL),HOOVV(IWBULL),HVVVV(IWBULL),IOOOO(4
$ ,IWBULL),IOOVV(4,IWBULL),IVVVV(4,IWBULL),NOOOOD,NOOVVD,NVVVVD
$ ,IUOOOO,IUOOVV,IUVVVV
C.. INCLUDE 'counters.h'
LOGICAL LF,LI,LJ,LA,LB,LII,LJJ,LAA,LBB
C
C read transformed two-electron integrals, and select MP2 (MP3) integrals
C
IREJ=0
NCOU=0
NEXC=0
NMP2=0
NOVOV=0
NMP3=0
NOOOO=0
NOOVV=0
NVVVV=0
NOOOOD=0
NOOVVD=0
NVVVVD=0
IF (LMP3) THEN
IUOOOO=55
IUOOVV=56
IUVVVV=57
OPEN(UNIT=IUOOOO,FILE='OOOO',STATUS='UNKNOWN',FORM='UNFORMATTED')
OPEN(UNIT=IUOOVV,FILE='OOVV',STATUS='UNKNOWN',FORM='UNFORMATTED')
OPEN(UNIT=IUVVVV,FILE='VVVV',STATUS='UNKNOWN',FORM='UNFORMATTED')
WRITE(6,*) ' OPENED UNIT 55 as OOOO'
WRITE(6,*) ' OPENED UNIT 56 as OOVV'
WRITE(6,*) ' OPENED UNIT 57 as VVVV'
END IF
LF=.FALSE.
NUMINT=0
IF (LUNFOR) THEN
IPOS=1
IMODE=2
CALL WOPEN(IUNIT4,IMODE)
11 CONTINUE
CALL WREAD(LF)
NUMINT=NUMINT+IPOS
DO I=1,IPOS
INDI=IIND(I)
INDJ=JIND(I)
INDK=KIND(I)
INDL=LIND(I)
HHH=HVAL(I)
CALL HRANGE(INDI,INDJ,INDK,INDL,HHH)
END DO
IF (.NOT.LF) GO TO 11
12 CONTINUE
CALL WCLOS(IMODE)
ELSE
OPEN(UNIT=IUNIT4,FILE='FILE53',FORM='FORMATTED',STATUS='OLD')
111 CONTINUE
READ(IUNIT4,*,IOSTAT=KK) INDI,INDJ,INDK,INDL,HHH
IF (KK.NE.0) GO TO 112
NUMINT=NUMINT+1
CALL HRANGE(INDI,INDJ,INDK,INDL,HHH)
IF (.NOT.LF) GO TO 111
112 CONTINUE
END IF
C
C we have to close the files OOOO, OOVV, VVVV
IF (LMP3) THEN
IF (NOOOOD+1.EQ.IWBULL) THEN
WRITE(IUOOOO) ((IOOOO(JDUM,IDUM),JDUM=1,4),HOOOO(IDUM),IDUM=1
$ ,IWBULL)
NOOOOD=0
END IF
HOOOO(IWBULL)=-10000.D0
IOOOO(1,IWBULL)=NOOOOD
WRITE(IUOOOO) ((IOOOO(JDUM,IDUM),JDUM=1,4),HOOOO(IDUM),IDUM=1
$ ,IWBULL)
CLOSE(IUOOOO)
C
IF (NOOVVD+1.EQ.IWBULL) THEN
WRITE(IUOOVV) ((IOOVV(JDUM,IDUM),JDUM=1,4),HOOVV(IDUM),IDUM=1
$ ,IWBULL)
NOOVVD=0
END IF
HOOVV(IWBULL)=-10000.D0
IOOVV(1,IWBULL)=NOOVVD
WRITE(IUOOVV) ((IOOVV(JDUM,IDUM),JDUM=1,4),HOOVV(IDUM),IDUM=1
$ ,IWBULL)
CLOSE(IUOOVV)
C
IF(NVVVVD+1.EQ.IWBULL) THEN
WRITE(IUVVVV) ((IVVVV(JDUM,IDUM),JDUM=1,4),HVVVV(IDUM),IDUM=1
$ ,IWBULL)
NVVVVD=0
END IF
HVVVV(IWBULL)=-10000.D0
IVVVV(1,IWBULL)=NVVVVD
WRITE(IUVVVV) ((IVVVV(JDUM,IDUM),JDUM=1,4),HVVVV(IDUM),IDUM=1
$ ,IWBULL)
CLOSE(IUVVVV)
END IF
C
C remap the indices
C
IF (LFROZ) THEN
IF (NFROZ.NE.0) THEN
DO I=NFROZ+1,NBAS
II=I-NFROZ
ORBEN(II)=ORBEN(I)
DO J=NFROZ+1,NBAS
JJ=J-NFROZ
F(II,JJ)=F(I,J)
HCOU(II,JJ)=HCOU(I,J)
HEXC(II,JJ)=HEXC(I,J)
END DO
END DO
C
C these have been remapped already in HRANGE
C
c$$$ DO I=1,NMP2+NMP3
c$$$ DO J=1,4
c$$$ II=IH0(J,I)
c$$$ IFU=II
c$$$ IFU=IFU-NFROZ
c$$$ II=IFU
c$$$ IH0(J,I)=II
c$$$ END DO
c$$$ END DO
NOCC=NOCC-NFROZ
NBAS=NBAS-NFROZ
WRITE(6,*) ' NEW NBAS, NOCC: ',NBAS,NOCC
END IF
END IF
IF (LMP2PR) THEN
DO I=1,NMP2+NMP3
WRITE(6,*) ' SEL: ',(IH0(J,I),J=1,4),H0(I)
END DO
END IF
WRITE(6,*)
WRITE(6,*) ' READ ',NUMINT,' BI- + MONOELECTRONIC INTEGRALS'
IF (LMP3) THEN
WRITE(6,9803) NMP2,NMP3,NOOOO,NOVOV,NOOVV,NVVVV
9803 FORMAT(' KEPT ',I10,' MP2 INTEGRALS',/,' KEPT ',I10
$ ,' MP3 INTEGRALS',//,' KEPT ',I10
$ ,' OOOO INTEGRALS',/,' KEPT ',I10
$ ,' OVOV INTEGRALS',/,' KEPT ',I10
$ ,' OOVV INTEGRALS',/,' KEPT ',I10
$ ,' VVVV INTEGRALS',/)
ELSE
WRITE(6,*) ' FOUND ',NMP2+IREJ,' MP2 INTEGRALS'
WRITE(6,*) ' KEPT ',NOVOV, ' OVOV INTEGRALS'
WRITE(6,*) ' REJECTED ',IREJ, ' MP2 INTEGRALS'
END IF
WRITE(6,*) ' FOUND ',NCOU ,' COULOMB INTEGRALS'
WRITE(6,*) ' FOUND ',NEXC ,' EXCHANGE INTEGRALS'
WRITE(6,*)
C
CALL TIMING('RD2I')
C
C now all integrals are in-core
C
NUMINT=NMP2
C
C sort lexically (heapsort)
C
CALL INTSRT(IH0,H0,NUMINT)
CD WRITE(6,*) 'NUMINT :',NUMINT
CALL TIMING('SRT ')
C
C get nonzero MP2 Integrals
C
NN0=0
DO I=1,NUMINT
IF (H0(I).NE.0.D0) NN0=NN0+1
END DO
WRITE(6,*) ' NON-ZERO INTEGRALS: ',NN0
C
RETURN
END
C
SUBROUTINE TOTCAL(E1,E2,EN)
INCLUDE 'param.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
C
ETWO=0.D0
EONE=0.D0
DO I=1,NOCC
DO J=1,NOCC
ETWO=ETWO+2.D0*HCOU(I,J)-HEXC(I,J)
END DO
EONE=EONE+ORBEN(I)
END DO
EONE=EONE*2.D0
ETOT=EN+EONE-ETWO
EONE=ETOT-EN-ETWO
C
WRITE(6,*)
WRITE(6,*) ' RECALCULATED TOTAL ENERGY ',ETOT
WRITE(6,*) ' RECALCULATED ONE-ELECTRON ENERGY ',EONE
WRITE(6,*) ' RECALCULATED TWO-ELECTRON ENERGY ',ETWO
WRITE(6,*) ' NUCLEAR REPULSION ',EN
WRITE(6,*)
WRITE(6,*) ' READ TOTAL ENERGY ',E1+E2+EN
WRITE(6,*) ' READ ONE-ELECTRON ENERGY ',E1
WRITE(6,*) ' READ TWO-ELECTRON ENERGY ',E2
WRITE(6,*)
CALL TIMING('ETOT')
C
RETURN
END
C
SUBROUTINE DC2(I,J,K,L,HHH,EDEN,E,ES)
INCLUDE 'param.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
C ii -> kl
EDE=EDEN-HCOU(I,I)-HCOU(K,L)
- +2.D0*HCOU(I,K)+2.D0*HCOU(I,L)
- -HEXC(I,K)-HEXC(I,L)
E=HHH*HHH/EDE
ES=HHH*HHH/(EDE-HEXC(K,L))
RETURN
END
SUBROUTINE EDC3(I,J,K,L,HHH,EDEN,E,ES)
INCLUDE 'param.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
C ij -> kk
EDE=EDEN-HCOU(I,J)-HCOU(K,K)
- +2.D0*HCOU(I,K)+2.D0*HCOU(J,K)
- -HEXC(I,K)-HEXC(J,K)
E=HHH*HHH/EDE
ES=HHH*HHH/(EDE-HEXC(I,J))
RETURN
END
C
SUBROUTINE EDC4(I,J,K,L,H1,H2,EDEN,E,ES)
INCLUDE 'param.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
C
CHI=EDEN-HCOU(I,J)-HCOU(K,L)+HCOU(I,K)+HCOU(I,L)
- +HCOU(J,K)+HCOU(J,L)
C
C alternative, direct calculation of ES
C
H11=H1*H1
H22=H2*H2
H12=H1*H2
C
XNOM2=(H11+H22-H12)
XNRM=2.D0*SQRT(H11+H22-H12)
XKIJ=HEXC(I,J)
XKRS=HEXC(K,L)
XKIR=HEXC(I,K)
XKJR=HEXC(J,K)
XKIS=HEXC(I,L)
XKJS=HEXC(J,L)
C
XDEN2=-CHI*XNOM2*4.D0+2.D0*(XKIJ+XKRS)*(4.D0*H12-H11-H22)+2.D0
$ *(XKIS+XKJR)*(H11+4.D0*H22-4.D0*H12)+2.D0*(XKJS+XKIR)*(H22+4
$ .D0*H11-4.D0*H12)
XDEN2=-XDEN2/XNRM/XNRM
C
ES=XNOM2/XDEN2
E1=CHI+HEXC(I,J)+HEXC(K,L)-HEXC(I,K)-HEXC(I,L)
$ -HEXC(J,K)-HEXC(J,L)
E2=CHI-HEXC(I,K)-HEXC(J,L)
E3=CHI-HEXC(I,L)-HEXC(J,K)
E=(H11+H22-H12-H12)/E1+H11/E2+H22/E3
E=E*.5D0
C
RETURN
END
C
SUBROUTINE EDC2N(I,K,L,EDEN,E,ES)
INCLUDE 'param.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
C ii -> kl
C only Epstein-Nesbet denominator
E=EDEN-HCOU(I,I)-HCOU(K,L)
- +2.D0*HCOU(I,K)+2.D0*HCOU(I,L)
- -HEXC(I,K)-HEXC(I,L)
C the spin-adapted one
ES=E-HEXC(K,L)
C WRITE(6,*) ' EDC2N ',EDEN,E,ES
RETURN
END
SUBROUTINE EDC3N(I,J,K,EDEN,E,ES)
INCLUDE 'param.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
C ij -> kk
C only Epstein-Nesbet denominator
E=EDEN-HCOU(I,J)-HCOU(K,K)
- +2.D0*HCOU(I,K)+2.D0*HCOU(J,K)
- -HEXC(I,K)-HEXC(J,K)
ES=E-HEXC(I,J)
C WRITE(6,*) ' EDC3N ',EDEN,E,ES
RETURN
END
C
SUBROUTINE EDC4N(I,J,K,L,H1,H2,EDEN,E,ES)
INCLUDE 'param.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
C only Epstein-Nesbet denominator
C WRITE(6,*) ' EDC4N: ',I,J,K,L,H1,H2
CHI=EDEN-HCOU(I,J)-HCOU(K,L)+HCOU(I,K)+HCOU(I,L)
- +HCOU(J,K)+HCOU(J,L)
H11=H1*H1
H22=H2*H2
H12=H1*H2
C
IF (H1.NE.0.D0.OR.H2.NE.0.D0) THEN
XNOM2=(H11+H22-H12)
XNRM=2.D0*SQRT(H11+H22-H12)
XKIJ=HEXC(I,J)
XKRS=HEXC(K,L)
XKIR=HEXC(I,K)
XKJR=HEXC(J,K)
XKIS=HEXC(I,L)
XKJS=HEXC(J,L)
C
XDEN2=-CHI*XNOM2*4.D0+2.D0*(XKIJ+XKRS)*(4.D0*H12-H11-H22)
- +2.D0*(XKIS+XKJR)*(H11+4.D0*H22-4.D0*H12)
- +2.D0*(XKJS+XKIR)*(H22+4.D0*H11-4.D0*H12)
XDEN2=-XDEN2/XNRM/XNRM
C
ES=XDEN2
ELSE
ES=1.D0
END IF
C
E1=CHI+HEXC(I,J)+HEXC(K,L)-HEXC(I,K)-HEXC(I,L)
- -HEXC(J,K)-HEXC(J,L)
E2=CHI-HEXC(I,K)-HEXC(J,L)
E3=CHI-HEXC(I,L)-HEXC(J,K)
E=(H11+H22-H12-H12)/E1+H11/E2+H22/E3
XNN=H11+H22-H12
IF (E.NE.0.D0) THEN
E=2.D0*XNN/E
ELSE
E=1.D0
END IF
C WRITE(6,*) ' EDC4N ',I,J,K,L,IVJ,IVK,IVL,H1,H2,EDEN,E,ES
C
RETURN
END
C
SUBROUTINE RDINP
INCLUDE 'param.h'
COMMON /FLOW/ THRESH,ICUT,NFROZ,LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX
$ ,L3RDF,L4THF,LINFF,LCUT,LTHRE,LMART,L34EPS,LMP3
LOGICAL LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX,L3RDF,L4THF,LINFF,LCUT
$ ,LTHRE,LMART,L34EPS,LMP3
C.. INCLUDE 'flow.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
CHARACTER*4 KEYW(2),STR3
CHARACTER*6 KEYOPT(15),STR6
CHARACTER*80 LINE
DATA KEYW /'*EPS','*END'/
DATA KEYOPT /'MP2PRI','FOCKPR','COUEXC','XXXXXX','3RD F ',
- '4TH F ','INF F ','FORMAT','THRESH','FREEZE',
- 'MARTIN','MP3 ','EPS34 ','XXXXXX','XXXXXX'/
C
C DEFAULTS
C
THRESH=1.D-10
NFROZ =0
LMP3=.FALSE.
LFROZ =.FALSE.
LUNFOR=.TRUE.
LMP2PR=.FALSE.
LPRIFO=.FALSE.
LPRCEX=.FALSE.
L3RDF =.FALSE.
L4THF =.FALSE.
LINFF =.FALSE.
LCUT =.FALSE.
LTHRE =.FALSE.
LMART =.FALSE.
L34EPS=.FALSE.
C structure of input:
C keyword for each program
IOINP=83
OPEN(UNIT=IOINP,FILE='INPUT.EPS',ERR=2217,FORM='FORMATTED',
- STATUS='OLD')
C first, look for keyword '*EPS'
1 CONTINUE
READ(IOINP,'(A4)',END=2,ERR=921) STR3
IF (STR3.EQ.KEYW(1)) THEN
C look for Keyoptions
11 CONTINUE
READ(IOINP,'(A6)',END=920,ERR=921) STR6
C this is the keyword *END
IF (STR6(1:4).EQ.KEYW(2)) RETURN
IF (STR6.EQ.KEYOPT(1)) THEN
C print MP2 integrals
LMP2PR=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(2)) THEN
C print Fock matrix
LPRIFO=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(3)) THEN
C print Coulomb and Exchange integrals
LPRCEX=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(5)) THEN
C 3rd-order corrections for off-diagonal Fock matrix elements
L3RDF=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(6)) THEN
C 4th-order corrections for off-diagonal Fock matrix elements
L4THF=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(7)) THEN
C infinte summation of 4th order diagrams
LINFF=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(8)) THEN
C read integrals formatted
LUNFOR=.FALSE.
ELSE IF (STR6.EQ.KEYOPT(9)) THEN
C define integrals threshold
READ(IOINP,*) IDUM
THRESH=10.D0**(-DBLE(IDUM))
ELSE IF (STR6.EQ.KEYOPT(10)) THEN
C freeze occupied orbitals
LFROZ=.TRUE.
NFROZ=0
READ(5,*) NFROZ
IF (NFROZ.GE.NOCC) STOP ' TOO MANY ORBITALS FROZEN '
ELSE IF (STR6.EQ.KEYOPT(11)) THEN
C Martin Albrecht's summation of diagrams
LMART=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(12)) THEN
C MP3
LMP3=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(13)) THEN
C order 3 and 4 in Epstein-Nesbet
L34EPS=.TRUE.
ELSE
WRITE(6,*) ' POSSIBLE OPTIONS IN THE BLOCK *EPS ... *END ARE:'
WRITE(6,*)
WRITE(6,*) ' ',KEYOPT(1)
WRITE(6,*) ' ',KEYOPT(2)
WRITE(6,*) ' ',KEYOPT(3)
WRITE(6,*) ' ',KEYOPT(5)
WRITE(6,*) ' ',KEYOPT(6)
WRITE(6,*) ' ',KEYOPT(7)
WRITE(6,*) ' ',KEYOPT(8)
WRITE(6,*) ' ',KEYOPT(9)
WRITE(6,*) ' ',KEYOPT(10)
WRITE(6,*) ' ',KEYOPT(11)
WRITE(6,*) ' ',KEYOPT(12)
WRITE(6,*) ' ',KEYOPT(13)
STOP ' CHOOSE CORRECT OPTION '
END IF
GO TO 11
END IF
GO TO 1
C
921 CONTINUE
CLOSE(IOINP)
STOP ' ERROR IN INPUT'
2 CONTINUE
CLOSE(IOINP)
STOP ' NO KEYWORD *EPS FOUND, USING DEFAULTS '
920 CONTINUE
CLOSE(IOINP)
STOP ' YOU SHOULD TERMINATE THE BLOCK EPS BY <*END> '
2217 CONTINUE
WRITE(6,*) ' NO FILE , USING THE DEFAULT VALUES'
RETURN
END
C
SUBROUTINE HRANGE(INDI,INDJ,INDK,INDL,HHH)
INCLUDE 'param.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
PARAMETER (NBULL=43000000)
COMMON /TWOI/ H0(NBULL),IH0(4,NBULL),NUMINT
C.. INCLUDE 'nbull.h'
COMMON /FLOW/ THRESH,ICUT,NFROZ,LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX
$ ,L3RDF,L4THF,LINFF,LCUT,LTHRE,LMART,L34EPS,LMP3
LOGICAL LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX,L3RDF,L4THF,LINFF,LCUT
$ ,LTHRE,LMART,L34EPS,LMP3
C.. INCLUDE 'flow.h'
COMMON /CNTERS/ NMP2,NMP3,NOOOO,NOVOV,NOOVV,NVVVV,NCOU,NEXC,IREJ
COMMON /BIBUF/HOOOO(IWBULL),HOOVV(IWBULL),HVVVV(IWBULL),IOOOO(4
$ ,IWBULL),IOOVV(4,IWBULL),IVVVV(4,IWBULL),NOOOOD,NOOVVD,NVVVVD
$ ,IUOOOO,IUOOVV,IUVVVV
C.. INCLUDE 'counters.h'
LOGICAL LF,LHOLD
C occupied orbitals
LOGICAL LIO,LJO,LKO,LLO
C core orbitals
LOGICAL LIC,LJC,LKC,LLC
C
IF (ABS(HHH).LE.THRESH) HHH=0.D0
IF (INDK.NE.0) THEN
IF (LFROZ) THEN
LIC=INDI.LE.NFROZ
LJC=INDJ.LE.NFROZ
LKC=INDK.LE.NFROZ
LLC=INDL.LE.NFROZ
ELSE
LIC=.FALSE.
LJC=.FALSE.
LKC=.FALSE.
LLC=.FALSE.
END IF
C
C finding Coulomb and Exchange integrals
C
IF (INDI.EQ.INDJ.AND.INDK.EQ.INDL) THEN
NCOU=NCOU+1
HCOU(INDI,INDK)=HHH
HCOU(INDK,INDI)=HHH
IF (LPRCEX) WRITE(6,*) ' COULOMB ',NCOU,INDI,INDJ,INDK,INDL
$ ,HHH
END IF
IF (INDI.EQ.INDK.AND.INDJ.EQ.INDL) THEN
NEXC=NEXC+1
HEXC(INDI,INDJ)=HHH
HEXC(INDJ,INDI)=HHH
IF (LPRCEX) WRITE(6,*) ' EXCHANGE ',NEXC,INDI,INDJ,INDK,INDL
$ ,HHH
END IF
C
C classifying bielectronic integral
C
C if there is a core orbital we return
IF (LIC) RETURN
IF (LJC) RETURN
IF (LKC) RETURN
IF (LLC) RETURN
C
C otherwise we substract from all indices the number of core orbitals
INDI=INDI-NFROZ
INDJ=INDJ-NFROZ
INDK=INDK-NFROZ
INDL=INDL-NFROZ
C
LIO=INDI.LE.(NOCC-NFROZ)
LJO=INDJ.LE.(NOCC-NFROZ)
LKO=INDK.LE.(NOCC-NFROZ)
LLO=INDL.LE.(NOCC-NFROZ)
C
LHOLD=.FALSE.
C we don't store integrals with core orbitals
IF (LMP3) THEN
C (oo|vv), (oo|vv), (vv|vv)
IF ((LIO.EQV.LJO).AND.(LKO.EQV.LLO)) THEN
NMP3=NMP3+1
CALL ORDIND(INDI,INDJ,INDK,INDL,I,J,K,L)
C
C (aa|aa) 1/8
C (aa|ab) 1/2
C (aa|bb) 1/4
C (ab|ab) 1/2
C (aa|bc) 1/2
C (ab|bc) 1
C (ab|cd) 1
C see extract.r
FACTOR=1.D0
IF (I.EQ.J) FACTOR=FACTOR*0.5D0
IF (K.EQ.L) FACTOR=FACTOR*0.5D0
IF ((I.EQ.K).AND.(I.EQ.J).AND.(K.EQ.L)) FACTOR=FACTOR*0.5D0
IF ((I.EQ.K).AND.(J.EQ.L).AND..NOT.(I.EQ.J)) FACTOR=0.50D0
HHH=HHH*FACTOR
C
IF (LIO.NEQV.LKO) THEN
NOOVV=NOOVV+1
IF (NOOVVD.EQ.IWBULL) THEN
WRITE(IUOOVV) ((IOOVV(JDUM,IDUM),JDUM=1,4),HOOVV(IDUM),IDUM=1
$ ,IWBULL)
NOOVVD=0
END IF
NOOVVD=NOOVVD+1
IOOVV(1,NOOVVD)=I
IOOVV(2,NOOVVD)=J
IOOVV(3,NOOVVD)=K
IOOVV(4,NOOVVD)=L
HOOVV(NOOVVD)=HHH
ELSE
IF (LIO) THEN
NOOOO=NOOOO+1
IF (NOOOOD.EQ.IWBULL) THEN
WRITE(IUOOOO) ((IOOOO(JDUM,IDUM),JDUM=1,4),HOOOO(IDUM),IDUM
$ =1,IWBULL)
NOOOOD=0
END IF
NOOOOD=NOOOOD+1
IOOOO(1,NOOOOD)=I
IOOOO(2,NOOOOD)=J
IOOOO(3,NOOOOD)=K
IOOOO(4,NOOOOD)=L
HOOOO(NOOOOD)=HHH
ELSE
NVVVV=NVVVV+1
IF (NVVVVD.EQ.IWBULL) THEN
WRITE(IUVVVV) ((IVVVV(JDUM,IDUM),JDUM=1,4),HVVVV(IDUM),IDUM
$ =1,IWBULL)
NVVVVD=0
END IF
NVVVVD=NVVVVD+1
IVVVV(1,NVVVVD)=I
IVVVV(2,NVVVVD)=J
IVVVV(3,NVVVVD)=K
IVVVV(4,NVVVVD)=L
HVVVV(NVVVVD)=HHH
END IF
END IF
END IF
END IF
C (ov|ov)
IF ((LIO.NEQV.LJO).AND.(LKO.NEQV.LLO)) THEN
LHOLD=.TRUE.
NMP2=NMP2+1
NOVOV=NOVOV+1
END IF
C
IF (LHOLD) THEN
INDX=NMP2
IF (INDX.GT.NBULL) THEN
WRITE(6,*) ' REACHED BUFFER SIZE OF NBULL=',NBULL
STOP ' INCREASE NBULL '
END IF
C put the indices in the right order, without destroying indi, indj,
C indk or indl
CALL ORDIND(INDI,INDJ,INDK,INDL,I,J,K,L)
IH0(1,INDX)=I
IH0(2,INDX)=J
IH0(3,INDX)=K
IH0(4,INDX)=L
H0(INDX) =HHH
ELSE
IREJ=IREJ+1
END IF
END IF
RETURN
END
C
SUBROUTINE FABINX(II,JJ,KK,LL,INDX)
INCLUDE 'param.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
C
C I,J are assumed to be occupied, K,L virtual
C
I=II
J=JJ
K=KK
L=LL
CT IF (I.GT.NOCC) STOP ' FABINX: I.GT.NOCC '
CT IF (J.GT.NOCC) STOP ' FABINX: J.GT.NOCC '
CT IF (K.LE.NOCC) STOP ' FABINX: K.LE.NOCC '
CT IF (L.LE.NOCC) STOP ' FABINX: L.LE.NOCC '
IF (I.GT.J) THEN
IDUM=I
I=J
J=IDUM
END IF
IF (K.GT.L) THEN
IDUM=K
K=L
L=IDUM
END IF
NVO=NVIRT*(NVIRT+1)/2
INDX=(I*(I-1)/2+J)*NVO+K*(K-1)/2+L
RETURN
END
C
SUBROUTINE ORDIND(I,J,K,L,II,JJ,KK,LL)
C
C canonical order for i,j,k,l, i.e. i<=k and j<=l if i=k
C without destroying the variables i,j,k,l
C
II=I
JJ=J
KK=K
LL=L
IF (II.GT.KK) THEN
IDUM=II
II=KK
KK=IDUM
IDUM=JJ
JJ=LL
LL=IDUM
ELSE IF (II.EQ.KK) THEN
IF (JJ.GT.LL) THEN
IDUM=JJ
JJ=LL
LL=IDUM
END IF
END IF
RETURN
END
C
SUBROUTINE MP3BIS
INCLUDE 'param.h'
PARAMETER (NBULL=43000000)
COMMON /TWOI/ H0(NBULL),IH0(4,NBULL),NUMINT
C.. INCLUDE 'nbull.h'
COMMON /FLOW/ THRESH,ICUT,NFROZ,LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX
$ ,L3RDF,L4THF,LINFF,LCUT,LTHRE,LMART,L34EPS,LMP3
LOGICAL LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX,L3RDF,L4THF,LINFF,LCUT
$ ,LTHRE,LMART,L34EPS,LMP3
C.. INCLUDE 'flow.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCCS(NBASM),IOCC(NBASM)
C.. INCLUDE 'common_vect.h'
LOGICAL LIJ,LAB
C
C loop over (all!) indices
C
C we have to evaluate 12 diagrams
C
NOV=NOCC+1
EMP3=0.D0
E1 =0.D0
E2 =0.D0
E3 =0.D0
E4 =0.D0
E5 =0.D0
E6 =0.D0
C
DO I=1,NOCC
EI=ORBEN(I)
DO J=1,NOCC
EJ=ORBEN(J)
DO IA=NOV,NBAS
EA=ORBEN(IA)
DO IB=NOV,NBAS
EB=ORBEN(IB)
C
C the common term we calculate at the end
C
C first: loops over c and d
C
TCD=0.D0
DO IC=NOV,NBAS
EC=ORBEN(IC)
DO ID=NOV,NBAS
ED=ORBEN(ID)
EIJCD=EI+EJ-EC-ED
C (vv|vv)
HACBD=HFIND(IA,IC,IB,ID)
C (ov|ov)
HICJD=HFIND(I,IC,J,ID)
C
TCD=TCD+HACBD/EIJCD*HICJD
C
END DO
END DO
C
C next: loops over k and l
C
TKL=0.D0
DO K=1,NOCC
EK=ORBEN(K)
DO L=1,NOCC
EL=ORBEN(L)
EKLAB=EK+EL-EA-EB
C (oo|oo)
HIKJL=HFIND(I,K,J,L)
C (ov|ov)
HKALB=HFIND(K,IA,L,IB)
C
TKL=TKL+HIKJL/EKLAB*HKALB
C
END DO
END DO
C
C last: loops over k and c
C
TKC=0.D0
TKC1=0.D0
TKC2=0.D0
TKC3=0.D0
TKC4=0.D0
DO K=1,NOCC
EK=ORBEN(K)
DO IC=NOV,NBAS
EC=ORBEN(IC)
EIKAC=EI+EK-EA-EC
EIKBC=EI+EK-EB-EC
C
C (oo|vv)
HACJK=HFIND(J,K,IA,IC)
HBCJK=HFIND(J,K,IB,IC)
C (ov|ov)
HICKB=HFIND(I,IC,K,IB)
HIAKC=HFIND(I,IA,K,IC)
C
C diagrams 5-8
C
TERM1=-HACJK*HICKB/EIKBC
TERM2=-HBCJK*HIAKC/EIKAC
C
C diagrams 9-12
C
C only (ov|ov)
HIAKC=HFIND(I,IA,K,IC)
HICKA=HFIND(I,IC,K,IA)
HJBKC=HFIND(J,IB,K,IC)
C
TERM3=2.D0*HIAKC-HICKA
TERM3=TERM3*HJBKC/EIKAC
C
TKC=TKC+2.D0*(TERM1+TERM2+TERM3)
C
TKC1=TKC1+2.D0*TERM1
TKC2=TKC2+2.D0*TERM2
TKC3=TKC3+2.D0*TERM3
C
END DO
END DO
C
C the sum of all third-order diagrams for common i,j,a,b
C
THIRD=TKL+TCD+TKC
C
C we have always the common factor: [2 (ia|jb) - (ib|ja)]/(ei+ej-ea-eb)
C
EIJAB=EI+EJ-EA-EB
HIAJB=HFIND(I,IA,J,IB)
HIBJA=HFIND(I,IB,J,IA)
TIJAB=(2.D0*HIAJB-HIBJA)/EIJAB
EMP3=EMP3+TIJAB*THIRD
C
E1=E1+TKL*TIJAB
E2=E2+TCD*TIJAB
E3=E3+TKC1*TIJAB
E4=E4+TKC2*TIJAB
E5=E5+TKC3*TIJAB
C
C this was the common term
C
C Close loops over i,j,a,b
C
END DO
END DO
END DO
END DO
C
WRITE(6,*)
WRITE(6,*) ' DIAGRAMS 1+2 : ',E2
WRITE(6,*) ' DIAGRAMS 3+4 : ',E1
WRITE(6,*) ' DIAGRAMS 9+10+11+12: ',E5+E6
WRITE(6,*) ' DIAGRAMS 5+6+7+8 : ',E3+E4
WRITE(6,*)
WRITE(6,9913) '(MP3 LOC.) ',EMP3
C
9913 FORMAT(/,' CORRELATION ENERGY ',A10,':',F20.12,//)
C
RETURN
END
C
SUBROUTINE MP3
INCLUDE 'param.h'
PARAMETER (NBULL=43000000)
COMMON /TWOI/ H0(NBULL),IH0(4,NBULL),NUMINT
C.. INCLUDE 'nbull.h'
COMMON /FLOW/ THRESH,ICUT,NFROZ,LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX
$ ,L3RDF,L4THF,LINFF,LCUT,LTHRE,LMART,L34EPS,LMP3
LOGICAL LFROZ,LUNFOR,LMP2PR,LPRIFO,LPRCEX,L3RDF,L4THF,LINFF,LCUT
$ ,LTHRE,LMART,L34EPS,LMP3
C.. INCLUDE 'flow.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ HCOU(NBASM,NBASM),HEXC(NBASM,NBASM),
- F(NBASM,NBASM),ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /VECT/ CI(NBASM,NBASM),IOCCS(NBASM),IOCC(NBASM)
C.. INCLUDE 'common_vect.h'
COMMON /CNTERS/ NMP2,NMP3,NOOOO,NOVOV,NOOVV,NVVVV,NCOU,NEXC,IREJ
COMMON /BIBUF/HOOOO(IWBULL),HOOVV(IWBULL),HVVVV(IWBULL),IOOOO(4
$ ,IWBULL),IOOVV(4,IWBULL),IVVVV(4,IWBULL),NOOOOD,NOOVVD,NVVVVD
$ ,IUOOOO,IUOOVV,IUVVVV
C.. INCLUDE 'counters.h'
LOGICAL LIJ,LAB,LF
C
WRITE(6,*)
WRITE(6,*) ' THE MP3 CALCULATION '
WRITE(6,*)
$ ' (ov|ov) in memory, (oo|oo), (oo|vv), and (vv|vv) on disk'
WRITE(6,*)
C
C loop over (all!) indices
C
C we have to evaluate 12 diagrams
C
NOV=NOCC+1
C
C for the integrals (ov|ov) we loop as well over the list
C
D99=0.D0
DO INDX=1,NUMINT
I =IH0(1,INDX)
IA=IH0(2,INDX)
J =IH0(3,INDX)
IB=IH0(4,INDX)
HHH=H0(INDX)
C
C we have to divide by two the integrals (ia|ia)
C
IF (I.EQ.J.AND.IA.EQ.IB) HHH=0.5D0*HHH
C
EI=ORBEN(I)
EJ=ORBEN(J)
EA=ORBEN(IA)
EB=ORBEN(IB)
EIA=EI-EA
EJB=EJ-EB
C
DO K=1,NOCC
EK=ORBEN(K)
EIKA=EIA+EK
EJKB=EJB+EK
DO IC=NOV,NBAS
EC=ORBEN(IC)
EIKAC=EIKA-EC
EJKBC=EJKB-EC
DENOM=EIKAC*EJKBC
HIAKC=HFIND(I,IA,K,IC)
HICKA=HFIND(I,IC,K,IA)
HJBKC=HFIND(J,IB,K,IC)
HJCKB=HFIND(J,IC,K,IB)
T1=2.D0*HIAKC-HICKA
T2=2.D0*HJBKC-HJCKB
D99=D99+HHH*T1*T2/DENOM
END DO
END DO
END DO
D99=4.D0*D99
EMP3=D99
c$$$C
c$$$C for diagrams 9-12 we loop over ijabkc
c$$$C
c$$$ D99=0.D0
c$$$ DO I=1,NOCC
c$$$ EI=ORBEN(I)
c$$$ DO J=1,NOCC
c$$$ EJ=ORBEN(J)
c$$$ EIJ=EI+EJ
c$$$ LIJ=I.EQ.J
c$$$ DO IA=NOV,NBAS
c$$$ EA=ORBEN(IA)
c$$$ EIJA=EIJ-EA
c$$$ DO IB=NOV,NBAS
c$$$ LAB=IA.EQ.IB
c$$$ EB=ORBEN(IB)
c$$$ EIJAB=EIJA-EB
c$$$ EAB=-EA-EB
c$$$C
c$$$C the common term we calculate at the end
c$$$C
c$$$ DIA99=0.D0
c$$$ DO K=1,NOCC
c$$$ EK=ORBEN(K)
c$$$ EIK=EI+EK
c$$$ EIKA=EI+EK-EA
c$$$ EIKB=EI+EK-EB
c$$$ DO IC=NOV,NBAS
c$$$ EC=ORBEN(IC)
c$$$C
c$$$C diagrams 9-12
c$$$C
c$$$ EIKAC=EIKA-EC
c$$$ HIAKC=HFIND(I,IA,K,IC)
c$$$ HJBKC=HFIND(J,IB,K,IC)
c$$$ HICKA=HFIND(I,IC,K,IA)
c$$$ DIA99=DIA99+HJBKC/EIKAC*(2.D0*HIAKC-HICKA)
c$$$ END DO
c$$$ END DO
c$$$C
c$$$C we have always the common factor: [2 (ia|jb) - (ib|ja)]/(ei+ej-ea-eb)
c$$$C
c$$$ IF (LIJ) THEN
c$$$ IF (LAB) THEN
c$$$ TIJAB=HEXC(I,IA)/EIJAB
c$$$ D99=D99+TIJAB*DIA99
c$$$ ELSE
c$$$ TIJAB=HFIND(I,IA,I,IB)/EIJAB
c$$$ D99=D99+TIJAB*DIA99
c$$$ END IF
c$$$ ELSE
c$$$ IF (LAB) THEN
c$$$ TIJAB=HFIND(I,IA,J,IA)/EIJAB
c$$$ D99=D99+TIJAB*DIA99
c$$$ ELSE
c$$$ H1=HFIND(I,IA,J,IB)
c$$$ H2=HFIND(I,IB,J,IA)
c$$$ TIJAB=(2.D0*H1-H2)/EIJAB
c$$$ D99=D99+TIJAB*DIA99
c$$$ END IF
c$$$ END IF
c$$$C
c$$$C this was the common term
c$$$C
c$$$C Close loops over i,j,a,b
c$$$C
c$$$ END DO
c$$$ END DO
c$$$ END DO
c$$$ END DO
c$$$ D99 =2.D0*D99
c$$$ EMP3=D99
WRITE(6,*) ' diagrams 9 - 12 ',D99
CALL TIMING('9-12')
C
C diagrams 1+2 need integrals (vv|vv)
C
D12=0.D0
NRD=0
OPEN(UNIT=IUVVVV,FILE='VVVV',STATUS='OLD',FORM='UNFORMATTED',ERR
$ =802)
LF=.FALSE.
1002 CONTINUE
READ(IUVVVV) ((IVVVV(JDUM,IDUM),JDUM=1,4),HVVVV(IDUM),IDUM
$ =1,IWBULL)
IF (HVVVV(IWBULL).EQ.-10000.D0) THEN
MXIND=IVVVV(1,IWBULL)
LF=.TRUE.
ELSE
MXIND=IWBULL
END IF
NRD=NRD+MXIND
DO INDX=1,MXIND
IA=IVVVV(1,INDX)
IB=IVVVV(2,INDX)
IC=IVVVV(3,INDX)
ID=IVVVV(4,INDX)
HHH=HVVVV(INDX)
C
EA=ORBEN(IA)
EB=ORBEN(IB)
EC=ORBEN(IC)
ED=ORBEN(ID)
EAC=-EA-EC
EBC=-EB-EC
EAD=-EA-ED
EBD=-EB-ED
DO I=1,NOCC
EI=ORBEN(I)
EIAC=EAC+EI
EIAD=EAD+EI
EIBC=EBC+EI
EIBD=EBD+EI
DO J=1,NOCC
EJ=ORBEN(J)
EIJAC=EIAC+EJ
EIJAD=EIAD+EJ
EIJBC=EIBC+EJ
EIJBD=EIBD+EJ
DENOM1=EIJAC*EIJBD
DENOM2=EIJBC*EIJAD
HIAJC=HFIND(I,IA,J,IC)
HIAJD=HFIND(I,IA,J,ID)
HIBJC=HFIND(I,IB,J,IC)
HIBJD=HFIND(I,IB,J,ID)
HICJA=HFIND(I,IC,J,IA)
HIDJA=HFIND(I,ID,J,IA)
HICJB=HFIND(I,IC,J,IB)
HIDJB=HFIND(I,ID,J,IB)
T1A=HIAJC*(2.D0*HIBJD-HIDJB)
T1B=HICJA*(2.D0*HIDJB-HIBJD)
T1=(T1A+T1B)/DENOM1
T2A=HIBJC*(2.D0*HIAJD-HIDJA)
T2B=HICJB*(2.D0*HIDJA-HIAJD)
T2=(T2A+T2B)/DENOM2
D12=D12+HHH*(T1+T2)
END DO
END DO
C
END DO
IF (.NOT.LF) GO TO 1002
CLOSE(IUVVVV,STATUS='DELETE')
D12=2.D0*D12
EMP3=EMP3+D12
WRITE(6,*) ' read ',NRD,' integrals '
WRITE(6,*) ' diagrams 1 and 2 ',D12
CALL TIMING('1+2 ')
C
C diagrams 3+4 need integrals (oo|oo)
C
NRD=0
D34=0.D0
OPEN(UNIT=IUOOOO,FILE='OOOO',STATUS='OLD',FORM='UNFORMATTED',ERR
$ =801)
LF=.FALSE.
1001 CONTINUE
READ(IUOOOO) ((IOOOO(JDUM,IDUM),JDUM=1,4),HOOOO(IDUM),IDUM
$ =1,IWBULL)
IF (HOOOO(IWBULL).EQ.-10000.D0) THEN
MXIND=IOOOO(1,IWBULL)
LF=.TRUE.
ELSE
MXIND=IWBULL
END IF
NRD=NRD+MXIND
DO INDX=1,MXIND
C well, what do we have to do?
I=IOOOO(1,INDX)
J=IOOOO(2,INDX)
K=IOOOO(3,INDX)
L=IOOOO(4,INDX)
HHH=HOOOO(INDX)
C
C (ij|kl) can contribute to ijab, jiab, klab or lkab
C
EI=ORBEN(I)
EJ=ORBEN(J)
EK=ORBEN(K)
EL=ORBEN(L)
EIK=EI+EK
EJK=EJ+EK
EIL=EI+EL
EJL=EJ+EL
DO IA=NOV,NBAS
EA=ORBEN(IA)
EIKA=EIK-EA
EJKA=EJK-EA
EILA=EIL-EA
EJLA=EJL-EA
DO IB=NOV,NBAS
EB=ORBEN(IB)
EIKAB=EIKA-EB
EJKAB=EJKA-EB
EILAB=EILA-EB
EJLAB=EJLA-EB
DENOM1=EIKAB*EJLAB
DENOM2=EJKAB*EILAB
HIAKB=HFIND(I,IA,K,IB)
HIALB=HFIND(I,IA,L,IB)
HIBKA=HFIND(I,IB,K,IA)
HIBLA=HFIND(I,IB,L,IA)
HJAKB=HFIND(J,IA,K,IB)
HJALB=HFIND(J,IA,L,IB)
HJBKA=HFIND(J,IB,K,IA)
HJBLA=HFIND(J,IB,L,IA)
T1A=HIAKB*(2.D0*HJALB-HJBLA)
T1B=HIBKA*(2.D0*HJBLA-HJALB)
T1=(T1A+T1B)/DENOM1
T2A=HIALB*(2.D0*HJAKB-HJBKA)
T2B=HIBLA*(2.D0*HJBKA-HJAKB)
T2=(T2A+T2B)/DENOM2
D34=D34+HHH*(T1+T2)
END DO
END DO
C
END DO
IF (.NOT.LF) GO TO 1001
C
C we do not need the file any more
CLOSE(IUOOOO,STATUS='DELETE')
D34=2.D0*D34
EMP3=EMP3+D34
WRITE(6,*) ' read ',NRD,' integrals '
WRITE(6,*) ' diagrams 3 and 4 ',D34
CALL TIMING('3+4 ')
C
C
C diagrams 5-8 need integrals (oo|vv)
C
OPEN(UNIT=IUOOVV,FILE='OOVV',STATUS='OLD',FORM='UNFORMATTED',ERR
$ =803)
LF=.FALSE.
D56=0.D0
D78=0.D0
NRD=0
1003 CONTINUE
READ(IUOOVV) ((IOOVV(JDUM,IDUM),JDUM=1,4),HOOVV(IDUM),IDUM
$ =1,IWBULL)
IF (HOOVV(IWBULL).EQ.-10000.D0) THEN
MXIND=IOOVV(1,IWBULL)
LF=.TRUE.
ELSE
MXIND=IWBULL
END IF
NRD=NRD+MXIND
DO INDX=1,MXIND
I =IOOVV(1,INDX)
J =IOOVV(2,INDX)
IA=IOOVV(3,INDX)
IB=IOOVV(4,INDX)
HHH=HOOVV(INDX)
C
EI=ORBEN(I)
EJ=ORBEN(J)
EA=ORBEN(IA)
EB=ORBEN(IB)
EIA=EI-EA
EIB=EI-EB
EJA=EJ-EA
EJB=EJ-EB
DO K=1,NOCC
EK=ORBEN(K)
EIKA=EK+EIA
EIKB=EK+EIB
EJKA=EK+EJA
EJKB=EK+EJB
DO IC=NOV,NBAS
EC=ORBEN(IC)
EIKAC=EIKA-EC
EIKBC=EIKB-EC
EJKAC=EJKA-EC
EJKBC=EJKB-EC
C
DENOM1=EIKBC*EJKAC
DENOM2=EIKAC*EJKBC
C
HIAKC=HFIND(I,IA,K,IC)
HICKA=HFIND(I,IC,K,IA)
HIBKC=HFIND(I,IB,K,IC)
HICKB=HFIND(I,IC,K,IB)
HJAKC=HFIND(J,IA,K,IC)
HJCKA=HFIND(J,IC,K,IA)
HJBKC=HFIND(J,IB,K,IC)
HJCKB=HFIND(J,IC,K,IB)
C
T1A=(4.D0*HICKB-HIBKC)*HJCKA
T1B=HJAKC*HICKB
T1=(T1A-T1B)/DENOM1
T2A=(4.D0*HICKA-HIAKC)*HJCKB
T2B=HJBKC*HICKA
T2=(T2A-T2B)/DENOM2
D56=D56+HHH*(T1+T2)
C
T1A=(4.D0*HIBKC-HICKB)*HJAKC
T1B=HJCKA*HIBKC
T1=(T1A-T1B)/DENOM1
T2A=(4.D0*HIAKC-HICKA)*HJBKC
T2B=HJCKB*HIAKC
T2=(T2A-T2B)/DENOM2
D78=D78+HHH*(T1+T2)
END DO
END DO
C
END DO
IF (.NOT.LF) GO TO 1003
CLOSE(IUOOVV,STATUS='DELETE')
C
D56=-2.D0*D56
D78=-2.D0*D78
EMP3=EMP3+D56+D78
WRITE(6,*) ' read ',NRD,' integrals '
WRITE(6,*) ' diagrams 5 + 6 ',D56
WRITE(6,*) ' diagrams 7 + 8 ',D78
CALL TIMING('5..8')
C
WRITE(6,9913) '(MP3 LOC.) ',EMP3
C
9913 FORMAT(/,' CORRELATION ENERGY ',A10,':',F20.12,//)
C
C diagram 1+2: loops over cd, direct+exchange
C diagram 3+4: loops over kl, direct+exchange
C diagram 5+6: loops over kc, between the 2 bubbles
C diagram 5+6: loops over kc, inside 1 bubble
C diagram 9+10: 3 bubbles
C diagram 11+12: 1 small bubble
C
c$$$ WRITE(6,9914) D12,D34,D5678,D99
c$$$ 9914 FORMAT(/,' Decomposition in diagrams:',/,' Diagrams 1+2: ',F12.7,
c$$$ $ /,' Diagrams 3+4: ',F12.7,/,' Diagrams 5+6: ',F12.7,/
c$$$ $ ,' Diagrams 7+8: ',F12.7,/,' Diagrams 9-12: ',F12.7,/)
c$$$C
RETURN
C
801 CONTINUE
WRITE(6,*) ' what happened to file OOOO?'
STOP
802 CONTINUE
WRITE(6,*) ' what happened to file VVVV?'
STOP
803 CONTINUE
WRITE(6,*) ' what happened to file OOVV?'
STOP
C
RETURN
END
C
SUBROUTINE LEXSRT(IH0,H0,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
INTEGER IH0(4,N)
DIMENSION H0(N)
INTEGER IRA(4)
LOGICAL LCOMP4,LACT
C
C SORT BY IH0
C
IF (N.LT.2) GO TO 200
L=N/2+1
IR=N
10 CONTINUE
IF(L.GT.1)THEN
L=L-1
DO III=1,4
IRA(III)=IH0(III,L)
END DO
RH0=H0(L)
ELSE
DO III=1,4
IRA(III)=IH0(III,IR)
IH0(III,IR)=IH0(III,1)
END DO
RH0=H0(IR)
H0(IR)=H0(1)
IR=IR-1
IF (IR.EQ.1) THEN
DO III=1,4
IH0(III,1)=IRA(III)
END DO
H0(1)=RH0
GO TO 200
END IF
END IF
I=L
J=L+L
20 IF(J.LE.IR)THEN
IF (J.LT.IR)THEN
LACT=LCOMP4(IH0(1,J),IH0(1,J+1))
IF (LACT) THEN
J=J+1
END IF
END IF
LACT=LCOMP4(IRA,IH0(1,J))
IF (LACT) THEN
DO III=1,4
IH0(III,I)=IH0(III,J)
END DO
H0(I)=H0(J)
I=J
J=J+J
ELSE
J=IR+1
END IF
GOTO 20
END IF
DO III=1,4
IH0(III,I)=IRA(III)
END DO
H0(I)=RH0
GOTO 10
C
200 CONTINUE
RETURN
END
C
FUNCTION LCOMP4(IARR,JARR)
LOGICAL LCOMP4
INTEGER IARR(4),JARR(4)
C
IF (IARR(1).LT.JARR(1)) THEN
LCOMP4=.TRUE.
ELSE IF (IARR(1).GT.JARR(1)) THEN
LCOMP4=.FALSE.
ELSE
IF (IARR(2).LT.JARR(2)) THEN
LCOMP4=.TRUE.
ELSE IF (IARR(2).GT.JARR(2)) THEN
LCOMP4=.FALSE.
ELSE
IF (IARR(3).LT.JARR(3)) THEN
LCOMP4=.TRUE.
ELSE IF (IARR(3).GT.JARR(3)) THEN
LCOMP4=.FALSE.
ELSE
IF (IARR(4).LT.JARR(4)) THEN
LCOMP4=.TRUE.
ELSE IF (IARR(4).GT.JARR(4)) THEN
LCOMP4=.FALSE.
ELSE
WRITE(6,*) 'IARR: ',(IARR(JJJ),JJJ=1,4),' JARR: ',(JARR(KKK)
$ ,KKK=1,4)
STOP ' EQUAL INDICES ENCOUNTERED! '
END IF
END IF
END IF
END IF
C
RETURN
END
C