SUBROUTINE ReadPropInts(nBasis, PropFile, CoreVal, OneElInts)
integer, intent(in) :: nBasis
HElement_t(dp) :: OneElInts(nBasis, nBasis)
HElement_t(dp) z
real(dp) :: CoreVal
integer :: i, j, k, l, iunit
integer :: NORB, NELEC, MS2, ISYM, SYML(1000), IUHF
integer(int64) :: ORBSYM(1000)
integer :: iSpins, ispn, SYMLZ(1000), ST, III
integer(int64) :: ZeroedInt
real(dp) :: diff, core
character(len=100) :: file_name, PropFile
logical :: TREL, UHF
real(dp), allocatable :: FOCK(:)
character(*), parameter :: t_r = 'ReadPropInts'
NAMELIST /FCI/ NORB, NELEC, MS2, ORBSYM, ISYM, IUHF, UHF, TREL, SYML, &
SYMLZ, PROPBITLEN, NPROP, ST, III, FOCK
ZeroedInt = 0
UHF = .false.
allocate(FOCK(nbasis/2), source=0.0_dp)
if (iProcIndex == 0) then
iunit = get_free_unit()
file_name = PropFile
write(stdout, *) 'Reading integral from the file:', trim(file_name)
open(iunit, FILE=file_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
!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)
core = 0.0d0
iSpins = 2
IF ((UHF .and. (.not. tROHF)) .or. tReltvy) ISPINS = 1
if (iProcIndex == 0) then
101 continue
read(iunit, *, END=199) z, i, j, k, l
! Remove integrals that are too small
if (abs(z) < UMatEps) 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 (i == 0) then
! Reading the zero-electron part of the integrals
core = real(z, dp)
else if (k == 0) then
! Reading the one-electron part of the integrals
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(OneElInts(iSpins * i - ispn + 1, iSpins * j - ispn + 1) - z)
if (abs(OneElInts(iSpins * i - ispn + 1, iSpins * j - ispn + 1)) > 0.0d-6 .and. diff > 1.0e-7_dp) then
write(stdout, *) i, j, z, OneElInts(iSpins * i - ispn + 1, iSpins * j - ispn + 1)
call Stop_All(t_R, "Error filling OneElInts - different values for same orbitals")
end if
OneElInts(iSpins * I - ispn + 1, iSpins * J - ispn + 1) = z
#ifdef CMPLX_
OneElInts(iSpins * J - ispn + 1, iSpins * I - ispn + 1) = conjg(z)
#else
OneElInts(iSpins * J - ispn + 1, iSpins * I - ispn + 1) = z
#endif
end do
else
call stop_all(t_r, 'Cannot currently deal with two-electron properties')
end if
goto 101
end if
199 continue
if (iProcIndex == 0) close(iunit)
CoreVal = core
END SUBROUTINE ReadPropInts