SUBROUTINE GETFCIBASIS(NBASISMAX, ARR, BRR, G1, LEN, TBIN)
integer, intent(in) :: LEN
integer, intent(inout) :: nBasisMax(5, *)
integer, intent(out) :: BRR(LEN)
real(dp), intent(out) :: ARR(LEN, 2)
type(BasisFN), intent(out) :: G1(LEN)
HElement_t(dp) Z
#ifndef CMPLX_
COMPLEX(dp) :: CompInt
#endif
integer(int64) IND, MASK
INTEGER I, J, K, L, I1
INTEGER ISYMNUM, ISNMAX, SYMLZ(1000), iunit
INTEGER NORB, NELEC, MS2, ISYM, ISPINS, ISPN, SYML(1000), ST, III
integer OCC(nIrreps), CLOSED(nIrreps), FROZEN(nIrreps)
integer(int64) ORBSYM(1000)
INTEGER IUHF
character(len=*), parameter :: t_r = 'GETFCIBASIS'
LOGICAL TBIN
logical :: uhf, tRel
integer :: orbsPerIrrep(nIrreps)
real(dp), allocatable :: FOCK(:)
#ifdef 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
iunit = 0
UHF = .FALSE.
PROPBITLEN = 0
NPROP = 0
IUHF = 0
TREL = .false.
SYMLZ(:) = 0
allocate(FOCK(len/2), source=0.0_dp)
IF (iProcIndex == 0) THEN
iunit = get_free_unit()
IF (TBIN) THEN
open(iunit, FILE='FCISYM', STATUS='OLD', FORM='FORMATTED')
read(iunit, FCI)
close(iunit)
open(iunit, FILE=FCIDUMP_name, STATUS='OLD', FORM='UNFORMATTED')
ELSE
open(iunit, FILE=FCIDUMP_name, STATUS='OLD', FORM='FORMATTED')
read(iunit, FCI)
end if
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
ISYMNUM = 0
ISNMAX = 0
ISPINS = 2
IF ((UHF .and. (.not. tROHF)) .or. tReltvy) ISPINS = 1
IF (tROHF) THEN
if (.not. tMolpro) then
write(stdout, *) "Reading in Spin-orbital FCIDUMP but storing as spatial orbitals."
! write(stdout,*) "*** Warning - Fock energy of orbitals will be "&
! & //"incorrect ***"
!We are reading in the symmetry of the orbitals in spin-orbitals - we need to change this to spatial orbitals
DO i = 1, NORB - 1, 2
IF (ORBSYM(i) /= ORBSYM(i + 1)) THEN
CALL Stop_All("GetFCIBASIS", "Spin-orbitals are not ordered in symmetry pairs")
end if
ORBSYM((i + 1) / 2) = ORBSYM(i)
SYML((i + 1) / 2) = SYML(i)
SYMLZ((i + 1) / 2) = SYMLZ(i)
end do
!NORB is the length in spatial orbitals, LEN is spin orbitals (both spin for UHF)
NORB = NORB / 2
DO i = NORB + 1, 1000
ORBSYM(i) = 0
SYML(i) = 0
SYMLZ(i) = 0
end do
else
!Molpro goes here. tROHF is handled exactly the same way as RHF.
!Therefore, the arrays should be the correct length anyway
write(stdout, *) "Reading in ROHF FCIDUMP, and storing as spatial orbitals."
end if
end if
! Count the number of orbs per irrep
if (any(ORBSYM(1:NORB) == 0)) then
write(stdout, *) "WARNING: Invalid ORBSYM in FCIDUMP, are you sure you know what you are doing?"
else
orbsPerIrrep = 0
do i = 1, NORB
orbsPerIrrep(ORBSYM(i)) = orbsPerIrrep(ORBSYM(i)) + ISPINS
end do
irrepOrbOffset(1) = 0
do i = 2, nIrreps
irrepOrbOffset(i) = irrepOrbOffset(i - 1) + orbsPerIrrep(i - 1)
end do
end if
!For Molpro, ISPINS should always be 2, and therefore NORB is spatial, and len is spin orbtials
IF (LEN /= ISPINS * NORB) call stop_all(t_r, 'LEN .NE. NORB in GETFCIBASIS')
G1(1:LEN) = NullBasisFn
ARR = 0.0_dp
IF (iProcIndex == 0 .and. (.not. tMolpro)) THEN
!Just read in integrals on head node. This is only trying to read in fock energies, which aren't written out by molpro anyway
IF (TBIN) THEN
if (tMolpro) call stop_all(t_r, 'UHF Bin read not functional through molpro')
IF (UHF .or. tROHF) call stop_all(t_r, 'UHF Bin read not functional')
MASK = (2**16) - 1
!IND contains all the indices in an integer(int64) - use mask of 16bit to extract them
2 read(iunit, END=99) Z, IND
L = int(iand(IND, MASK))
IND = Ishft(IND, -16)
K = int(iand(IND, MASK))
IND = Ishft(IND, -16)
J = int(iand(IND, MASK))
IND = Ishft(IND, -16)
I = int(iand(IND, MASK))
! I=Index(I)
! J=Index(J)
! K=Index(K)
! L=Index(L)
!.. Each orbital in the file corresponds to alpha and beta spinorbitals
!Fill ARR with the energy levels
IF (I /= 0 .AND. K == 0 .AND. I == J) THEN
IF (I > 1) THEN
IF (ORBSYM(I) /= ORBSYM(I - 1)) THEN
IF (ISYMNUM > ISNMAX) ISNMAX = ISYMNUM
ISYMNUM = 0
end if
end if
ISYMNUM = ISYMNUM + 1
ARR(2 * I - 1, 1) = real(Z, dp)
ARR(2 * I, 1) = real(Z, dp)
else if (I /= 0 .AND. K == 0 .AND. J == 0) THEN
ARR(2 * I - 1, 1) = real(Z, dp)
ARR(2 * I, 1) = real(Z, dp)
end if
!.. At the moment we're ignoring the core energy
IF (I /= 0) GOTO 2
ELSE !Reading in formatted FCIDUMP file
! 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...
#ifdef CMPLX_
1 if (t_complex_ints) then
read(iunit, *, END=99) Z, I, J, K, L
else
read(iunit, *, END=99) real_time_Z, I, J, K, L
Z = dcmplx(real_time_Z, 0.0_dp)
end if
#else
1 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=99) CompInt, I, J, K, L
Z = real(CompInt, dp)
if (abs(aimag(CompInt)) > 1.0e-7_dp) then
write(stdout, *) "Using the *real* neci compiled code, &
&however non-zero imaginary parts of &
&integrals found"
write(stdout, *) "If all k-points are at Gamma or BZ &
&boundary, orbitals should be able to be &
&rotated to be real"
write(stdout, *) "Check this is the case, or rerun in &
&complex mode to handle complex integrals."
call stop_all("GETFCIBASIS", "Real orbitals indicated,&
& but imaginary part of integrals larger&
& than 1.0e-7_dp")
end if
else
if (tMolpro .or. tReadFreeFormat) then
!If calling from within molpro, integrals are written out to greater precision
!Here is where we read integrals in Molcas/NECI interface
read(iunit, *, END=99) Z, I, J, K, L
else
read(iunit, '(1X,G20.12,4I3)', END=99) 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)
IF (tROHF .and. (.not. tMolpro)) THEN
!The FCIDUMP file is in spin-orbitals - we need to transfer them to spatial orbitals (unless from molpro, where already spatial).
IF (I /= 0) THEN
I1 = GTID(I)
ELSE
I1 = 0
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
I1 = I !Create new I index, since ROHF wants both spin and spatial indicies to get the fock energies right.
end if
!.. Each orbital in the file corresponds to alpha and beta spinorbitals
!Fill ARR with the energy levels
IF (I1 /= 0 .AND. K == 0 .AND. I1 == J) THEN
IF (I1 > 1) THEN
IF (ORBSYM(I1) /= ORBSYM(I1 - 1)) THEN
IF (ISYMNUM > ISNMAX) ISNMAX = ISYMNUM
ISYMNUM = 0
end if
end if
ISYMNUM = ISYMNUM + 1
!This fills the single particle energy array (ARR) with the diagonal one-electron integrals, so that if
!there are no fock energies printed out in the FCIDUMP, then we can still order the orbitals in some way.
!If the fock energies are printed out before the one-electron integrals, this will cause problems.
!These integrals must be real.
IF (ISPINS == 1) THEN
ARR(I1, 1) = real(Z, dp)
else if (tROHF) THEN
ARR(I, 1) = real(Z, dp)
ELSE
ARR(2 * I1, 1) = real(Z, dp)
ARR(2 * I1 - 1, 1) = real(Z, dp)
end if
else if (I1 /= 0 .AND. K == 0 .AND. J == 0) THEN
! write(stdout,*) I
IF (ISPINS == 1) THEN
ARR(I1, 1) = real(Z, dp)
else if (tROHF) THEN
ARR(I, 1) = real(Z, dp)
ELSE
ARR(2 * I1, 1) = real(Z, dp)
ARR(2 * I1 - 1, 1) = real(Z, dp)
end if
! DO ISPN=1,ISPINS
! ARR(ISPINS*I-ISPN+1,1)=real(Z,dp)
! end do
end if
!.. At the moment we're ignoring the core energy
IF (I1 /= 0) GOTO 1
end if
99 CONTINUE
close(iunit)
end if
if (tMolpro .and. (iProcIndex == 0)) close(iunit)
!We now need to broadcast all the information we've just read in...
CALL MPIBCast(ISNMAX, 1)
CALL MPIBCast(ISYMNUM, 1)
CALL MPIBCast(Arr, LEN * 2)
SYMMAX = 1
iMaxLz = 0
DO I = 1, NORB
DO ISPN = 1, ISPINS
BRR(ISPINS * I - ISPN + 1) = ISPINS * I - ISPN + 1
if (TwoCycleSymGens) then
G1(ISPINS * I - ISPN + 1)%Sym%s = ORBSYM(I) - 1
else
G1(ISPINS * I - ISPN + 1)%Sym%s = ORBSYM(I)
end if
IF (tFixLz) THEN
G1(ISPINS * I - ISPN + 1)%Ml = SYMLZ(I)
ELSE
G1(ISPINS * I - ISPN + 1)%Ml = 0
end if
!.. set momentum to 0
G1(ISPINS * I - ISPN + 1)%k(1) = 0
G1(ISPINS * I - ISPN + 1)%k(2) = 0
G1(ISPINS * I - ISPN + 1)%k(3) = 0
G1(ISPINS * I - ISPN + 1)%Ms = -MOD(ISPINS * I - ISPN + 1, 2) * 2 + 1
IF (SYMMAX < ORBSYM(I)) SYMMAX = int(ORBSYM(I))
IF (abs(SYMLZ(I)) > iMaxLz) iMaxLz = abs(SYMLZ(I))
end do
end do
IF (.not. tFixLz) iMaxLz = 0
if (.not. TwoCycleSymGens) then
SYMMAX = ISYM
else
! We use bit strings to store symmetry information.
! SYMMAX needs to be the smallest power of 2 greater or equal to
! the actual number of symmetry representations spanned by the basis.
SYMMAX = 2**ceiling(log(real(SYMMAX, dp)) / log(2.0_dp))
end if
IF (tFixLz) write(stdout, "(A,I3)") "Maximum Lz orbital: ", iMaxLz
write(stdout, "(A,I3)") " Maximum number of symmetries: ", SYMMAX
NBASISMAX(1, 1) = 0
NBASISMAX(1, 2) = 0
NBASISMAX(5, 1) = 0
NBASISMAX(5, 2) = SYMMAX - 1
NBASISMAX(2, 1) = 0
NBASISMAX(2, 2) = 0
END SUBROUTINE GETFCIBASIS