c-*- Fortran -*-
C
C DELCI formatted intermediates
C DELCO print info on output on temporary files
C DELCY simulation of reading, no actual I/O
C DELCD Debug
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Four-index transformation
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 'common_syst.h'
C..FILE 'common_intu.h'
C..FILE 'common_vectco.h'
C..FILE 'common_two.h'
C..FILE 'common_intcount.h'
C
C..FILE 'common_twosplit.h'
C
C..FILE 'flow.h'
C..FILE 'common_wrk.h'
C
PROGRAM MAIN
C
C 4-index transformation of two-electron integrals
C out-of-core, integral-driven
C
INCLUDE 'param.h'
COMMON /FLOW/ THRESH1,THRESH2,THRESH3,LINFO,LOUTFO,LDIAGF
- ,LMONO,LBIEL,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE
$ ,LFROZ
C LSTH: store half-transformed integrals somewhere
LOGICAL LINFO,LOUTFO,LDIAGF,LMONO,LBIEL
$ ,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE,LFROZ
C.. INCLUDE 'flow.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
COMMON /ENERG/ EONE,ETWO,ECOUL,EEXCH
C
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ S(NBASM,NBASM),F(NBASM,NBASM),HONE(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /VECTCO/ CI(NBASM,NBASM),IOCCS(NBASM),IOCC(NBASM)
C.. INCLUDE 'common_vectco.h'
COMMON /INTCNT/ XNAO,XNAOTMP,XNHFTR,XNMO,XNINTIN,XNINTFF,XIRD
$ ,XISTOR,XUNIQ
C.. INCLUDE 'common_intcount.h'
DIMENSION IDEC(NBASM*NBASM)
DIMENSION HHH(NBASM,NBASM)
DIMENSION FREF(NBASM*(NBASM+1)/2),EW(NBASM),EV(NBASM,NBASM)
$ ,P1(NBASM),P2(NBASM)
DIMENSION CTMP(NBASM)
C
INCLUDE 'compiler_stamp'
X=CPTIME(3)
CALL DATING('VIND ',1)
WRITE(6,*)
WRITE(6,*) ' V I N D - simple molecule'
WRITE(6,*) ' 4-index transformation, integral-driven '
WRITE(6,*)
WRITE(6,*) ' P.Reinhardt Dresden 8/1999, Paris 12/2000'
WRITE(6,*) ' with integral prefactors 4/2001 '
WRITE(6,*) ' direct algorithm 10/2003 '
WRITE(6,*)
C
CALL RDINP
C
LSTH=.FALSE.
WRITE(6,*) ' Scanning bielectronic integrals ','on input against'
$ ,THRESH1
WRITE(6,*) ' Scanning half-transformed integrals ',' against '
$ ,THRESH2
IF (LALL) THEN
WRITE(6,*)
$ ' GENERATING ALL TWO-ELECTRON INTEGRALS in MOs'
WRITE(6,*) ' including all zeroes '
THRESH3=0.D0
ELSE
WRITE(6,*) ' Scanning bielectronic integrals '
$ ,'on output against',THRESH3
END IF
IF (LBIEL) THEN
WRITE(6,*) ' BIELECTRONIC INTEGRALS WILL BE GENERATED'
ELSE
WRITE(6,*) ' NO BIELECTRONIC INTEGRALS WILL BE GENERATED'
END IF
IF (LMONO) THEN
WRITE(6,*)
$ ' MONOELECTRONIC INTEGRALS WILL BE GENERATED'
ELSE
WRITE(6,*)
$ ' NO MONOELECTRONIC INTEGRALS WILL BE GENERATED'
END IF
IF (LINFO) WRITE(6,*) ' EXPECTING FORMATTED INTEGRALS TO READ '
IF (LOUTFO) THEN
WRITE(6,*) ' at the end of the transformation : '
WRITE(6,*) ' READING FILE53 and writing formatted file MOTWO_FORM
$ATTED'
C folded 1 (fixf $Revision: 1.3 $)
ELSE
WRITE(6,*) ' WRITING UNFORMATTED INTEGRALS, BUFFER LENGTH '
$ ,IWBULL
END IF
IF (LDIAGF) WRITE(6,*)
$ ' DIAGONALIZING THE FOCK OPERATOR IN THE REFERENCE CELL'
WRITE(6,*)
IF (LNORM) THEN
WRITE(6,*) ' Normalizing the molecular orbitals '
ELSE
WRITE(6,*) ' No normalization of the molecular orbitals '
END IF
IF (LOVOV) THEN
WRITE(6,*) ' generating only (ov|ov) integrals '
END IF
IF (LDELE) THEN
WRITE(6,*) ' deleting virtual orbitals from the transformation '
WRITE(6,*) ' not yet operational, work harder '
STOP
END IF
IF (LFROZ) THEN
WRITE(6,*) ' freezing occupied orbitals '
WRITE(6,*) ' not yet operational, work harder '
STOP
END IF
C
C two-electron integrals, output
IUNIT4=53
C read ftn26
IUNITR=26
C read information on system
C
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
C
WRITE(6,*)
WRITE(6,*) ' READING INFORMATION ON SYSTEM FROM UNIT',IUNITR
WRITE(6,9946) NATOM,NATMX,NBAS,NBASM
9946 FORMAT(' ACTUAL MAXIMUM ',/,
- ' NATOMS ',2I8,/,
- ' NBAS ',2I8)
IF (NBAS.GT.NBASM) THEN
WRITE(6,*) ' NBAS (',NBAS,' ) > NBASM (',NBASM,' )'
STOP 'INCREASE NBASM'
END IF
C
C read vector
C
IUNITO=31
OPEN(UNIT=IUNITO,FILE='VECTOR',STATUS='OLD',
- FORM='FORMATTED',ERR=993)
WRITE(6,*)
WRITE(6,*) ' READING VECTOR FROM UNIT',IUNITO
$ ,' FILE '
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
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,*)
C
C get overlap matrix
C
WRITE(6,*) ' reading OVERLAP '
IUNITR=14
OPEN(UNIT=IUNITR,FILE='OVERLAP',STATUS='OLD',ERR=992,
- FORM='FORMATTED')
721 CONTINUE
READ(IUNITR,*,IOSTAT=KK) I,J,HDUM
IF (KK.NE.0) GO TO 722
S(I,J)=HDUM
S(J,I)=HDUM
GO TO 721
722 CONTINUE
CLOSE(IUNITR)
IF (LNORM) THEN
CALL TRANSF(S,CI,NBAS)
DO I=1,NBAS
SFACT=1.D0/SQRT(S(I,I))
DO IALPH=1,NBAS
CI(IALPH,I)=CI(IALPH,I)*SFACT
END DO
END DO
END IF
C
C read Fock matrix
C
IUNITF=23
OPEN(UNIT=IUNITF,STATUS='OLD',FORM='FORMATTED',
- FILE='FOCK',ERR=7811)
READ(IUNITF,*) IDUM1
WRITE(6,*) ' READING FOCK MATRIX FROM UNIT',IUNITF
READ(IUNITF,*) ((F(I,J),I=1,NBAS),J=1,NBAS)
READ(IUNITF,*) EDUM,EN,E1,E2
CLOSE(IUNITF)
OPEN(UNIT=IUNITF,FILE='ENUC',FORM='FORMATTED',STATUS='UNKNOWN')
WRITE(IUNITF,*) EN
CLOSE(IUNITF)
C
C F in MOs
C
CALL TRANSF(F,CI,NBAS)
write(6,*) ' transforming '
C
WRITE(6,*)
WRITE(6,*) ' TRANSFORMED FOCK MATRIX '
WRITE(6,*) ' ======================= '
WRITE(6,*)
OPEN (UNIT=17,FILE='FOCK.MO',STATUS='UNKNOWN',FORM='FORMATTED')
DO I=1,NBAS
DO J=I,NBAS
IF (ABS(F(I,J)).GT.1.D-12) WRITE(17,'(2I6,E20.12)') I,J,F(I,J)
END DO
END DO
CLOSE(16)
C CALL PFUNC(F,NBAS,NBAS)
WRITE(6,*)
DO I=1,NBAS
ORBEN(I)=F(I,I)
END DO
WRITE(6,*) ' ORBEN: '
WRITE(6,'(4(I4,F14.5))') (I,ORBEN(I),I=1,NBAS)
WRITE(6,*)
FMX=0.D0
DO I=1,NOCC
DO J=NOCC+1,NBAS
FF=ABS(F(I,J))
IF (FF.GT.FMX) FMX=FF
END DO
END DO
WRITE(6,8102) FMX
8102 FORMAT(' LARGEST ELEMENT BETWEEN OCC/VIRT: ',E12.6,/)
WRITE(6,*)
CALL TIMING('FOCK')
GO TO 7812
7811 CONTINUE
WRITE(6,*)
WRITE(6,*) ' NO FOCK MATRIX FOUND . . . . '
WRITE(6,*)
EN=0.D0
E1=0.D0
E2=0.D0
C
C if we do not find a Fock matrix, we renormalize the vectors
C they might have come from a Slater calculation and want to be
C symmetry-adapted
C
7812 CONTINUE
C
IF (LMONO) THEN
C
C ONE-ELECTRON INTEGRALS
C read one-electron integrals
C
C Integrals are on files KINETIC, HAMILTON
C we accumulate in HONE
C
C KINETIC
C
CALL AOTOMO('OVERLAP',CI,S,NBAS)
CALL AOTOMO('KINETIC',CI,HONE,NBAS)
C
C on HAMILTO is E-N + E_kin
C
CALL AOTOMO('HAMILTO',CI,HONE,NBAS)
C
WRITE(6,*)
WRITE(6,*) ' ONE-ELECTRON INTEGRALS IN MOs'
EONE=0.D0
DO I=1,NOCC
WRITE(6,*) ' ',I,HONE(I,I)
EONE=EONE+HONE(I,I)
END DO
DO I=1+NOCC,NBAS
WRITE(6,*) ' ',I,HONE(I,I)
END DO
EONE=EONE*2.D0
WRITE(6,*) ' RECALCULATED ONE-ELECTRON ENERGY:',EONE
WRITE(6,*)
CALL TIMING('ONE ')
END IF
C
IF (LBIEL) THEN
C
IUNITK=55
OPEN(UNIT=IUNITK,STATUS='UNKNOWN',FORM='FORMATTED',FILE
$ ='CORR.INPDAT')
WRITE(IUNITK,'(4I8)') NBAS,NOCC
CLOSE(IUNITK)
C
C count integrals
C
NBAO=NBAS
NBMO=NBAS
C
WRITE(6,*)
WRITE(6,*) ' SIZE PARAMETERS '
WRITE(6,*) ' NBAO = ',NBAO
WRITE(6,*) ' NBMO = ',NBMO
C
NTRIA=NBAO*(NBAO+1)/2
NTRIB=NBMO*(NBMO+1)/2
C in AOs
XNAO=DBLE(NTRIA)*(DBLE(NTRIA)+1.D0)*0.5D0
XNAOTMP=DBLE(NTRIA)*DBLE(NTRIA)
C half-transformed
XNHFTR=DBLE(NTRIA)*DBLE(NTRIB)
C in MOs
XNMO=DBLE(NTRIB)*(DBLE(NTRIB)+1.D0)*0.5D0
C
WRITE(6,*) ' WE MIGHT HAVE ',XNAO,' UNIQUE AO INTEGRALS'
WRITE(6,*) ' FROM THESE WE ARRIVE AT ',XNAOTMP
$ ,' INTERMEDIATE AO INTEGRALS'
WRITE(6,*) ' WE MIGHT HAVE ',XNHFTR,' HALFTRANSFORMED INTEGRALS'
WRITE(6,*) ' WE MIGHT HAVE ',XNMO,' MO INTEGRALS'
C
IF (LDAF) THEN
CALL TRANSDA
ELSE
CALL TRANSBI
END IF
CALL TIMING('TWOI')
END IF
WRITE(6,9908)
9908 FORMAT(/,' Integrals statistics',//
$ ,' type max. real ratio (in %)')
9909 FORMAT(3X,A7,2F13.0,F9.2)
WRITE(6,9909) ' AO ',XNAO,XIRD,XIRD/XNAO
WRITE(6,9909) ' AOTMP ',XNAOTMP,XIRD,XIRD/XNAO
WRITE(6,9909) 'HALFTRF',XNHFTR,XNINTIN,XNINTIN/XNHFTR
WRITE(6,*)
CALL TIMING('STAT')
WRITE(6,*)
C
ETOT=EONE+ETWO+EN
C
WRITE(6,9917) ETOT,EONE,ETWO,EN,E1+E2+EN,E1,E2
9917 FORMAT(/,' RECALCULATED TOTAL ENERGY ',F20.12,/
$ ,' RECALCULATED ONE-ELECTRON ENERGY ',F20.12,/
$ ,' RECALCULATED TWO-ELECTRON ENERGY ',F20.12,/
$ ,' NUCLEAR REPULSION ',F20.12,//
$ ,' READ TOTAL ENERGY ',F20.12,/
$ ,' READ ONE-ELECTRON ENERGY ',F20.12,/
$ ,' READ TWO-ELECTRON ENERGY ',F20.12)
WRITE(6,*) ' DIFFERENCES '
WRITE(6,*) ' ONE EL. : ',EONE-E1
WRITE(6,*) ' TWO EL. : ',ETWO-E2
WRITE(6,*) ' TOTAL : ',ETOT-E1-E2-EN
WRITE(6,*)
C
CALL TIMING('ENER')
WRITE(6,*)
IF (LOUTFO) THEN
WRITE(6,*) ' we transform the bielectronic integrals '
WRITE(6,*) ' from unformatted to a formatted file'
CALL READ53(53)
END IF
C
WRITE(6,*)
X=CPTIME(4)
CALL DATING('VIND ',2)
STOP
901 CONTINUE
WRITE(6,*) ' NO FILE FOUND ... '
STOP ' PLEASE SUPPLY FILE! '
991 CONTINUE
WRITE(6,*) ' NO FILE FOUND ... '
STOP ' PLEASE SUPPLY FILE! '
992 CONTINUE
WRITE(6,*) ' NO FILE FOUND ... '
STOP ' PLEASE SUPPLY FILE! '
993 CONTINUE
WRITE(6,*) ' NO FILE FOUND ... '
STOP ' PLEASE SUPPLY FILE! '
994 CONTINUE
WRITE(6,*) ' NO FILE FOUND ... '
STOP ' PLEASE SUPPLY FILE! '
995 CONTINUE
WRITE(6,*) ' NO FILE FOUND ... '
STOP ' PLEASE SUPPLY FILE! '
END
C
SUBROUTINE TRANSF(A,CI,NBAS)
INCLUDE 'param.h'
DIMENSION AJKLM(NBASM,NBASM)
DIMENSION A(NBASM,NBASM),CI(NBASM,NBASM)
C
DO J=1,NBAS
DO IALPH=1,NBAS
AJKLM(IALPH,J)=0.D0
END DO
END DO
C
DO J=1,NBAS
DO IBETA=1,NBAS
BB=0.D0
DO IALPH=1,NBAS
BB=BB+CI(IALPH,J)*A(IBETA,IALPH)
END DO
AJKLM(IBETA,J)=BB
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 IBETA=1,NBAS
SSS=SSS+AJKLM(IBETA,J)*CI(IBETA,I)
END DO
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,NBAS)
INCLUDE 'param.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 IPOS(I,J,K,L,NBAS)
ICJ=NBAS*NBAS
ICVJ=ICJ*NBAS
ICI=ICVJ
ICVL=NBAS
ICK=ICVL
ICVK=ICK*NBAS
ICONST=L-ICVL-ICK-ICVK-ICJ-ICVJ-ICI
IPOS=ICONST+ICVL+ICK*K+ICVK+ICJ*J+ICVJ+ICI*I
RETURN
END
C
SUBROUTINE RDINP
INCLUDE 'param.h'
COMMON /FLOW/ THRESH1,THRESH2,THRESH3,LINFO,LOUTFO,LDIAGF
- ,LMONO,LBIEL,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE
$ ,LFROZ
C LSTH: store half-transformed integrals somewhere
LOGICAL LINFO,LOUTFO,LDIAGF,LMONO,LBIEL
$ ,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE,LFROZ
C.. INCLUDE 'flow.h'
CHARACTER*4 KEYW(2),STR3
CHARACTER*6 KEYOPT(15),STR6
CHARACTER*80 LINE
DATA KEYW /'*VIN','*END'/
DATA KEYOPT /'NONORM','INFORM','OUTFOR','DIAGFO','THRESH',
- 'NOMONO','NOBIEL','DIRECT','ALL ','THRE1 ',
- 'THRE2 ','THRE3 ','OVOV ','DELETE','FREEZE'/
C
C DEFAULTS
THRESH1=1.D-10
THRESH2=1.D-10
THRESH3=1.D-10
LINFO=.FALSE.
LOUTFO=.FALSE.
LDIAGF=.FALSE.
LMONO=.TRUE.
LBIEL=.TRUE.
LDAF=.FALSE.
LALL=.FALSE.
LNORM=.TRUE.
LOVOV=.FALSE.
LDELE=.FALSE.
LFROZ=.FALSE.
IOINP=83
OPEN(UNIT=IOINP,FILE='INPUT.VIN',ERR=2217,FORM='FORMATTED',
- STATUS='OLD')
C first, look for keyword '*VIN'
1 CONTINUE
READ(IOINP,'(A4)',END=22,ERR=921) STR3
IF (STR3.EQ.KEYW(1)) THEN
11 CONTINUE
READ(IOINP,'(A6)',END=23,ERR=921) STR6
c$$$ WRITE(6,*) ' READ STRING: |',STR6,'|'
C keyword *END
IF (STR6(1:4).EQ.KEYW(2)) GO TO 2
IF (STR6.EQ.KEYOPT(1)) THEN
C no norm of the vector
LNORM=.FALSE.
ELSE IF (STR6.EQ.KEYOPT(2)) THEN
C we try to read bielectronic integrals formatted
LINFO=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(3)) THEN
C we will write transformed integrals formatted
LOUTFO=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(4)) THEN
C Fock operator will be diagonalized in the unit cell
LDIAGF=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(5)) THEN
C a global threshold for integral scanning will be read
READ(IOINP,*) ITHRE
THRESH1=10.D0**DBLE(-ITHRE)
THRESH2=10.D0**DBLE(-ITHRE)
THRESH3=10.D0**DBLE(-ITHRE)
ELSE IF (STR6.EQ.KEYOPT(6)) THEN
C no monoelectronic integrals
LMONO=.FALSE.
ELSE IF (STR6.EQ.KEYOPT(7)) THEN
C no bielectronic integrals
LBIEL=.FALSE.
ELSE IF (STR6.EQ.KEYOPT(8)) THEN
C transformation via direct-access file
LDAF=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(9)) THEN
C all integrals to MO file
LALL=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(10)) THEN
C integral threshold for input
READ(IOINP,*) ITHRE
THRESH1=10.D0**DBLE(-ITHRE)
ELSE IF (STR6.EQ.KEYOPT(11)) THEN
C integral threshold for half-transformed integrals
READ(IOINP,*) ITHRE
THRESH2=10.D0**DBLE(-ITHRE)
ELSE IF (STR6.EQ.KEYOPT(12)) THEN
C integral threshold for output
READ(IOINP,*) ITHRE
THRESH3=10.D0**DBLE(-ITHRE)
ELSE IF (STR6.EQ.KEYOPT(13)) THEN
C only (ov|ov) integrals
LOVOV=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(14)) THEN
C delete virtual orbitals
LDELE=.TRUE.
ELSE IF (STR6.EQ.KEYOPT(15)) THEN
C freeze occupied orbitals
LFROZ=.TRUE.
ELSE
WRITE(6,*) ' POSSIBLE OPTIONS IN THE BLOCK *VIN ... *END ARE:'
WRITE(6,*)
WRITE(6,*) ' ',KEYOPT(2)
WRITE(6,*) ' ',KEYOPT(3)
WRITE(6,*) ' ',KEYOPT(4)
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)
WRITE(6,*) ' ',KEYOPT(14)
WRITE(6,*) ' ',KEYOPT(15)
STOP ' CHOOSE CORRECT OPTION '
END IF
GO TO 11
END IF
GO TO 1
2 CONTINUE
CLOSE(IOINP)
RETURN
921 CONTINUE
CLOSE(IOINP)
STOP ' ERROR IN INPUT'
23 CONTINUE
CLOSE(IOINP)
WRITE(6,*) ' YOU SHOULD TERMINATE THE BLOCK VIN BY <*END> '
RETURN
22 CONTINUE
WRITE(6,*) ' NO KEYWORD <*VIN>, USING THE DEFAULT VALUES'
RETURN
2217 CONTINUE
WRITE(6,*) ' NO FILE , USING THE DEFAULT VALUES'
RETURN
END
C
FUNCTION IMAP(I,IS)
IF (IS.EQ.1) THEN
IF (I.EQ.1) THEN
IMAP=1
ELSE
IMAP=3
END IF
ELSE
IF (I.EQ.2) THEN
IMAP=2
ELSE
IMAP=4
END IF
END IF
RETURN
END
C
SUBROUTINE TRANSBI
INCLUDE 'param.h'
C the buffer for transformed bielectronic integrals (oo|oo)
C we split files into maximally 20 smaller files
PARAMETER (NSPLIT=20,NBUL2=36000000,NLEVEL=5)
C NBUFL is the buffer length NBUL2/NSPLIT
COMMON /TWOSPL/ H0,IH0,NBUFW,NINDXB,INCORE,ISTBUF,IOFILE,NINDXF
$ ,NINDXG,NBUFL
C the buffer for integrals
DIMENSION H0(NBUL2)
C the buffer for indices
DIMENSION IH0(4,NBUL2)
C the two pointers for each file - how many buffers are on disk, how
C many elements are still in-core, we may use maximally 10 different
C levels of these pointers. ISTBUF contains the start addresses of the
C buffers
C NINDXB
DIMENSION NBUFW(NSPLIT,NLEVEL),INCORE(NSPLIT,NLEVEL)
- ,ISTBUF(NSPLIT),NINDXB(0:NLEVEL+1),IOFILE(NSPLIT)
- ,IPAIRS(NSPLIT,NLEVEL),IPAIRF(NSPLIT,NLEVEL),NINDXF(NLEVEL)
$ ,NINDXG(NLEVEL)
C.. INCLUDE 'common_twosplit.h'
C
COMMON /WBUF/ HBUF(IWBULL),IBUF(IWBULL),JBUF(IWBULL)
$ ,KBUF(IWBULL),LBUF(IWBULL),INTRD,IWFILE,IBUFL
COMMON /CNTRS/ IWCNT,IMOCNT,IMOREJ,LREAD
LOGICAL LREAD
C.. INCLUDE 'common_two.h'
COMMON /VECTCO/ CI(NBASM,NBASM),IOCCS(NBASM),IOCC(NBASM)
C.. INCLUDE 'common_vectco.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
COMMON /ENERG/ EONE,ETWO,ECOUL,EEXCH
C
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ S(NBASM,NBASM),F(NBASM,NBASM),HONE(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /FLOW/ THRESH1,THRESH2,THRESH3,LINFO,LOUTFO,LDIAGF
- ,LMONO,LBIEL,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE
$ ,LFROZ
C LSTH: store half-transformed integrals somewhere
LOGICAL LINFO,LOUTFO,LDIAGF,LMONO,LBIEL
$ ,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE,LFROZ
C.. INCLUDE 'flow.h'
LOGICAL LIK,LCD
INTEGER*2 IDT(4)
C
CHARACTER*1 FN(0:9)
DATA FN /'0','1','2','3','4','5','6','7','8','9'/
C
ETWO=0.D0
C
C the two-electron integrals
C
C for the moment square matrices
NBAO=NBAS
NBMO=NBAS
C how many MO integrals can we have?
C
WRITE(6,*)
WRITE(6,*) ' SIZE PARAMETERS '
WRITE(6,*) ' NBAO = ',NBAO
WRITE(6,*) ' NBMO = ',NBMO
WRITE(6,*) ' NBAS = ',NBAS
C
NTRIA=NBAO*(NBAO+1)/2
IF (LOVOV) THEN
NTRIB=NOCC*(NOCC+1)/2+NOCC*(NBAS-NOCC)
IAFIN=NOCC
IBFIN=NBMO
ELSE
NTRIB=NBMO*(NBMO+1)/2
IAFIN=NBMO
IBFIN=NBMO
END IF
C in AOs
XNAO=DBLE(NTRIA)*(DBLE(NTRIA)+1.D0)*0.5D0
XNAOTMP=DBLE(NTRIA)*DBLE(NTRIA)
C half-transformed
XNHFTR=DBLE(NTRIA)*DBLE(NTRIB)
C in MOs
XNMO=DBLE(NTRIB)*(DBLE(NTRIB)+1.D0)*0.5D0
C
c$$$ WRITE(6,*) ' WE MIGHT HAVE ',XNAO,' UNIQUE AO INTEGRALS'
C
C the main buffer has a length of NBUL2, we split it into NSPLIT
C files; each segment must hold at least NBAS*NBAS entries
C
NBUFL=INTDIV(NBUL2,NSPLIT)
IF (NBUFL.LT.NBAO*NBAO) THEN
WRITE(6,*)
$ ' TRANSBI: the buffer NBUL2 is too short: for this problem '
$
WRITE(6,*) ' we need ',NBAO*NBAO,' entries in each segment,'
WRITE(6,*) ' but we have only ',NBUFL,' entries available'
WRITE(6,*) ' either you augment NBUL2 or you diminish NSPLIT '
STOP ' buffer too small '
ELSE
C we divide the buffer in NSPLIT parts
ISTBUF(1)=NBUFL+1
C
C we start at logical unit 54, as 53 is used for the fully transformed integrals
C
IOFILE(1)=54
IONE=1
C and associate I/O units to the buffers
DO IFIL=2,NSPLIT-1
ISTBUF(IFIL)=ISTBUF(IFIL-1)+NBUFL
WRITE(6,9811) IFIL-1,ISTBUF(IFIL-1),ISTBUF(IFIL)-1
9811 FORMAT(' buffer No ',I2,' starts at ',I8,' and ends at address '
$ ,I8)
IOFILE(IFIL)=IOFILE(IFIL-1)+1
END DO
WRITE(6,9811) NSPLIT-1,ISTBUF(NSPLIT-1),NBUL2
C
END IF
C we have to store nbas*(nbas+1)/2 indices (gamma,delta) on NSPLIT-1
c files making it
C
NFILCD=INTDIV(NTRIA,NSPLIT-1)
NFTMP=INTDIV(NTRIA,NFILCD)
WRITE(6,*)
WRITE(6,*) ' First decomposition: we store ',NFILCD
$ ,' couples (gamma,delta) on one file'
WRITE(6,*) ' we need to open ',NFTMP
$ ,' intermediate files '
WRITE(6,*)
IFII =1
IFIJ =1
IFILE=1
INDX=0
INDX2=0
INDX3=0
DO I=1,NBAO
DO J=I,NBAO
IFILE2=(ITRANG(I,J,NBAO)-1)/NFILCD+1
IF (IFILE.NE.IFILE2) THEN
WRITE(6,9821) IFILE,IFII,IFIJ,IFII2,IFIJ2,INDX
9821 FORMAT(' File No ',I3,' has indices starting at ',2I4
$ ,' and ending at ',2I4,'; total of ',I8,' index couples')
C these are the new starting indices
IFII=I
IFIJ=J
IFILE=IFILE2
INDX3=INDX3+INDX
INDX=0
END IF
INDX=INDX+1
INDX2=INDX2+1
IFII2=I
IFIJ2=J
END DO
END DO
C the last file
INDX3=INDX3+INDX
WRITE(6,9821) IFILE,IFII,IFIJ,IFII2,IFIJ2,INDX
WRITE(6,*) ' NTRIA = ',NTRIA,' INDX2 = ',INDX2,' INDX3 = ',INDX3
WRITE(6,*)
C how deep do we have to decompose the bielectronic integrals file?
C
C this we should test before starting to work
C
NCD=NTRIA
N1ST=INTDIV(NCD ,NSPLIT-1)
N2ND=INTDIV(N1ST,NSPLIT-1)
N3RD=INTDIV(N2ND,NSPLIT-1)
N4TH=INTDIV(N3RD,NSPLIT-1)
N5TH=INTDIV(N4TH,NSPLIT-1)
c$$$ WRITE(6,*) ' N1ST = ',N1ST
c$$$ WRITE(6,*) ' N2ND = ',N2ND
c$$$ WRITE(6,*) ' N3RD = ',N3RD
c$$$ WRITE(6,*) ' N4TH = ',N4TH
c$$$ WRITE(6,*) ' N5TH = ',N5TH
NINDXF(1)=N1ST
NINDXF(2)=N2ND
NINDXF(3)=N3RD
NINDXF(4)=N4TH
NINDXF(5)=N5TH
C
IF (N1ST.EQ.1) THEN
WRITE(6,*) ' WE NEED ONE SINGLE DECOMPOSITION OF'
$ ,' THE BIELECTRONIC INTEGRALS FILE '
WRITE(6,*) ' each primary file contains up to ',N1ST
$ ,' index couples'
NDECO=0
ELSE
IF (N2ND.EQ.1) THEN
WRITE(6,*) ' WE NEED TWO SUCCESSIVE DECOMPOSITIONS OF'
$ ,' THE BIELECTRONIC INTEGRALS FILE '
NDECO=1
WRITE(6,*) ' each primary file contains up to ',N1ST
- ,' index couples'
WRITE(6,*) ' each secondary file contains up to ',N2ND
$ ,' index couples'
ELSE
IF (N3RD.EQ.1) THEN
WRITE(6,*) ' WE NEED THREE SUCCESSIVE DECOMPOSITIONS OF'
$ ,' THE BIELECTRONIC INTEGRALS FILE '
NDECO=2
WRITE(6,*) ' each primary file contains up to ',N1ST
$ ,' index couples'
WRITE(6,*) ' each secondary file contains up to ',N2ND
$ ,' index couples'
WRITE(6,*) ' each third-level file contains up to ',N3RD
$ ,' index couples'
ELSE
IF (N4TH.EQ.1) THEN
WRITE(6,*) ' WE NEED FOUR SUCCESSIVE DECOMPOSITIONS OF'
$ ,' THE BIELECTRONIC INTEGRALS FILE '
NDECO=3
WRITE(6,*) ' each primary file contains up to ',N1ST
$ ,' index couples'
WRITE(6,*) ' each secondary file contains up to ',N2ND
$ ,' index couples'
WRITE(6,*) ' each third-level file contains up to ',N3RD
$ ,' index couples'
WRITE(6,*) ' each fourth-level file contains up to ',N4TH
$ ,' index couples'
ELSE
IF (N5TH.EQ.1) THEN
WRITE(6,*) ' WE NEED FIVE SUCCESSIVE DECOMPOSITIONS OF'
$ ,' THE BIELECTRONIC INTEGRALS FILE '
NDECO=4
WRITE(6,*) ' each primary file contains up to ',N1ST
$ ,' index couples'
WRITE(6,*) ' each secondary file contains up to ',N2ND
$ ,' index couples'
WRITE(6,*) ' each third-level file contains up to ',N3RD
$ ,' index couples'
WRITE(6,*) ' each fourth-level file contains up to ',N4TH
$ ,' index couples'
WRITE(6,*) ' each fifth-level file contains up to ',N5TH
$ ,' index couples'
ELSE
WRITE(6,*) ' a buffers stacking of depth',NLEVEL
$ ,' is not sufficent for '
WRITE(6,*)
$ ' this problem - you have to add more possible levels '
$
WRITE(6,*) ' to the algorithm (PARAMETER NLEVEL)'
STOP
END IF
END IF
END IF
END IF
END IF
C
NIJ=NTRIB
N1ST=INTDIV(NIJ ,NSPLIT-1)
N2ND=INTDIV(N1ST,NSPLIT-1)
N3RD=INTDIV(N2ND,NSPLIT-1)
N4TH=INTDIV(N3RD,NSPLIT-1)
N5TH=INTDIV(N4TH,NSPLIT-1)
NINDXG(1)=N1ST
NINDXG(2)=N2ND
NINDXG(3)=N3RD
NINDXG(4)=N4TH
NINDXG(5)=N5TH
CALL OPENTMP('SRTCD_',1,NFTMP)
C
C intermediate files are open
C
CO WRITE(6,*) ' NFILCD is ',NFILCD
CO IFI=1
CO IFI0=0
CO DO IGAMM=1,NBAS
CO DO IDELT=IGAMM,NBAS
CO IFI0=IFI0+1
CO IF (IFI0.GT.NFILCD) THEN
CO IFI=IFI+1
CO IFI0=0
CO END IF
CO IFI2=ITRANG(IGAMM,IDELT,NBAS)
CO WRITE(6,*) IGAMM,IDELT,' = ',IFI2,' GO TO FILE ',(IFI2-1)/NFILCD
CO $ +1
CO END DO
CO END DO
C
XNJET=0
XISTOR=0.D0
XIRD=0.D0
NBRD=0
C
WRITE(6,*)
WRITE(6,*) ' WE MIGHT HAVE ',XNMO,' MOLECULAR INTEGRALS'
WRITE(6,*)
C
IUNITR=4
CALL WOPEN(IUNITR,2)
100 CONTINUE
NBRD=NBRD+1
C--SKIPCY
CALL WREAD(LREAD)
C WRITE(6,*) ' read buffer ',NBRD
c$$$ READ(IWFILE) IBUF,JBUF,KBUF,LBUF,HBUF
c$$$ IF (IBUF(IBUFL).EQ.-1) THEN
c$$$ INTRD=JBUF(IBUFL)
c$$$ LREAD=.FALSE.
c$$$ ELSE
c$$$ INTRD=IBUFL
c$$$ END IF
C//ENDCY
CY INTRD=0
CY LREAD=.FALSE.
C
DO III=1,INTRD
XIRD=XIRD+1.D0
HHH=HBUF(III)
C
IF (ABS(HHH).LT.THRESH1) THEN
XNJET=XNJET+1.D0
ELSE
I1=IBUF(III)
I2=JBUF(III)
I3=KBUF(III)
I4=LBUF(III)
C
CALL OCAN(I1,I2,I3,I4)
IFI34=(ITRANG(I3,I4,NBAS)-1)/NFILCD+1
IFI12=(ITRANG(I1,I2,NBAS)-1)/NFILCD+1
C
C this index goes to which file?
CALL ADDTMP(IFI34,1,HHH,I1,I2,I3,I4)
CALL ADDTMP(IFI12,1,HHH,I3,I4,I1,I2)
XISTOR=XISTOR+2.D0
C
END IF
C
C close loop over integrals in one buffer
C
END DO
c$$$ IF (LREAD) GO TO 100
IF (.NOT.LREAD) GO TO 100
C
c$$$ CLOSE(IWFILE)
CALL WCLOS(2)
WRITE(6,*) ' READ ',XIRD,' BIELECTRONIC INTEGRALS '
WRITE(6,*)
C
WRITE(6,*) ' and stored them in ',NSPLIT-1,' intermediate files'
WRITE(6,*)
XUNIQ=XIRD-XNJET
WRITE(6,9221) XNJET,THRESH1,XISTOR,XUNIQ
9221 FORMAT(' ',F15.0,' Integrals were smaller than ',E12.5,/
$ ,' we will continue with ',F20.0,' AO integrals ',/
$ ,' from',F20.0,' unique integrals' )
WRITE(6,*)
C
C we close the intermediate files, but we do not delete them
C 4,1 corresponds to 4 integer indices and 1st level
C
CALL CLOSETMP(4,1,NFTMP)
WRITE(6,*)
WRITE(6,*) ' FIRST LECTURE OF BIELECTRONIC INTEGRALS COMPLETE '
WRITE(6,*)
CALL TIMING('TWOR')
C
C
CALL SRTASS('SRTCD_','AO_SRT',NBAS)
C
CALL TIMING('SRT1')
C
C open temporary files for couples ij
C
NFILIJ=INTDIV(NTRIB,NSPLIT-1)
NFTMP=INTDIV(NTRIB,NFILIJ)
N1ST=INTDIV(NTRIB ,NSPLIT-1)
N2ND=INTDIV(N1ST,NSPLIT-1)
N3RD=INTDIV(N2ND,NSPLIT-1)
N4TH=INTDIV(N3RD,NSPLIT-1)
N5TH=INTDIV(N4TH,NSPLIT-1)
C WRITE(6,*) ' NTRIB= ',NTRIB
C WRITE(6,*) ' NSPLIT= ',NSPLIT
C WRITE(6,*) ' N1ST = ',N1ST
C WRITE(6,*) ' N2ND = ',N2ND
C WRITE(6,*) ' N3RD = ',N3RD
C WRITE(6,*) ' N4TH = ',N4TH
C WRITE(6,*) ' N5TH = ',N5TH
NINDXF(1)=N1ST
NINDXF(2)=N2ND
NINDXF(3)=N3RD
NINDXF(4)=N4TH
NINDXF(5)=N5TH
WRITE(6,*)
WRITE(6,*) ' Second decomposition: we store ',NFILIJ
$ ,' couples (i,j) on one file'
WRITE(6,*) ' we need to open ',NFTMP
$ ,' intermediate files '
WRITE(6,*)
IFII =1
IFIJ =1
IFILE=1
INDX=0
INDX2=0
INDX3=0
DO I=1,IAFIN
DO J=I,IBFIN
C we have to call another function for counting the position
IFILE2=(ITRANG(I,J,IBFIN)-1)/NFILIJ+1
IF (IFILE.NE.IFILE2) THEN
WRITE(6,9821) IFILE,IFII,IFIJ,IFII2,IFIJ2,INDX
C these are the new starting indices
IFII=I
IFIJ=J
IFILE=IFILE2
INDX3=INDX3+INDX
INDX=0
END IF
INDX=INDX+1
INDX2=INDX2+1
IFII2=I
IFIJ2=J
END DO
END DO
C the last file
INDX3=INDX3+INDX
WRITE(6,9821) IFILE,IFII,IFIJ,IFII2,IFIJ2,INDX
WRITE(6,*) ' NTRIB = ',NTRIB,' INDX2 = ',INDX2,' INDX3 = ',INDX3
WRITE(6,*)
C
CALL OPENTMP('SRTIJ_',1,NFTMP)
C we open again the file of the presorted integrals
CI OPEN(UNIT=93,FILE='AO_SRT.TMP',STATUS='OLD',FORM='FORMATTED')
C--SKIPCI
OPEN(UNIT=93,FILE='AO_SRT.TMP',STATUS='OLD',FORM='UNFORMATTED')
C//ENDCI
C
DO IGAMM=1,NBAS
DO IDELT=IGAMM,NBAS
C--SKIPCI
READ(93) IGG,IDD,NBLOCK,NBDUM,NICOR
C//ENDCI
CI READ(93,*) IGG,IDD,NBLOCK,NBDUM,NICOR
IF (IGG.NE.IGAMM) THEN
WRITE(6,*) ' IGG, IGAMM: ',IGG,IGAMM
STOP ' INCONSISTENCY '
END IF
IF (IDD.NE.IDELT) THEN
WRITE(6,*) ' IDD, IDELT: ',IDD,IDELT
STOP ' INCONSISTENCY '
END IF
IF (NBUFL.NE.NBDUM) THEN
WRITE(6,*) ' NBUFL, NBDUM: ',NBUFL,NBDUM
STOP ' INCONSISTENCY '
END IF
NINTIN=NBLOCK*NBUFL+NICOR
C the first 2 quarter-transformations
CO WRITE(6,9769) IGAMM,IDELT,NBLOCK,NICOR
9769 FORMAT(' CALLING TRABIJ with IGAMM,IDELT,NBLOCK,NICOR: ',4I6)
CALL TRABIJ(93,NBLOCK,NICOR,IGAMM,IDELT,NFILIJ)
END DO
END DO
C the file 93 can be deleted
CLOSE(93,STATUS='DELETE')
C
CALL CLOSETMP(4,1,NFTMP)
C
WRITE(6,*)
WRITE(6,*) ' FIRST HALF-TRANSFORMATION COMPLETE '
WRITE(6,*)
CALL TIMING('TRF1')
CALL SRTAS2('SRTIJ_')
CALL TIMING('SRT2')
C
WRITE(6,*)
WRITE(6,*) ' SECOND HALF-TRANSFORMATION COMPLETE '
WRITE(6,*)
CALL TIMING('TRF2')
C
C we have to dump the last block of the transformed integrals
C
C we add the monoelectronic integrals
C
IZERO=0
DO I=1,NBAS
DO J=1,I
HHH=HONE(I,J)
IF (ABS(HHH).GE.THRESH3) CALL WADD(I,J,IZERO,IZERO,HHH)
END DO
END DO
C
CALL WCLOS(1)
C
WRITE(6,*)
WRITE(6,*) ' TRANSFORMED ALL BIELECTRONIC INTEGRALS '
WRITE(6,*)
RETURN
8100 CONTINUE
WRITE(6,*) ' NO FILE FOUND '
STOP
END
C
SUBROUTINE ADDTMP(IFILE,LEVEL,HHH,I1,I2,I3,I4)
INCLUDE 'param.h'
C we split files into maximally 20 smaller files
PARAMETER (NSPLIT=20,NBUL2=36000000,NLEVEL=5)
C NBUFL is the buffer length NBUL2/NSPLIT
COMMON /TWOSPL/ H0,IH0,NBUFW,NINDXB,INCORE,ISTBUF,IOFILE,NINDXF
$ ,NINDXG,NBUFL
C the buffer for integrals
DIMENSION H0(NBUL2)
C the buffer for indices
DIMENSION IH0(4,NBUL2)
C the two pointers for each file - how many buffers are on disk, how
C many elements are still in-core, we may use maximally 10 different
C levels of these pointers. ISTBUF contains the start addresses of the
C buffers
C NINDXB
DIMENSION NBUFW(NSPLIT,NLEVEL),INCORE(NSPLIT,NLEVEL)
- ,ISTBUF(NSPLIT),NINDXB(0:NLEVEL+1),IOFILE(NSPLIT)
- ,IPAIRS(NSPLIT,NLEVEL),IPAIRF(NSPLIT,NLEVEL),NINDXF(NLEVEL)
$ ,NINDXG(NLEVEL)
C.. INCLUDE 'common_twosplit.h'
C
ISTART=ISTBUF(IFILE)
INCORE(IFILE,LEVEL)=INCORE(IFILE,LEVEL)+1
IF (INCORE(IFILE,LEVEL).GT.NBUFL) THEN
IUNIT=IOFILE(IFILE)
NBUFW(IFILE,LEVEL)=NBUFW(IFILE,LEVEL)+1
C--SKIPCI
WRITE(IUNIT) ((IH0(INDX,I),INDX=1,4),H0(I),I=ISTART,ISTART
$ +NBUFL-1)
C//ENDCI
CI WRITE(IUNIT,'(4I4,F20.12)') ((IH0(INDX,I),INDX=1,4),
CI - H0(I),I=ISTART,ISTART+NBUFL-1)
INCORE(IFILE,LEVEL)=1
END IF
IADDR=ISTART+INCORE(IFILE,LEVEL)-1
H0(IADDR)=HHH
IH0(1,IADDR)=I1
IH0(2,IADDR)=I2
IH0(3,IADDR)=I3
IH0(4,IADDR)=I4
C
RETURN
END
C
SUBROUTINE CLOSETMP(NINDX,LEVEL,NFTMP)
INCLUDE 'param.h'
C we split files into maximally 20 smaller files
PARAMETER (NSPLIT=20,NBUL2=36000000,NLEVEL=5)
C NBUFL is the buffer length NBUL2/NSPLIT
COMMON /TWOSPL/ H0,IH0,NBUFW,NINDXB,INCORE,ISTBUF,IOFILE,NINDXF
$ ,NINDXG,NBUFL
C the buffer for integrals
DIMENSION H0(NBUL2)
C the buffer for indices
DIMENSION IH0(4,NBUL2)
C the two pointers for each file - how many buffers are on disk, how
C many elements are still in-core, we may use maximally 10 different
C levels of these pointers. ISTBUF contains the start addresses of the
C buffers
C NINDXB
DIMENSION NBUFW(NSPLIT,NLEVEL),INCORE(NSPLIT,NLEVEL)
- ,ISTBUF(NSPLIT),NINDXB(0:NLEVEL+1),IOFILE(NSPLIT)
- ,IPAIRS(NSPLIT,NLEVEL),IPAIRF(NSPLIT,NLEVEL),NINDXF(NLEVEL)
$ ,NINDXG(NLEVEL)
C.. INCLUDE 'common_twosplit.h'
C we close the intermediate files
DO IFIL=1,NFTMP
IUNIT =IOFILE(IFIL)
ISTART=ISTBUF(IFIL)
IEND =ISTART+INCORE(IFIL,LEVEL)-1
NTOT =NBUFW(IFIL,LEVEL)*NBUFL+INCORE(IFIL,LEVEL)
C write the remainders of the buffers to file and close the files
C--SKIPCI
WRITE(IUNIT) ((IH0(INDX,I),INDX=1,NINDX),H0(I),I=ISTART,IEND)
C//ENDCI
CI WRITE(IUNIT,'(4I4,F20.12)') ((IH0(INDX,I),INDX=1,NINDX),H0(I)
CI - ,I=ISTART,IEND)
C
CO WRITE(6,*) ' dumped buffer and CLOSED UNIT ',IUNIT
CO WRITE(6,*) ' on this unit (level',LEVEL,') we have now ',NTOT
$ ,' integrals '
CLOSE(IUNIT)
END DO
C
RETURN
END
C
SUBROUTINE OCAN(I,J,K,L)
INCLUDE 'param.h'
C put indices into canonical order
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
RETURN
END
C
C
SUBROUTINE OCAN2(I,J,K,L,WEIGHT)
IMPLICIT NONE
DOUBLE PRECISION WEIGHT
INTEGER I,J,K,L,IDUM
LOGICAL LIJ,LKL,LIK,LJL
C
C canonical ordering:
C
C i>>> WROTE ',IDREC,' RECORDS OF LENGTH',NCDMX
WRITE(6,*)
CALL TIMING('TRF1')
C
WRITE(6,*) ' second transformation '
CLOSE(IUNIT1)
OPEN(UNIT=IUNIT1,FILE='DAFILE.TMP',STATUS='OLD',
- ACCESS='DIRECT',RECL=NCDMX*8)
C
IUNITW=53
CALL WOPEN(IUNITW,1)
C
C again as many as possible in H0
C
C there are NTRIB blocks of size NTRIA in total
NABMX=NBUL2/NTRIA
WRITE(6,*) ' we have ',NABMX,' triangles of size ',NTRIA
- ,' in the main buffer '
C and NCHNK buffer loads
NCHNK=NTRIB/NABMX
WRITE(6,*) ' we have ',NCHNK,' buffer loads '
C
C how many pieces of length NCDMX do we have to read
NPIECE=NTRIB/NCDMX
NREST=NTRIB-NPIECE*NCDMX
WRITE(6,*) ' each triangle gamma/delta is cut in '
- ,NPIECE,' pieces '
WRITE(6,*) ' there remain ',NREST,' index couples '
WRITE(6,*) ' this makes ',NCDMX*NPIECE+NREST
- ,' different index couples '
C
IAB=0
IREC0=1
INDX0=0
IAST=1
IBST=1
ETWO=0.D0
IF (LOVOV) THEN
IAFIN=NOCC
IBFIN=NBAS
ELSE
IAFIN=NBAS
IBFIN=NBAS
END IF
DO IA=1,IAFIN
DO IB=IA,IBFIN
IF (IAB.EQ.NABMX) THEN
WRITE(6,*)
WRITE(6,*) ' transforming and dumping the buffer: '
WRITE(6,*) ' indices start at a=',IAST,' and b=',IBST
WRITE(6,*) ' indices end at a=',IAFI,' and b=',IBFI
C in IAST, IBST we have the beginning of the chunk, and in IAFI, IBFI
C the end
WRITE(6,*) ' H0 is full up to position ',INDX0
INDX0=0
C transformation
DO IAA=IAST,IAFI
LA=IAA.LE.NOCC
IF (IAA.EQ.IAST) THEN
IBEG=IBST
ELSE
IBEG=IAA
END IF
IF (IAA.EQ.IAFI) THEN
IBND=IBFI
ELSE
IBND=NBMO
END IF
DO IBB=IBEG,IBND
c$$$ WRITE(6,*) ' A=',IAA,' B= ',IBB
c$$$ - ,' treating elements ',INDX0+1,' to ',INDX0+NTRIA
LB=IBB.LE.NOCC
LAB=IAA.EQ.IBB
C
INDX1=INDX0
DO IGAMM=1,NBAO
DO IDELT=IGAMM,NBAO
INDX1=INDX1+1
WRKTMP(IGAMM,IDELT)=H0(INDX1)
WRKTMP(IDELT,IGAMM)=H0(INDX1)
END DO
END DO
CALL TRANS(WRKTMP,CI,NBMO,NBAO)
DO IC=IAA,IAFIN
IF (IC.EQ.IAA) THEN
IDST=IBB
ELSE
IDST=IC
END IF
DO ID=IDST,IBFIN
HHH=WRKTMP(IC,ID)
IF (ABS(HHH).GT.THRESH3.OR.LALL) THEN
CALL WADD(IAA,IBB,IC,ID,HHH)
END IF
END DO
END DO
C extract integrals for the total energy
IF (LA.AND.LB) THEN
C
IF (LAB) THEN
DO IC=IAA,NOCC
C Coulomb integrals
ETWO=ETWO+WRKTMP(IC,IC)
c$$$ WRITE(6,*) ' IA,IB,IC,IC, HHH ',IAA,IBB,IC,IC, WRKTMP(IC
c$$$ $ ,IC)
IF (IC.NE.IAA) THEN
ETWO=ETWO+3.D0*WRKTMP(IC,IC)
END IF
END DO
ELSE
C exchange integrals
ETWO=ETWO-2.D0*WRKTMP(IAA,IBB)
c$$$ WRITE(6,*) ' IA,IB,IC,ID, HHH ',IAA,IBB,IAA,IBB,WRKTMP(IAA
c$$$ $ ,IBB)
END IF
END IF
INDX0=INDX0+NTRIA
END DO
END DO
WRITE(6,*) ' H0 has been transformed until position ',INDX0
IAST=IAFI
IBST=IBFI+1
IF (IBST.EQ.NBMO+1) THEN
IAST=IAST+1
IBST=IAST
END IF
IAB=0
INDX0=0
WRITE(6,*) ' ... finished '
END IF
C
IAB=IAB+1
c$$$ WRITE(6,*) ' read triangle No',IAB,' for IA= ',IA,' and IB=',IB
c$$$ WRITE(6,*) ' elements ',INDX0+1,' to ',INDX0+NTRIA
IREC1=IREC0
DO IPI=1,NPIECE
READ(IUNIT1,REC=IREC1,ERR=8110) (H0(INDX0+I),I=1,NCDMX)
INDX0=INDX0+NCDMX
IREC1=IREC1+NTRIB
END DO
READ(IUNIT1,REC=IREC1,ERR=8110) (H0(INDX0+I),I=1,NREST)
INDX0=INDX0+NREST
IREC0=IREC0+1
IBFI=IB
IAFI=IA
END DO
END DO
C
CLOSE(IUNIT1,STATUS='DELETE')
WRITE(6,*)
WRITE(6,*) ' deleting DA-File '
WRITE(6,*)
WRITE(6,*) ' last transformation: ',INDX0/NTRIA,' triangles '
WRITE(6,*) ' indices start at a= ',IAST,' and b= ',IBST
WRITE(6,*) ' indices end at a= ',NBMO,' and b= ',NBMO
WRITE(6,*) ' H0 is full up to position ',INDX0
C
C the last chunk to be transformed
C
INDX0=0
DO IAA=IAST,IAFIN
LA=IAA.LE.NOCC
IF (IAA.EQ.IAST) THEN
IBEG=IBST
ELSE
IBEG=IAA
END IF
DO IBB=IBEG,IBFIN
LB=IBB.LE.NOCC
LAB=IAA.EQ.IBB
INDX1=INDX0
DO IGAMM=1,NBAO
DO IDELT=IGAMM,NBAO
INDX1=INDX1+1
WRKTMP(IGAMM,IDELT)=H0(INDX1)
WRKTMP(IDELT,IGAMM)=H0(INDX1)
END DO
END DO
CALL TRANS(WRKTMP,CI,NBMO,NBAO)
DO IC=IAA,IBFIN
IF (IC.EQ.IAA) THEN
IDST=IBB
ELSE
IDST=IC
END IF
DO ID=IDST,IBFIN
HHH=WRKTMP(IC,ID)
IF (ABS(HHH).GT.THRESH3.OR.LALL) THEN
CALL WADD(IAA,IBB,IC,ID,HHH)
END IF
END DO
END DO
C extract integrals for the total energy
IF (LA.AND.LB) THEN
C
IF (LAB) THEN
DO IC=1,NOCC
C Coulomb integrals
ETWO=ETWO+WRKTMP(IC,IC)
IF (IC.NE.IAA) ETWO=ETWO+3.D0*WRKTMP(IC,IC)
END DO
ELSE
C exchange integrals
ETWO=ETWO-2.D0*WRKTMP(IAA,IBB)
END IF
END IF
INDX0=INDX0+NTRIA
END DO
END DO
WRITE(6,*) ' H0 has been transformed until position ',INDX0
WRITE(6,*) ' ... finished '
WRITE(6,*)
C
CALL TIMING('TRF2')
END IF
C
C we add the monoelectronic integrals
C
IZERO=0
DO I=1,NBAS
DO J=1,I
HHH=HONE(I,J)
IF (ABS(HHH).GE.THRESH3) THEN
CALL WADD(I,J,IZERO,IZERO,HHH)
END IF
END DO
END DO
C
CALL WCLOS(1)
C
RETURN
8100 CONTINUE
WRITE(6,*) ' NO FILE FOUND '
STOP
8110 CONTINUE
WRITE(6,*) ' trying to read record No ',IREC1
WRITE(6,*) ' IREC0, IREC1, INDX0',IREC0,IREC1,INDX0
STOP
END
C
SUBROUTINE READ53(IUNIT4)
INCLUDE 'param.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
COMMON /ENERG/ EONE,ETWO,ECOUL,EEXCH
C
C.. INCLUDE 'common_syst.h'
COMMON /INTU/ S(NBASM,NBASM),F(NBASM,NBASM),HONE(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /FLOW/ THRESH1,THRESH2,THRESH3,LINFO,LOUTFO,LDIAGF
- ,LMONO,LBIEL,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE
$ ,LFROZ
C LSTH: store half-transformed integrals somewhere
LOGICAL LINFO,LOUTFO,LDIAGF,LMONO,LBIEL
$ ,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE,LFROZ
C.. INCLUDE 'flow.h'
COMMON /WBUF/ HVAL(IWBULL),IIND(IWBULL),JIND(IWBULL),KIND(IWBULL),
- LIND(IWBULL),IPOS,IUNIT,IBLE
LOGICAL LF
C
WRITE(6,*)
WRITE(6,*) ' WE READ THE INTEGRALS FROM FILE ',IUNIT4
WRITE(6,*)
C
OPEN(UNIT=26,FILE='MOTWO_FORMATTED',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
IMODE=2
CALL WOPEN(IUNIT4,IMODE)
NUMINT=0
NJET=0
100 CONTINUE
CALL WREAD(LF)
DO I=1,IPOS
I1=IIND(I)
J1=JIND(I)
K1=KIND(I)
L1=LIND(I)
HHH=HVAL(I)
IF (ABS(HHH).GE.THRESH3)
- WRITE(26,'(4I5,E24.16)') I1,J1,K1,L1,HHH
END DO
IF (.NOT.LF) GO TO 100
C
200 CONTINUE
C
CALL WCLOS(IMODE)
CLOSE(26)
C
C dump Fock matrix in molecular orbitals
OPEN(UNIT=26,FILE='FOCK.MO',STATUS='UNKNOWN',FORM
$ ='FORMATTED')
DO I=1,NBAS
DO J=1,I
HHH=F(I,J)
WRITE(26,'(2I5,E24.16)') I,J,HHH
END DO
END DO
CLOSE(26)
C
RETURN
END
SUBROUTINE AOTOMO(FILEN,CI,AMAT,NBAS)
INCLUDE 'param.h'
COMMON /FLOW/ THRESH1,THRESH2,THRESH3,LINFO,LOUTFO,LDIAGF
- ,LMONO,LBIEL,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE
$ ,LFROZ
C LSTH: store half-transformed integrals somewhere
LOGICAL LINFO,LOUTFO,LDIAGF,LMONO,LBIEL
$ ,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE,LFROZ
C.. INCLUDE 'flow.h'
DIMENSION AMAT(NBASM,NBASM)
DIMENSION CI(NBASM,NBASM)
CHARACTER*7 FILEN
DO I=1,NBAS
DO J=1,NBAS
AMAT(I,J)=0.D0
END DO
END DO
IUNIT=12
OPEN(UNIT=IUNIT,FILE=FILEN,STATUS='OLD',ERR=994,
- FORM='FORMATTED')
220 CONTINUE
READ(IUNIT,*,IOSTAT=KK) I,J,HH
IF (KK.NE.0) GO TO 221
AMAT(I,J)=HH
AMAT(J,I)=HH
GO TO 220
221 CONTINUE
CLOSE (IUNIT)
C
CALL TRANSF(AMAT,CI,NBAS)
OPEN(UNIT=17,FILE=FILEN//'.MO',FORM='FORMATTED',STATUS='UNKNOWN')
DO I=1,NBAS
DO J=I,NBAS
IF (ABS(AMAT(I,J)).GT.THRESH3) WRITE(17,'(2I6,F20.12)') I,J
$ ,AMAT(I,J)
END DO
END DO
WRITE(6,*) ' AO file ',FILEN,' transformed '
RETURN
994 CONTINUE
WRITE(6,*) ' FILE <',FILEN,'> not found '
STOP ' FILE missing '
END
FUNCTION INTDIV(I,J)
IMPLICIT NONE
INTEGER I,J,K,INTDIV
C
C integer division, which gives the next higher integer
C
C we have for instance 5/6=1 instead of zero
C
K=I/J
IF (K*J.NE.I) K=K+1
INTDIV=K
RETURN
END
SUBROUTINE TRFCOR
INCLUDE 'param.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
COMMON /ENERG/ EONE,ETWO,ECOUL,EEXCH
C
C.. INCLUDE 'common_syst.h'
COMMON /FLOW/ THRESH1,THRESH2,THRESH3,LINFO,LOUTFO,LDIAGF
- ,LMONO,LBIEL,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE
$ ,LFROZ
C LSTH: store half-transformed integrals somewhere
LOGICAL LINFO,LOUTFO,LDIAGF,LMONO,LBIEL
$ ,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE,LFROZ
C.. INCLUDE 'flow.h'
COMMON /VECTCO/ CI(NBASM,NBASM),IOCCS(NBASM),IOCC(NBASM)
C.. INCLUDE 'common_vectco.h'
C we split files into maximally 20 smaller files
PARAMETER (NSPLIT=20,NBUL2=36000000,NLEVEL=5)
C NBUFL is the buffer length NBUL2/NSPLIT
COMMON /TWOSPL/ H0,IH0,NBUFW,NINDXB,INCORE,ISTBUF,IOFILE,NINDXF
$ ,NINDXG,NBUFL
C the buffer for integrals
DIMENSION H0(NBUL2)
C the buffer for indices
DIMENSION IH0(4,NBUL2)
C the two pointers for each file - how many buffers are on disk, how
C many elements are still in-core, we may use maximally 10 different
C levels of these pointers. ISTBUF contains the start addresses of the
C buffers
C NINDXB
DIMENSION NBUFW(NSPLIT,NLEVEL),INCORE(NSPLIT,NLEVEL)
- ,ISTBUF(NSPLIT),NINDXB(0:NLEVEL+1),IOFILE(NSPLIT)
- ,IPAIRS(NSPLIT,NLEVEL),IPAIRF(NSPLIT,NLEVEL),NINDXF(NLEVEL)
$ ,NINDXG(NLEVEL)
C.. INCLUDE 'common_twosplit.h'
COMMON /WBUF/ HBUF(IWBULL),IBUF(IWBULL),JBUF(IWBULL)
$ ,KBUF(IWBULL),LBUF(IWBULL),INTRD,IWFILE,IBUFL
COMMON /CNTRS/ IWCNT,IMOCNT,IMOREJ,LREAD
LOGICAL LREAD
C.. INCLUDE 'common_two.h'
COMMON /INTU/ S(NBASM,NBASM),F(NBASM,NBASM),HONE(NBASM,NBASM)
$ ,ORBEN(NBASM)
C.. INCLUDE 'common_intu.h'
COMMON /WRK/ WRKTMP(NBASM,NBASM)
C.. INCLUDE 'common_wrk.h'
LOGICAL LCOUL,LTOTE
C
C the completely-in-core transformation
C
C
C count integrals
C
NBAO=NBAS
NBMO=NBAS
C
WRITE(6,*)
c$$$ WRITE(6,*) ' SIZE PARAMETERS '
c$$$ WRITE(6,*) ' NBAO = ',NBAO
c$$$ WRITE(6,*) ' NBMO = ',NBMO
C
NTRIA=NBAO*(NBAO+1)/2
IF (LOVOV) THEN
NTRIB = NOCC*(NOCC+1)/2 + NOCC*(NBAS-NOCC)
IAFIN=NOCC
IBFIN=NBAS
ELSE
NTRIB = NBMO*(NBMO+1)/2
IAFIN=NBAS
IBFIN=NBAS
END IF
C in AOs
NAO=NTRIA*(NTRIA+1)/2
NAOTMP=NTRIA*NTRIA
C half-transformed
NHFTR=NTRIA*NTRIB
C in MOs
NMO=NTRIB*(NTRIB+1)/2
C
C we use only the buffer NBUL2, and one temporary array for the blocks
C
NMX=MAX(NAO,NAOTMP)
NMX=MAX(NMX,NHFTR)
NMX=MAX(NMO,NMX)
WRITE(6,*)
WRITE(6,*) ' DOING THE WHOLE TRANSFORMATION IN-CORE'
WRITE(6,*)
DO I=1,NAOTMP
H0(I)=0.D0
END DO
C
IUNITR=4
CALL WOPEN(IUNITR,2)
IRD=0
NJET=0
100 CONTINUE
CALL WREAD(LREAD)
C
DO III=1,INTRD
IRD=IRD+1
HHH=HBUF(III)
C
IF (ABS(HHH).LT.THRESH1) THEN
NJET=NJET+1
ELSE
I1=IBUF(III)
I2=JBUF(III)
I3=KBUF(III)
I4=LBUF(III)
C
CALL OCAN2(I1,I2,I3,I4,WEIGHT)
C
C IDBTRI gives the position in a double triangle
C
IFI34=IDBTRI(I1,I2,I3,I4,NBAO)
IFI12=IDBTRI(I3,I4,I1,I2,NBAO)
C
H0(IFI34)=HHH*WEIGHT
H0(IFI12)=HHH*WEIGHT
C
END IF
C
C close loop over integrals in one buffer
C
END DO
IF (.NOT.LREAD) GO TO 100
C
CALL WCLOS(2)
WRITE(6,*) ' READ ',IRD,' BIELECTRONIC INTEGRALS '
WRITE(6,*) ' REJECTED ',NJET,' INTEGRALS '
WRITE(6,*)
C
C all AO integrals are in-core, in blocks of length NTRIA
C
INDX0=0
DO IALPH=1,NBAO
DO IBETA=IALPH,NBAO
INDX1=INDX0
DO IGAMM=1,NBAO
DO IDELT=IGAMM,NBAO
INDX1=INDX1+1
HDUM=H0(INDX1)
WRKTMP(IGAMM,IDELT)=HDUM
WRKTMP(IDELT,IGAMM)=HDUM
END DO
END DO
C
CALL TRANS(WRKTMP,CI,NBMO,NBAO)
C
C copy the result into H0, if necessary with some spacing
C
INDX1=INDX0
DO IC=1,IAFIN
DO ID=IC,IBFIN
INDX1=INDX1+1
H0(INDX1)=WRKTMP(IC,ID)
END DO
END DO
INDX0=INDX0+NTRIA
END DO
END DO
C
CALL TIMING('TRF1')
C
C now the second half-transformation
C
INDX0=1
DO IC=1,IAFIN
DO ID=IC,IBFIN
INDX1=INDX0
DO IALPH=1,NBAO
DO IBETA=IALPH,NBAO
HDUM=H0(INDX1)
WRKTMP(IALPH,IBETA)=HDUM
WRKTMP(IBETA,IALPH)=HDUM
INDX1=INDX1+NTRIA
END DO
END DO
C
CALL TRANS(WRKTMP,CI,NBMO,NBAO)
C
C we could store the transformed integrals right away onto file
C
INDX1=INDX0
DO IA=1,IAFIN
DO IB=IA,IBFIN
H0(INDX1)=WRKTMP(IA,IB)
INDX1=INDX1+NTRIA
END DO
END DO
C
INDX0=INDX0+1
END DO
END DO
C
IUNITW=53
CALL WOPEN(IUNITW,1)
C
INDX0=1
DO IC=1,IAFIN
DO ID=IC,IBFIN
INDX1=INDX0
DO IA=1,IAFIN
DO IB=IA,IBFIN
HDUM=H0(INDX1)
IF (ABS(HDUM).GE.THRESH3) THEN
IF (IA.GT.IC) THEN
CALL WADD(IC,ID,IA,IB,HDUM)
ELSE IF (IA.EQ.IC) THEN
IF (IB.GE.ID) CALL WADD(IC,ID,IA,IB,HDUM)
END IF
END IF
INDX1=INDX1+NTRIA
END DO
END DO
INDX0=INDX0+1
END DO
END DO
C
ETWO=0.D0
INDX0=1
DO IC=1,NOCC
DO ID=IC,IBFIN
IF (ID.LE.NOCC) THEN
IF (IC.EQ.ID) THEN
LCOUL=.TRUE.
ELSE
LCOUL=.FALSE.
END IF
INDX1=INDX0
DO IA=1,NOCC
DO IB=IA,IBFIN
IF (IB.LE.NOCC) THEN
HDUM=H0(INDX1)
C
IF (LCOUL) THEN
IF (IA.EQ.IB) THEN
c$$$ WRITE(6,*) ' IA,IB,IC,ID, HHH ',IA,IB,IC,ID, HDUM
ETWO=ETWO+HDUM
IF (IA.NE.IC) ETWO=ETWO+HDUM
END IF
C is it an exchange (IJ|IJ)?
ELSE
IF (IA.EQ.IC.AND.IB.EQ.ID) THEN
c$$$ WRITE(6,*) ' IA,IB,IC,ID, HHH ',IA,IB,IC,ID, HDUM
ETWO=ETWO-2.D0*HDUM
END IF
END IF
END IF
INDX1=INDX1+NTRIA
END DO
END DO
END IF
INDX0=INDX0+1
END DO
END DO
RETURN
END
C
SUBROUTINE SRTAS2(FILIN)
INCLUDE 'param.h'
COMMON /FLOW/ THRESH1,THRESH2,THRESH3,LINFO,LOUTFO,LDIAGF
- ,LMONO,LBIEL,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE
$ ,LFROZ
C LSTH: store half-transformed integrals somewhere
LOGICAL LINFO,LOUTFO,LDIAGF,LMONO,LBIEL
$ ,LEXOUT,LMULLI,LSTH,LDAF,LALL,LNORM,LOVOV,LDELE,LFROZ
C.. INCLUDE 'flow.h'
COMMON /VECTCO/ CI(NBASM,NBASM),IOCCS(NBASM),IOCC(NBASM)
C.. INCLUDE 'common_vectco.h'
C we split files into maximally 20 smaller files
PARAMETER (NSPLIT=20,NBUL2=36000000,NLEVEL=5)
C NBUFL is the buffer length NBUL2/NSPLIT
COMMON /TWOSPL/ H0,IH0,NBUFW,NINDXB,INCORE,ISTBUF,IOFILE,NINDXF
$ ,NINDXG,NBUFL
C the buffer for integrals
DIMENSION H0(NBUL2)
C the buffer for indices
DIMENSION IH0(4,NBUL2)
C the two pointers for each file - how many buffers are on disk, how
C many elements are still in-core, we may use maximally 10 different
C levels of these pointers. ISTBUF contains the start addresses of the
C buffers
C NINDXB
DIMENSION NBUFW(NSPLIT,NLEVEL),INCORE(NSPLIT,NLEVEL)
- ,ISTBUF(NSPLIT),NINDXB(0:NLEVEL+1),IOFILE(NSPLIT)
- ,IPAIRS(NSPLIT,NLEVEL),IPAIRF(NSPLIT,NLEVEL),NINDXF(NLEVEL)
$ ,NINDXG(NLEVEL)
C.. INCLUDE 'common_twosplit.h'
COMMON /WRK/ WRKTMP(NBASM,NBASM)
C.. INCLUDE 'common_wrk.h'
COMMON /SYST/ NBAS,NOCC,NVIRT
COMMON /ENERG/ EONE,ETWO,ECOUL,EEXCH
C
C.. INCLUDE 'common_syst.h'
COMMON /INTCNT/ XNAO,XNAOTMP,XNHFTR,XNMO,XNINTIN,XNINTFF,XIRD
$ ,XISTOR,XUNIQ
C.. INCLUDE 'common_intcount.h'
LOGICAL LCOUL,LTOTE
C
CHARACTER*1 FN(0:9)
DATA FN /'0','1','2','3','4','5','6','7','8','9'/
C
CHARACTER*6 FILIN,FILOUT,FILTMP,FILT1,FILT2
CHARACTER*5 FIL5
C
C sort and assemble
C
C the final file stores 2 indices plus the value of each integral
C
C we need the 6 characters of the input files, the name of the output
C file, the range of the indices
C
C open the file for the fully transformed integrals
C
XNINTFF=0.D0
XNINTIN=0.D0
IUNITW=53
CALL WOPEN(IUNITW,1)
C
C we use only levels 1 .. NLEVEL
C
FIL5=FILIN(1:5)
IF (LOVOV) THEN
NINDX=NOCC*(NOCC+1)/2+NOCC*(NBAS-NOCC)
IAFIN=NOCC
IBFIN=NBAS
ELSE
NINDX=NBAS*(NBAS+1)/2
IAFIN=NBAS
IBFIN=NBAS
END IF
WRITE(6,*)
WRITE(6,*) ' sorting of halftransformed integrals: '
WRITE(6,*) ' we have ',NINDX,' different indices '
NINDXB(1)=NINDX
NFIL1=INTDIV(NINDX,NSPLIT-1)
NINDXB(2)=NFIL1
NSPL1=NSPLIT-1
NFILS=(NINDX-1)/NINDXB(2)+1
CO WRITE(6,*) 'NINDXB(1),NINDXB(2),NFIL1,NSPL1,NFILS',NINDXB(1)
CO $ ,NINDXB(2),NFIL1,NSPL1,NFILS
II=1
DO I=1,NFILS
IPAIRS(I,1)=II
II=MIN(II+NINDXB(2)-1,NINDX)
IPAIRF(I,1)=II
II=II+1
C
I1=MOD(I,10)
I10=I/10
CO WRITE(6,*) ' FILE ',FILIN//FN(I10)//FN(I1)/
CO $ /'.TMP contains indices ',IPAIRS(I,1),' to ',IPAIRF(I,1)
END DO
C
IF (NFIL1.EQ.1) THEN
NDEEP=1
WRITE(6,*) ' ONE SINGLE DECOMPOSITION SUFFICIENT '
ELSE
NFIL2=INTDIV(NFIL1,NSPLIT-1)
NINDXB(3)=NFIL2
IF (NFIL2.EQ.1) THEN
NDEEP=2
WRITE(6,*) ' TWO SUCCESSIVE DECOMPOSITIONS NECESSARY '
ELSE
NFIL3=INTDIV(NFIL2,NSPLIT-1)
NINDXB(4)=NFIL3
IF (NFIL3.EQ.1) THEN
NDEEP=3
WRITE(6,*) ' THREE SUCCESSIVE DECOMPOSITIONS NECESSARY '
ELSE
NFIL4=INTDIV(NFIL3,NSPLIT-1)
NINDXB(5)=NFIL4
IF (NFIL4.EQ.1) THEN
NDEEP=4
WRITE(6,*) ' FOUR SUCCESSIVE DECOMPOSITIONS NECESSARY '
ELSE
NFIL5=INTDIV(NFIL4,NSPLIT-1)
NINDXB(6)=NFIL5
IF (NFIL5.EQ.1) THEN
NDEEP=5
WRITE(6,*) ' FIVE SUCCESSIVE DECOMPOSITIONS NECESSARY '
ELSE
WRITE(6,*) ' FIVE SUCCESSIVE DECOMPOSITIONS NOT SUFFICIENT '
WRITE(6,*) ' ADD MORE LEVELS .... '
STOP ' ADD MORE LEVELS '
END IF
END IF
END IF
END IF
END IF
WRITE(6,*)
C
DO I=0,NDEEP
WRITE(6,*) ' level ',I,' we have ',NINDXB(I+1)
$ ,' indices on one file '
END DO
C
WRITE(6,*)
CALL FLUSH(6)
C
IF (NDEEP.EQ.1) THEN
NCNT1=1
C NSPL1 we should have already
CO WRITE(6,*) ' KL K L NCNT1 NSPL1 '
CO WRITE(6,*) ' ---------------------------'
DO K=1,IAFIN
DO L=K,IBFIN
KL=ITRANG(K,L,IBFIN)
C
CO WRITE(6,9503) KL,K,L,NCNT1,NSPL1
C
C the actions:
C
CO WRITE(6,*) ' we need file ',NCNT1,' at level 1'
CO WRITE(6,*) ' r e a d i n g File No ',NCNT1,' at level 1'
C we write how many blocks of buffer length NBUFL and how many
c integrals are on the file
FILTMP=FILIN
I1 =MOD(NCNT1,10)
I10=NCNT1/10
NBLOCKS=NBUFW(NCNT1,1)
NICOR=INCORE(NCNT1,1)
C we have all segments available
IF (NICOR+NBLOCKS.NE.0) THEN
C--SKIPCI
OPEN(UNIT=92,FILE=FILTMP//FN(I10)//FN(I1)//'.TMP',STATUS
$ ='OLD',FORM='UNFORMATTED')
C//ENDCI
CI OPEN(UNIT=92,FILE=FILTMP//FN(I10)//FN(I1)//'.TMP',STATUS
CI $ ='OLD',FORM='FORMATTED')
XNINTIN=XNINTIN+NBLOCKS*NBUFL+NICOR
DO INDK=1,IAFIN
DO IDELT=1,NBAS
WRKTMP(INDK,IDELT)=0.D0
END DO
END DO
DO IBL=1,NBLOCKS
C--SKIPCI
READ(92) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NBUFL)
C//ENDCI
CI READ(92,*) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NBUFL)
C..FILE 'rdint.inc'
C here we should transform and save right away
DO III=1,NBUFL
C third quarter
IGAMM=IH0(1,III)
IDELT=IH0(2,III)
HHH=H0(III)
C we need to transform INDK from 1 to I only
DO INDK=1,K
WRKTMP(INDK,IDELT)=WRKTMP(INDK,IDELT)+CI(IGAMM,INDK)*HHH
WRKTMP(INDK,IGAMM)=WRKTMP(INDK,IGAMM)+CI(IDELT,INDK)*HHH
END DO
END DO
C.. INCLUDE 'rdint.inc'
END DO
C--SKIPCI
READ(92) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NICOR)
C//ENDCI
CI READ(92,*) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NICOR)
C..FILE 'thrfour.inc'
9736 FORMAT(' TRABKL: stored ',I8,' fully transformed integrals')
DO III=1,NICOR
C first quarter
IGAMM=IH0(1,III)
IDELT=IH0(2,III)
HHH=H0(III)
DO INDK=1,K
WRKTMP(INDK,IDELT)=WRKTMP(INDK,IDELT)+CI(IGAMM,INDK)*HHH
WRKTMP(INDK,IGAMM)=WRKTMP(INDK,IGAMM)+CI(IDELT,INDK)*HHH
END DO
END DO
C
C we need only the indices (INDK INDL) .LEQ. (IJ)
IF (K.LE.NOCC.AND.L.LE.NOCC) THEN
LTOTE=.TRUE.
IF (K.EQ.L) THEN
LCOUL=.TRUE.
ELSE
LCOUL=.FALSE.
END IF
ELSE
LTOTE=.FALSE.
END IF
C
CO NINTFF=0.D0
DO INDK=1,K
IF (K.EQ.INDK) THEN
LMAX=L
ELSE
LMAX=NBAS
END IF
DO INDL=INDK,LMAX
CUMUL=0.D0
DO IDELT=1,NBAS
HTRF=WRKTMP(INDK,IDELT)
CUMUL=CUMUL+CI(IDELT,INDL)*HTRF
END DO
C CUMUL is the completed integral
IF (ABS(CUMUL).GE.THRESH3.OR.LALL) THEN
C WRITE(6,*) ' CALLING WADD WITH ',INDK,INDL,K,L,CUMUL
CALL WADD(INDK,INDL,K,L,CUMUL)
XNINTFF=XNINTFF+1.D0
CO NINTFF=NINTFF+1
END IF
IF (LTOTE) THEN
IF (INDK.LE.NOCC.AND.INDL.LE.NOCC) THEN
IF (LCOUL.AND.INDK.EQ.INDL) THEN
IF (K.EQ.INDK) THEN
C (ii|ii) = 2 (ii|ii) - (ii|ii)
C WRITE(6,*) ' (II|II) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO+ CUMUL
ELSE
C 2 { (ii|jj) + (jj|ii) }
C WRITE(6,*) ' (II|JJ) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO+4.D0*CUMUL
END IF
END IF
IF (K.EQ.INDK.AND.L.EQ.INDL) THEN
C - { (ij|ji) + (ji|ij) }
IF (K.NE.L) THEN
C WRITE(6,*) ' (IJ|IJ) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO-2.D0*CUMUL
END IF
END IF
END IF
END IF
END DO
END DO
CO WRITE(6,9736) NINTFF
C.. INCLUDE 'thrfour.inc'
END IF
CO WRITE(6,*) ' we may delete File No ',NCNT1,' at level 1'
CLOSE(92,STATUS='DELETE')
CO WRITE(6,*) ' DELETED FILE ',FILTMP//FN(I10)//FN(I1)//'.TMP'
NCNT1=NCNT1+1
END DO
END DO
9503 FORMAT(I6,I6,I4,I7,I7)
C
C 222222222222222222222222222222222222222222222222222222222222222222222222222
C
ELSE IF (NDEEP.EQ.2) THEN
NCNT1=1
NCNT2=1
C NSPL1 we should have already
NSPL2=0
CO WRITE(6,*) ' KL K L NCNT1 NCNT2 NSPL1 NSPL2 '
CO WRITE(6,*) ' -----------------------------------------'
DO K=1,IAFIN
DO L=K,IBFIN
KL=ITRANG(K,L,IBFIN)
C
CO WRITE(6,9603) KL,K,L,NCNT1,NCNT2,NSPL1,NSPL2
C
C the actions:
C
CO WRITE(6,*) ' we need file ',NCNT2,' at level 2'
IF (NCNT2.GT.NSPL2) THEN
KLST1=IPAIRS(NCNT1,1)
KLFI1=IPAIRF(NCNT1,1)
IF (KLST1.NE.KL) THEN
WRITE(6,*) ' KL, KLST1 ',KL,KLST1
STOP ' KL .NE. KLST1 '
END IF
NSPL2=(KLFI1-KLST1)/NINDXB(3)+1
II=KLST1
DO I=1,NSPL2-1
IPAIRS(I,2)=II
II=II+NINDXB(3)
IPAIRF(I,2)=II-1
END DO
IPAIRS(NSPL2,2)=II
IPAIRF(NSPL2,2)=KLFI1
CO WRITE(6,*) ' we need file ',NCNT1,' at level 1'
CO WRITE(6,*) ' covering index pairs ',KLST1,' to ',KLFI1
CO WRITE(6,*) ' we have to decompose File No ',NCNT1
CO $ ,' at level 1 into',NSPL2,' files'
FILTMP=FIL5//FN(2)
DO I=1,NSPL2
I1=MOD(I,10)
I10=I/10
CO WRITE(6,*) ' FILE ',FILTMP//FN(I10)//FN(I1)
CO $ ,'.TMP contains indices ',IPAIRS(I,2),' to ',IPAIRF(I
CO $ ,2)
END DO
CALL DECOMF(NCNT1,NSPL2,KLST1,KLFI1,FILIN,FILTMP,1,2
$ ,NBAS)
CO WRITE(6,*) ' we may delete File No ',NCNT1,' at level 1'
NCNT1=NCNT1+1
NCNT2=1
END IF
CO WRITE(6,*) ' r e a d i n g File No ',NCNT2,' at level 2'
C we write how many blocks of buffer length NBUFL and how many
c integrals are on the file
FILTMP=FIL5//FN(2)
I1 =MOD(NCNT2,10)
I10=NCNT2/10
NBLOCKS=NBUFW(NCNT2,2)
NICOR=INCORE(NCNT2,2)
IF (NBLOCKS+NICOR.NE.0) THEN
C we have all segments available
C--SKIPCI
OPEN(UNIT=92,FILE=FILTMP//FN(I10)//FN(I1)//'.TMP',STATUS
$ ='OLD',FORM='UNFORMATTED')
C//ENDCI
CI OPEN(UNIT=92,FILE=FILTMP//FN(I10)//FN(I1)//'.TMP',STATUS
CI $ ='OLD',FORM='FORMATTED')
XNINTIN=XNINTIN+NBLOCKS*NBUFL+NICOR
DO INDK=1,IAFIN
DO IDELT=1,NBAS
WRKTMP(INDK,IDELT)=0.D0
END DO
END DO
DO IBL=1,NBLOCKS
C--SKIPCI
READ(92) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NBUFL)
C//ENDCI
CI READ(92,*) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NBUFL)
C here we should transform and save right away
DO III=1,NBUFL
C third quarter
IGAMM=IH0(1,III)
IDELT=IH0(2,III)
HHH=H0(III)
C we need to transform INDK from 1 to I only
DO INDK=1,K
WRKTMP(INDK,IDELT)=WRKTMP(INDK,IDELT)+CI(IGAMM,INDK)*HHH
WRKTMP(INDK,IGAMM)=WRKTMP(INDK,IGAMM)+CI(IDELT,INDK)*HHH
END DO
END DO
C.. INCLUDE 'rdint.inc'
END DO
C--SKIPCI
READ(92) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NICOR)
C//ENDCI
CI READ(92,*) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NICOR)
DO III=1,NICOR
C first quarter
IGAMM=IH0(1,III)
IDELT=IH0(2,III)
HHH=H0(III)
DO INDK=1,K
WRKTMP(INDK,IDELT)=WRKTMP(INDK,IDELT)+CI(IGAMM,INDK)*HHH
WRKTMP(INDK,IGAMM)=WRKTMP(INDK,IGAMM)+CI(IDELT,INDK)*HHH
END DO
END DO
C
C we need only the indices (INDK INDL) .LEQ. (IJ)
IF (K.LE.NOCC.AND.L.LE.NOCC) THEN
LTOTE=.TRUE.
IF (K.EQ.L) THEN
LCOUL=.TRUE.
ELSE
LCOUL=.FALSE.
END IF
ELSE
LTOTE=.FALSE.
END IF
C
CO NINTFF=0.D0
DO INDK=1,K
IF (K.EQ.INDK) THEN
LMAX=L
ELSE
LMAX=NBAS
END IF
DO INDL=INDK,LMAX
CUMUL=0.D0
DO IDELT=1,NBAS
HTRF=WRKTMP(INDK,IDELT)
CUMUL=CUMUL+CI(IDELT,INDL)*HTRF
END DO
C CUMUL is the completed integral
IF (ABS(CUMUL).GE.THRESH3.OR.LALL) THEN
C WRITE(6,*) ' CALLING WADD WITH ',INDK,INDL,K,L,CUMUL
CALL WADD(INDK,INDL,K,L,CUMUL)
XNINTFF=XNINTFF+1.D0
CO NINTFF=NINTFF+1
END IF
IF (LTOTE) THEN
IF (INDK.LE.NOCC.AND.INDL.LE.NOCC) THEN
IF (LCOUL.AND.INDK.EQ.INDL) THEN
IF (K.EQ.INDK) THEN
C (ii|ii) = 2 (ii|ii) - (ii|ii)
C WRITE(6,*) ' (II|II) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO+ CUMUL
ELSE
C 2 { (ii|jj) + (jj|ii) }
C WRITE(6,*) ' (II|JJ) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO+4.D0*CUMUL
END IF
END IF
IF (K.EQ.INDK.AND.L.EQ.INDL) THEN
C - { (ij|ji) + (ji|ij) }
IF (K.NE.L) THEN
C WRITE(6,*) ' (IJ|IJ) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO-2.D0*CUMUL
END IF
END IF
END IF
END IF
END DO
END DO
CO WRITE(6,9736) NINTFF
C.. INCLUDE 'thrfour.inc'
END IF
CO WRITE(6,*) ' we may delete File No ',NCNT2,' at level 2'
CLOSE(92,STATUS='DELETE')
CO WRITE(6,*) ' DELETED FILE ',FILTMP//FN(I10)//FN(I1)//'.TMP'
NCNT2=NCNT2+1
END DO
END DO
9603 FORMAT(I6,I6,I4,I7,I4,I7,I4)
C
C 3333333333333333333333333333333333333333333333333333333333333333333333333
C
ELSE IF (NDEEP.EQ.3) THEN
NCNT1=1
NCNT2=1
NCNT3=1
C NSPL1 we should have already
NSPL2=0
NSPL3=0
CO WRITE(6,*) ' KL K L NCNT1 NCNT2 NCNT3 '
CO $ ,'NSPL1 NSPL2 NSPL3 '
CO WRITE(6,*) ' -----------'
CO $ ,'-----------------------------------'
DO K=1,IAFIN
DO L=K,IBFIN
KL=ITRANG(K,L,IBFIN)
C
CO WRITE(6,9703) KL,K,L,NCNT1,NCNT2,NCNT3,NSPL1,NSPL2,NSPL3
C
C the actions:
C
CO WRITE(6,*) ' we need file ',NCNT3,' at level 3'
IF (NCNT3.GT.NSPL3) THEN
CO WRITE(6,*) ' we need file ',NCNT2,' at level 2'
IF (NCNT2.GT.NSPL2) THEN
KLST1=IPAIRS(NCNT1,1)
KLFI1=IPAIRF(NCNT1,1)
IF (KLST1.NE.KL) THEN
WRITE(6,*) ' KL, KLST1 ',KL,KLST1
STOP ' KL .NE. KLST1 '
END IF
NSPL2=(KLFI1-KLST1)/NINDXB(3)+1
II=KLST1
DO I=1,NSPL2-1
IPAIRS(I,2)=II
II=II+NINDXB(3)
IPAIRF(I,2)=II-1
END DO
IPAIRS(NSPL2,2)=II
IPAIRF(NSPL2,2)=KLFI1
CO WRITE(6,*) ' we need file ',NCNT1,' at level 1'
CO WRITE(6,*) ' covering index pairs ',KLST1,' to ',KLFI1
CO WRITE(6,*) ' we have to decompose File No ',NCNT1
CO $ ,' at level 1 into',NSPL2,' files'
FILTMP=FIL5//FN(2)
DO I=1,NSPL2
I1=MOD(I,10)
I10=I/10
CO WRITE(6,*) ' FILE ',FILTMP//FN(I10)//FN(I1)
CO $ ,'.TMP contains indices ',IPAIRS(I,2),' to ',IPAIRF(I
CO $ ,2)
END DO
CALL DECOMF(NCNT1,NSPL2,KLST1,KLFI1,FILIN,FILTMP,1,2
$ ,NBAS)
CO WRITE(6,*) ' we may delete File No ',NCNT1,' at level 1'
NCNT1=NCNT1+1
NCNT2=1
END IF
KLST2=IPAIRS(NCNT2,2)
KLFI2=IPAIRF(NCNT2,2)
IF (KLST2.NE.KL) STOP ' KL .NE. KLST2 '
NSPL3=(KLFI2-KLST2)/NINDXB(4)+1
II=KLST2
DO I=1,NSPL3-1
IPAIRS(I,3)=II
II=II+NINDXB(4)
IPAIRF(I,3)=II-1
END DO
IPAIRS(NSPL3,3)=II
IPAIRF(NSPL3,3)=KLFI2
CO WRITE(6,*) ' we need file ',NCNT2,' at level 2'
CO WRITE(6,*) ' covering index pairs ',KLST2,' to ',KLFI2
CO WRITE(6,*) ' we have to decompose File No ',NCNT2
CO $ ,' at level 2 into',NSPL3,' files'
FILTMP=FIL5//FN(3)
DO I=1,NSPL3
I1=MOD(I,10)
I10=I/10
CO WRITE(6,*) ' FILE ',FILTMP//FN(I10)//FN(I1)
CO $ ,'.TMP contains indices ',IPAIRS(I,3)
CO $ ,' to ',IPAIRF(I,3)
END DO
FILT1=FIL5//FN(2)
FILT2=FIL5//FN(3)
CALL DECOMF(NCNT2,NSPL3,KLST2,KLFI2,FILT1,FILT2,2,3
$ ,NBAS)
CO WRITE(6,*) ' we may delete File No ',NCNT2,' at level 2'
NCNT2=NCNT2+1
NCNT3=1
END IF
CO WRITE(6,*) ' r e a d i n g File No ',NCNT3,' at level 3'
C we write how many blocks of buffer length NBUFL and how many
c integrals are on the file
FILTMP=FIL5//FN(3)
I1 =MOD(NCNT3,10)
I10=NCNT3/10
NBLOCKS=NBUFW(NCNT3,3)
NICOR=INCORE(NCNT3,3)
IF (NICOR+NBLOCKS.NE.0) THEN
C we have all segments available
C--SKIPCI
OPEN(UNIT=92,FILE=FILTMP//FN(I10)//FN(I1)//'.TMP',STATUS
$ ='OLD',FORM='UNFORMATTED')
C//ENDCI
CI OPEN(UNIT=92,FILE=FILTMP//FN(I10)//FN(I1)//'.TMP',STATUS
CI $ ='OLD',FORM='FORMATTED')
XNINTIN=XNINTIN+NBLOCKS*NBUFL+NICOR
DO INDK=1,IAFIN
DO IDELT=1,NBAS
WRKTMP(INDK,IDELT)=0.D0
END DO
END DO
DO IBL=1,NBLOCKS
C--SKIPCI
READ(92) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NBUFL)
C//ENDCI
CI READ(92,*) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NBUFL)
C here we should transform and save right away
DO III=1,NBUFL
C third quarter
IGAMM=IH0(1,III)
IDELT=IH0(2,III)
HHH=H0(III)
C we need to transform INDK from 1 to I only
DO INDK=1,K
WRKTMP(INDK,IDELT)=WRKTMP(INDK,IDELT)+CI(IGAMM,INDK)*HHH
WRKTMP(INDK,IGAMM)=WRKTMP(INDK,IGAMM)+CI(IDELT,INDK)*HHH
END DO
END DO
C.. INCLUDE 'rdint.inc'
END DO
C--SKIPCI
READ(92) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NICOR)
C//ENDCI
CI READ(92,*) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NICOR)
DO III=1,NICOR
C first quarter
IGAMM=IH0(1,III)
IDELT=IH0(2,III)
HHH=H0(III)
DO INDK=1,K
WRKTMP(INDK,IDELT)=WRKTMP(INDK,IDELT)+CI(IGAMM,INDK)*HHH
WRKTMP(INDK,IGAMM)=WRKTMP(INDK,IGAMM)+CI(IDELT,INDK)*HHH
END DO
END DO
C
C we need only the indices (INDK INDL) .LEQ. (IJ)
IF (K.LE.NOCC.AND.L.LE.NOCC) THEN
LTOTE=.TRUE.
IF (K.EQ.L) THEN
LCOUL=.TRUE.
ELSE
LCOUL=.FALSE.
END IF
ELSE
LTOTE=.FALSE.
END IF
C
CO NINTFF=0.D0
DO INDK=1,K
IF (K.EQ.INDK) THEN
LMAX=L
ELSE
LMAX=NBAS
END IF
DO INDL=INDK,LMAX
CUMUL=0.D0
DO IDELT=1,NBAS
HTRF=WRKTMP(INDK,IDELT)
CUMUL=CUMUL+CI(IDELT,INDL)*HTRF
END DO
C CUMUL is the completed integral
IF (ABS(CUMUL).GE.THRESH3.OR.LALL) THEN
C WRITE(6,*) ' CALLING WADD WITH ',INDK,INDL,K,L,CUMUL
CALL WADD(INDK,INDL,K,L,CUMUL)
XNINTFF=XNINTFF+1.D0
CO NINTFF=NINTFF+1
END IF
IF (LTOTE) THEN
IF (INDK.LE.NOCC.AND.INDL.LE.NOCC) THEN
IF (LCOUL.AND.INDK.EQ.INDL) THEN
IF (K.EQ.INDK) THEN
C (ii|ii) = 2 (ii|ii) - (ii|ii)
C WRITE(6,*) ' (II|II) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO+ CUMUL
ELSE
C 2 { (ii|jj) + (jj|ii) }
C WRITE(6,*) ' (II|JJ) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO+4.D0*CUMUL
END IF
END IF
IF (K.EQ.INDK.AND.L.EQ.INDL) THEN
C - { (ij|ji) + (ji|ij) }
IF (K.NE.L) THEN
C WRITE(6,*) ' (IJ|IJ) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO-2.D0*CUMUL
END IF
END IF
END IF
END IF
END DO
END DO
CO WRITE(6,9736) NINTFF
C.. INCLUDE 'thrfour.inc'
END IF
CO WRITE(6,*) ' we may delete File No ',NCNT3,' at level 3'
CLOSE(92,STATUS='DELETE')
CO WRITE(6,*) ' DELETED FILE ',FILTMP//FN(I10)//FN(I1)//'.TMP'
NCNT3=NCNT3+1
END DO
END DO
9703 FORMAT(I6,I6,I4,I7,I4,I4,I7,I4,I4)
C
C 4444444444444444444444444444444444444444444444444444444444444444444444444444
C
ELSE IF (NDEEP.EQ.4) THEN
NCNT1=1
NCNT2=1
NCNT3=1
NCNT4=1
C NSPL1 we should have already
NSPL2=0
NSPL3=0
NSPL4=0
CO WRITE(6,*) ' KL K L NCNT1 NCNT2 NCNT3 '
CO $ ,'NCNT4 NSPL1 NSPL2 NSPL3 NSPL4 '
CO WRITE(6,*) ' -----------'
CO $ ,'--------------------------------------------'
DO K=1,IAFIN
DO L=K,IBFIN
KL=ITRANG(K,L,IBFIN)
C
CO WRITE(6,9803) KL,K,L,NCNT1,NCNT2,NCNT3,NCNT4,NSPL1,NSPL2,NSPL3
CO $ ,NSPL4
C
C the actions:
C
CO WRITE(6,*) ' we need file ',NCNT4,' at level 4'
IF (NCNT4.GT.NSPL4) THEN
CO WRITE(6,*) ' we need file ',NCNT3,' at level 3'
IF (NCNT3.GT.NSPL3) THEN
CO WRITE(6,*) ' we need file ',NCNT2,' at level 2'
IF (NCNT2.GT.NSPL2) THEN
KLST1=IPAIRS(NCNT1,1)
KLFI1=IPAIRF(NCNT1,1)
IF (KLST1.NE.KL) THEN
WRITE(6,*) ' KL, KLST1 ',KL,KLST1
STOP ' KL .NE. KLST1 '
END IF
NSPL2=(KLFI1-KLST1)/NINDXB(3)+1
II=KLST1
DO I=1,NSPL2-1
IPAIRS(I,2)=II
II=II+NINDXB(3)
IPAIRF(I,2)=II-1
END DO
IPAIRS(NSPL2,2)=II
IPAIRF(NSPL2,2)=KLFI1
CO WRITE(6,*) ' we need file ',NCNT1,' at level 1'
CO WRITE(6,*) ' covering index pairs ',KLST1,' to ',KLFI1
CO WRITE(6,*) ' we have to decompose File No ',NCNT1
CO $ ,' at level 1 into',NSPL2,' files'
FILTMP=FIL5//FN(2)
DO I=1,NSPL2
I1=MOD(I,10)
I10=I/10
CO WRITE(6,*) ' FILE ',FILTMP//FN(I10)//FN(I1)
CO $ ,'.TMP contains indices ',IPAIRS(I,2),' to ',IPAIRF(I
CO $ ,2)
END DO
CALL DECOMF(NCNT1,NSPL2,KLST1,KLFI1,FILIN,FILTMP,1,2
$ ,NBAS)
CO WRITE(6,*) ' we may delete File No ',NCNT1,' at level 1'
NCNT1=NCNT1+1
NCNT2=1
END IF
KLST2=IPAIRS(NCNT2,2)
KLFI2=IPAIRF(NCNT2,2)
IF (KLST2.NE.KL) STOP ' KL .NE. KLST2 '
NSPL3=(KLFI2-KLST2)/NINDXB(4)+1
II=KLST2
DO I=1,NSPL3-1
IPAIRS(I,3)=II
II=II+NINDXB(4)
IPAIRF(I,3)=II-1
END DO
IPAIRS(NSPL3,3)=II
IPAIRF(NSPL3,3)=KLFI2
CO WRITE(6,*) ' we need file ',NCNT2,' at level 2'
CO WRITE(6,*) ' covering index pairs ',KLST2,' to ',KLFI2
CO WRITE(6,*) ' we have to decompose File No ',NCNT2
CO $ ,' at level 2 into',NSPL3,' files'
FILTMP=FIL5//FN(3)
DO I=1,NSPL3
I1=MOD(I,10)
I10=I/10
CO WRITE(6,*) ' FILE ',FILTMP//FN(I10)//FN(I1)
CO $ ,'.TMP contains indices ',IPAIRS(I,3)
CO $ ,' to ',IPAIRF(I,3)
END DO
FILT1=FIL5//FN(2)
FILT2=FIL5//FN(3)
CALL DECOMF(NCNT2,NSPL3,KLST2,KLFI2,FILT1,FILT2,2,3
$ ,NBAS)
CO WRITE(6,*) ' we may delete File No ',NCNT2,' at level 2'
NCNT2=NCNT2+1
NCNT3=1
END IF
KLST3=IPAIRS(NCNT3,3)
KLFI3=IPAIRF(NCNT3,3)
IF (KLST3.NE.KL) STOP ' KL .NE. KLST3 '
NSPL4=(KLFI3-KLST3)/NINDXB(5)+1
II=KLST3
DO I=1,NSPL4-1
IPAIRS(I,4)=II
II=II+NINDXB(5)
IPAIRF(I,4)=II-1
END DO
IPAIRS(NSPL4,4)=II
IPAIRF(NSPL4,4)=KLFI3
CO WRITE(6,*) ' we need file ',NCNT3,' at level 3'
CO WRITE(6,*) ' covering index pairs ',KLST3,' to ',KLFI3
CO WRITE(6,*) ' we have to decompose File No ',NCNT3
CO $ ,' at level 3 into',NSPL4,' files'
FILTMP=FIL5//FN(4)
DO I=1,NSPL4
I1=MOD(I,10)
I10=I/10
CO WRITE(6,*) ' FILE ',FILTMP//FN(I10)//FN(I1)
CO $ ,'.TMP contains indices ',IPAIRS(I,4)
CO $ ,' to ',IPAIRF(I,4)
END DO
FILT1=FIL5//FN(3)
FILT2=FIL5//FN(4)
CALL DECOMF(NCNT3,NSPL4,KLST3,KLFI3,FILT1,FILT2,3,4
$ ,NBAS)
CO WRITE(6,*) ' we may delete File No ',NCNT3,' at level 3'
NCNT3=NCNT3+1
NCNT4=1
END IF
CO WRITE(6,*) ' r e a d i n g File No ',NCNT4,' at level 4'
C we write how many blocks of buffer length NBUFL and how many
c integrals are on the file
FILTMP=FIL5//FN(4)
I1 =MOD(NCNT4,10)
I10=NCNT4/10
NBLOCKS=NBUFW(NCNT4,4)
NICOR=INCORE(NCNT4,4)
IF (NBLOCKS+NICOR.NE.0) THEN
C we have all segments available
C--SKIPCI
OPEN(UNIT=92,FILE=FILTMP//FN(I10)//FN(I1)//'.TMP',STATUS
$ ='OLD',FORM='UNFORMATTED')
C//ENDCI
CI OPEN(UNIT=92,FILE=FILTMP//FN(I10)//FN(I1)//'.TMP',STATUS
CI $ ='OLD',FORM='FORMATTED')
XNINTIN=XNINTIN+NBLOCKS*NBUFL+NICOR
DO INDK=1,IAFIN
DO IDELT=1,NBAS
WRKTMP(INDK,IDELT)=0.D0
END DO
END DO
DO IBL=1,NBLOCKS
C--SKIPCI
READ(92) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NBUFL)
C//ENDCI
CI READ(92,*) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NBUFL)
C here we should transform and save right away
DO III=1,NBUFL
C third quarter
IGAMM=IH0(1,III)
IDELT=IH0(2,III)
HHH=H0(III)
C we need to transform INDK from 1 to I only
DO INDK=1,K
WRKTMP(INDK,IDELT)=WRKTMP(INDK,IDELT)+CI(IGAMM,INDK)*HHH
WRKTMP(INDK,IGAMM)=WRKTMP(INDK,IGAMM)+CI(IDELT,INDK)*HHH
END DO
END DO
C.. INCLUDE 'rdint.inc'
END DO
C--SKIPCI
READ(92) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NICOR)
C//ENDCI
CI READ(92,*) ((IH0(INDX,II),INDX=1,4),H0(II),II=1,NICOR)
DO III=1,NICOR
C first quarter
IGAMM=IH0(1,III)
IDELT=IH0(2,III)
HHH=H0(III)
DO INDK=1,K
WRKTMP(INDK,IDELT)=WRKTMP(INDK,IDELT)+CI(IGAMM,INDK)*HHH
WRKTMP(INDK,IGAMM)=WRKTMP(INDK,IGAMM)+CI(IDELT,INDK)*HHH
END DO
END DO
C
C we need only the indices (INDK INDL) .LEQ. (IJ)
IF (K.LE.NOCC.AND.L.LE.NOCC) THEN
LTOTE=.TRUE.
IF (K.EQ.L) THEN
LCOUL=.TRUE.
ELSE
LCOUL=.FALSE.
END IF
ELSE
LTOTE=.FALSE.
END IF
C
CO NINTFF=0.D0
DO INDK=1,K
IF (K.EQ.INDK) THEN
LMAX=L
ELSE
LMAX=NBAS
END IF
DO INDL=INDK,LMAX
CUMUL=0.D0
DO IDELT=1,NBAS
HTRF=WRKTMP(INDK,IDELT)
CUMUL=CUMUL+CI(IDELT,INDL)*HTRF
END DO
C CUMUL is the completed integral
IF (ABS(CUMUL).GE.THRESH3.OR.LALL) THEN
C WRITE(6,*) ' CALLING WADD WITH ',INDK,INDL,K,L,CUMUL
CALL WADD(INDK,INDL,K,L,CUMUL)
XNINTFF=XNINTFF+1.D0
CO NINTFF=NINTFF+1
END IF
IF (LTOTE) THEN
IF (INDK.LE.NOCC.AND.INDL.LE.NOCC) THEN
IF (LCOUL.AND.INDK.EQ.INDL) THEN
IF (K.EQ.INDK) THEN
C (ii|ii) = 2 (ii|ii) - (ii|ii)
C WRITE(6,*) ' (II|II) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO+ CUMUL
ELSE
C 2 { (ii|jj) + (jj|ii) }
C WRITE(6,*) ' (II|JJ) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO+4.D0*CUMUL
END IF
END IF
IF (K.EQ.INDK.AND.L.EQ.INDL) THEN
C - { (ij|ji) + (ji|ij) }
IF (K.NE.L) THEN
C WRITE(6,*) ' (IJ|IJ) ',K,L,INDK,INDL,CUMUL
ETWO=ETWO-2.D0*CUMUL
END IF
END IF
END IF
END IF
END DO
END DO
CO WRITE(6,9736) NINTFF
C.. INCLUDE 'thrfour.inc'
END IF
CO WRITE(6,*) ' we may delete File No ',NCNT4,' at level 4'
CLOSE(92,STATUS='DELETE')
CO WRITE(6,*) ' DELETED FILE ',FILTMP//FN(I10)//FN(I1)//'.TMP'
NCNT4=NCNT4+1
END DO
END DO
9803 FORMAT(I6,I6,I4,I7,I4,I4,I4,I7,I4,I4,I4)
ELSE
WRITE(6,*) ' no implementation beyond NDEEP = 4 '
END IF
C
WRITE(6,*) ' READ ',XNINTIN,' HALF-TRANSFORMED INTEGRALS '
WRITE(6,*) ' WROTE ',XNINTFF,' FULLY TRANSFORMED INTEGRALS '
C
RETURN
END