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