SUBROUTINE READFCIINT(UMAT, umat_win, NBASIS, ECORE)
integer, intent(in) :: NBASIS
real(dp), intent(out) :: ECORE
HElement_t(dp), intent(inout) :: UMAT(:)
integer(MPIArg) :: umat_win
HElement_t(dp) Z
#ifndef CMPLX_
COMPLEX(dp) :: CompInt
#endif
INTEGER(int64) :: ZeroedInt, NonZeroInt, LzDisallowed
INTEGER I, J, K, L, X, Y, iSpinType, iunit
INTEGER NORB, NELEC, MS2, ISYM, SYML(1000)
integer(int64) ORBSYM(1000)
LOGICAL LWRITE
logical :: uhf
INTEGER ISPINS, ISPN, SYMLZ(1000), ST, III !,IDI,IDJ,IDK,IDL
integer OCC(nIrreps), CLOSED(nIrreps), FROZEN(nIrreps)
INTEGER TMatSize, IUHF
integer(int64) :: UMatSize
character(len=*), parameter :: t_r = 'READFCIINT'
real(dp) :: diff
logical :: tbad, tRel
integer(int64) :: start_ind, end_ind
integer(int64), parameter :: chunk_size = 1000000
integer:: bytecount
real(dp), allocatable :: FOCK(:)
#if defined(CMPLX_)
real(dp) :: real_time_Z
#endif
NAMELIST /FCI/ NORB, NELEC, MS2, ORBSYM, OCC, CLOSED, FROZEN, &
ISYM, IUHF, UHF, TREL, SYML, SYMLZ, PROPBITLEN, NPROP, ST, III, &
FOCK
LWRITE = .FALSE.
UHF = .FALSE.
TREL = .false.
IUHF = 0
PROPBITLEN = 0
NPROP = 0
SYMLZ(:) = 0
ZeroedInt = 0
LzDisallowed = 0
NonZeroInt = 0
iunit = 0
allocate(FOCK(nbasis/2), source=0.0_dp)
IF (iProcIndex == 0) THEN
iunit = get_free_unit()
open(iunit, FILE=FCIDUMP_name, STATUS='OLD')
read(iunit, FCI)
end if
! Currently, FOCK array has been introduced to prevent crashing calculations
! It is deallocated here, as it is not used in the code anymore.
! If you want to use FOCK, don't forget to add MPIBCast(FOCK).
deallocate(FOCK)
! Re-order the orbitals symmetry labels if required
call reorder_sym_labels(ORBSYM, SYML, SYMLZ)
!Now broadcast these values to the other processors (the values are only read in on root)
CALL MPIBCast(NORB, 1)
CALL MPIBCast(NELEC, 1)
CALL MPIBCast(MS2, 1)
CALL MPIBCast(ORBSYM, 1000)
CALL MPIBCast(SYML, 1000)
CALL MPIBCast(SYMLZ, 1000)
CALL MPIBCast(ISYM, 1)
CALL MPIBCast(IUHF, 1)
CALL MPIBCast(UHF)
CALL MPIBCast(TREL)
CALL MPIBCast(PROPBITLEN, 1)
CALL MPIBCast(NPROP, 3)
! If PropBitLen has been set then assume we're not using an Abelian
! symmetry group which has two cycle generators (ie the group has
! complex representations).
TwoCycleSymGens = PropBitLen == 0
tReltvy = tRel
ISPINS = 2
IF (UHF .and. (.not. tROHF)) ISPINS = 1
if (tMolpro .and. tUHF) ISPINS = 1
if (tReltvy) ISPINS = 1
IF (tROHF .and. (.not. tMolpro)) THEN
!We are reading in the symmetry of the orbitals in spin-orbitals - we need to change this to spatial orbitals
!NORB is the length in spatial orbitals, LEN is spin orbitals (both spin for UHF)
NORB = NORB / 2
end if
IF (iProcIndex == 0) THEN
write (stdout, '("Two-electron integrals with a magnitude over ", &
&g16.7," are screened")') UMatEps
if (tMolpro .and. tUHF) then
!In molpro, UHF FCIDUMPs are written out as:
!1: aaaa
!2: bbbb
!3: aabb
!4: aa
!5: bb
!with delimiters of 0.000000 0 0 0 0
iSpinType = 1
else
iSpinType = 0
end if
TMAT2D(:, :) = (0.0_dp)
! Can't use * as need to be backward compatible with existing
! FCIDUMP files, some of which have more than 100 basis
! functions and so the integer labels run into each other.
! This means it won't work with more than 999 basis
! functions...
#if defined(CMPLX_)
101 if (t_complex_ints) then
read(iunit, *, END=199) Z, I, J, K, L
else
read(iunit, *, END=199) real_time_Z, I, J, K, L
Z = dcmplx(real_time_Z, 0.0_dp)
end if
#else
101 CONTINUE
!It is possible that the FCIDUMP can be written out in complex notation, but still only
!have real orbitals. This occurs with solid systems, where all kpoints are at the
!gamma point or BZ boundary and have been appropriately rotated. In this case, all imaginary
!components should be zero to numerical precision.
if (tRotatedOrbsReal) then
!Still need to read in integrals as complex numbers - this FCIDUMP will be from VASP.
read(iunit, *, END=199) CompInt, I, J, K, L
Z = real(CompInt, dp)
if (abs(aimag(CompInt)) > 1.0e-7_dp) then
call stop_all("READFCIINT", "Real orbitals indicated, but imaginary part of integrals larger than 1.0e-7_dp")
end if
else
if (tMolpro .or. tReadFreeFormat) then
read(iunit, *, END=199) Z, I, J, K, L
else
read(iunit, '(1X,G20.12,4I3)', END=199) Z, I, J, K, L
end if
end if
#endif
! If a permutation is loaded, apply it to the read indices
call reorder_orb_label(I)
call reorder_orb_label(J)
call reorder_orb_label(K)
call reorder_orb_label(L)
! Remove integrals that are too small
! NOTE we need I > 0 in case we have a UHF-type FCIDUMP, as the spin
! delimiter has Z == 0
if (abs(Z) < UMatEps .and. I > 0) then
if (ZeroedInt < 100) then
write(stdout, '(a,2i4,a,2i4,a)', advance='no') &
'Ignoring integral (chem. notation) (', i, j, '|', k, &
l, '): '
write(stdout, *) Z
else if (ZeroedInt == 100) then
write(stdout, *) 'Ignored more than 100 integrals.'
write(stdout, *) 'Further threshold truncations not reported explicitly'
end if
ZeroedInt = ZeroedInt + 1
goto 101
end if
! If we are fixing Lz symmetry, test if symmetry-zero elements
! are being included
if (tFixLz) then
tbad = .false.
if (i /= 0 .and. j /= 0 .and. k /= 0 .and. l /= 0) then
if (SymLz(i) + symLz(k) /= SymLz(j) + SymLz(l)) &
tbad = .true.
end if
if (i /= 0 .and. j /= 0 .and. k == 0 .and. l == 0) then
if (SymLz(i) /= SymLz(j)) &
tbad = .true.
end if
if (tbad) then
if (LzDisallowed < 100) then
write(stdout, '(a,2i4,a,2i4,a)', advance='no') &
'Ignoring Lz disallowed integral (chem. notation)&
& (', i, j, '|', k, l, '): '
write(stdout, *) Z
else if (LzDisallowed == 100) then
write(stdout, *) 'Ignored more than 100 integrals.'
write(stdout, *) 'Further threshold truncations not reported explicitly'
end if
LzDisallowed = LzDisallowed + 1
goto 101
end if
end if
IF (tROHF .and. (.not. tMolpro)) THEN
!The FCIDUMP file is in spin-orbitals - we need to transfer them to spatial orbitals (unless molpro).
IF (I /= 0) THEN
I = GTID(I)
end if
IF (J /= 0) THEN
J = GTID(J)
end if
IF (K /= 0) THEN
K = GTID(K)
end if
IF (L /= 0) THEN
L = GTID(L)
end if
else if (tMolpro .and. tUHF) then
if (i /= 0) then
!Need to transfer the orbital index into spin orbital notation
if ((iSpinType == 1) .or. (iSpinType == 4)) then
!aaaa/aa00
I = I * 2 - 1
J = J * 2 - 1
if (iSpinType == 1) K = K * 2 - 1 !(just so it doesn't give -1!)
if (iSpinType == 1) L = L * 2 - 1
else if ((iSpinType == 2) .or. (iSpinType == 5)) then
!bbbb/bb00
I = I * 2
J = J * 2
K = K * 2
L = L * 2
else if (iSpinType == 3) then
!aabb spin type (remember its still in chemical notation!)
I = I * 2 - 1
J = J * 2 - 1
K = K * 2
L = L * 2
end if
end if
end if
!.. Each orbital in the file corresponds to alpha and beta spinorbitalsa
IF (I == 0) THEN
if (tMolpro .and. tUHF .and. (iSpinType /= 6)) then
if (abs(real(z, dp)) > 1.0e-8_dp) then
call stop_all(t_r, "This is not a delimiter in the FCIDUMP reading")
end if
iSpinType = iSpinType + 1
else
!.. Core energy
ECORE = real(Z, dp)
end if
else if (J == 0) THEN
!C.. HF Eigenvalues
! ARR(I*2-1,2)=real(Z,dp)
! ARR(I*2,2)=real(Z,dp)
! ARR(BRR(I*2-1),1)=real(Z,dp)
! ARR(BRR(I*2),1)=real(Z,dp)
! LWRITE=.TRUE.
else if (K == 0) THEN
!.. 1-e integrals
!.. These are stored as spinorbitals (with elements between different spins being 0
DO ISPN = 1, ISPINS
! Have read in T_ij. Check it's consistent with T_ji
! (if T_ji has been read in).
diff = abs(TMAT2D(ISPINS * I - ISPN + 1, ISPINS * J - ISPN + 1) - Z)
if (.not. near_zero(TMAT2D(ISPINS * I - ISPN + 1, ISPINS * J - ISPN + 1)) &
.and. diff > 1.0e-7_dp &
.and. .not. t_non_hermitian_1_body) then
write(stdout, *) i, j, Z, TMAT2D(ISPINS * I - ISPN + 1, ISPINS * J - ISPN + 1)
call stop_all("ReadFCIInt", "Error filling TMAT - different values for same orbitals")
end if
if (.not. t_non_hermitian_1_body) then
TMAT2D(ISPINS * I - ISPN + 1, ISPINS * J - ISPN + 1) = Z
#ifdef CMPLX_
TMAT2D(ISPINS * J - ISPN + 1, ISPINS * I - ISPN + 1) = conjg(Z)
#else
TMAT2D(ISPINS * J - ISPN + 1, ISPINS * I - ISPN + 1) = Z
#endif
end if
end do
ELSE
!.. 2-e integrals
!.. UMAT is stored as just spatial orbitals (not spinorbitals)
!.. we're reading in (IJ|KL), but we store <..|..> which is <IK|JL>
#ifdef CMPLX_
Z = UMatConj(I, K, J, L, Z)
#endif
IF (TUMAT2D) THEN
IF (I == J .and. I == K .and. I == L) THEN
!<ii|ii>
UMAT2D(I, I) = Z
else if ((I == J .and. K == L)) THEN
!<ij|ij> - coulomb - 1st arg > 2nd arg
X = MAX(I, K)
Y = MIN(I, K)
UMAT2D(Y, X) = Z
else if (I == L .and. J == K) THEN
!<ij|ji> - exchange - 1st arg < 2nd arg
X = MIN(I, J)
Y = MAX(I, J)
UMAT2D(Y, X) = Z
else if (I == K .and. J == L) THEN
!<ii|jj> - equivalent exchange for real orbs
X = MIN(I, J)
Y = MAX(I, J)
UMAT2D(Y, X) = Z
else if (tRIIntegrals) THEN
CALL Stop_All("ReadFCIINTS", "we should not be " &
& //"reading in generic 2e integrals from " &
& //"the FCIDUMP file with ri.")
ELSE
NonZeroInt = NonZeroInt + 1
!Read in all integrals as normal.
UMAT(UMatInd(I, K, J, L)) = Z
end if
else if (tRIIntegrals) THEN
CALL Stop_All("ReadFCIInts", "TUMAT2D should be set")
ELSE
UMAT(UMatInd(I, K, J, L)) = Z
NonZeroInt = NonZeroInt + 1
end if
end if
! IF(I.NE.0) GOTO 101
GOTO 101
199 CONTINUE
close(iunit)
if (tMolpro .and. tUHF .and. (iSpinType /= 6)) then
call stop_all(t_r, "Error in reading UHF FCIDUMP from molpro")
end if
end if
!Now broadcast the data read in
CALL MPIBCast(ZeroedInt, 1)
CALL MPIBCast(LzDisallowed, 1)
CALL MPIBCast(NonZeroInt, 1)
CALL MPIBCast(ECore, 1)
!Need to find out size of TMAT before we can BCast
CALL CalcTMATSize(nBasis, TMATSize)
CALL MPIBCast(TMAT2D, TMATSize)
IF (TUMAT2D) THEN
!Broadcast TUMAT2D...
CALL MPIBCast(UMAT2D, nStates**2)
end if
IF (.not. tRIIntegrals) THEN
CALL GetUMATSize(nBasis, UMatSize)
! If we are on a 64bit system, the maximum dimensions for MPI are
! still limited by 32bit limits.
! --> We need to loop around this
start_ind = 1
end_ind = min(UMatSize, chunk_size)
do while (start_ind <= UMatSize)
!use MPI_BYTE for transfer to be independent of the data type of UMat
bytecount = int(end_ind - start_ind + 1_int64) * int(sizeof(UMat(1)))
call MPIBCast_inter_byte(UMat(start_ind), bytecount)
start_ind = end_ind + 1
end_ind = min(UMatSize, end_ind + chunk_size)
end do
!make sure the shared memory data is synchronized on all tasks
call shared_sync_mpi(umat_win)
end if
if (ZeroedInt /= 0 .and. iProcIndex == 0) then
write(stdout, *) 'Number of removed two-index integrals: ', zeroedint
end if
if (LzDisallowed /= 0 .and. iProcIndex == 0) then
write(stdout, *) 'Number of Lz disallowed two-index integrals: ', &
LzDisallowed
end if
write(stdout, *) 'Number of non-zero integrals: ', NonZeroInt
END SUBROUTINE READFCIINT