Subroutine DetPreFreezeInit()
integer ierr, ms, iEl, flagAlpha
integer i, j, Lz, OrbOrder(8, 2), FDetTemp(NEl), lmsMax
type(BasisFn) s
logical :: tGenFDet
HElement_t(dp) :: OrbE
character(25), parameter :: this_routine = 'DetPreFreezeInit'
allocate (FDet(nEl), stat=ierr)
LogAlloc(ierr, 'FDet', nEl, 4, tagFDet)
IF (tDefineDet) THEN
write (stdout, *) 'Defining FDet according to input'
do i = 1, NEl
FDet(i) = DefDet(i)
end do
! A quick check that we have defined a reasonable det.
ms = sum(get_spin_pn(fdet(1:nel)))
tRef_Not_HF = .true.
else
tGenFDet = .true.
lmsMax = sum(nOccOrbs) - sum(nClosedOrbs)
if ((sum(nOccOrbs) + sum(nClosedOrbs)) == nel &
.and. (.not. TSPN .or. abs(LMS) == lmsMax)) then
tGenFDet = .false.
if (LMS < 0) then
flagAlpha = 0
else
flagAlpha = 1
end if
iEl = 1
do i = 1, nIrreps
! nClosedOrbs is the number of alpha/beta orbs (the ones with minority spin)
do j = 1, nClosedOrbs(i)
FDet(iEl) = irrepOrbOffset(i) + 2 * j - flagAlpha
iEl = iEl + 1
end do
! nOccOrbs is the number of majority spin orbs (per irrep)
do j = 1, nOccOrbs(i)
FDet(iEl) = irrepOrbOffset(i) + 2 * j - (1 - flagAlpha)
iEl = iEl + 1
end do
end do
call sort(FDet)
end if
if (tGenFDet) then
CALL GENFDET(FDET)
IF (tUEGSpecifyMomentum) THEN
write (stdout, *) 'Defining FDet according to a momentum input'
CALL ModifyMomentum(FDET)
end if
tRef_Not_HF = .false.
end if
end if
write (stdout, "(A)", advance='no') " Fermi det (D0):"
call write_det(stdout, FDET, .true.)
Call GetSym(FDet, nEl, G1, nBasisMax, s)
write (stdout, "(A)", advance='no') " Symmetry: "
Call writesym(stdout, s%Sym, .true.)
IF (tFixLz) THEN
Call GetLz(FDet, nEl, Lz)
write (stdout, "(A,I5)") "Lz of Fermi det:", Lz
end if
CALL NECI_ICOPY(NEL, FDET, 1, NUHFDET, 1)
if (tMolpro) then
!Orbitals are ordered by occupation number from MOLPRO, and not reordered in NECI
!Therefore, we know HF determinant is first four occupied orbitals.
write (stdout, "(A)") "NECI called from MOLPRO, so assuming orbitals ordered by occupation number."
if (.not. tDefineDet) then
FDetTemp(:) = FDet(:)
else
!We have defined our own reference determinant, but still use the first orbitals for the calculation
!of 'orbital energies'
CALL GENFDET(FDETTEMP)
end if
write (stdout, "(A)") "Calculating orbital energies..."
do i = 1, nBasis
OrbE = CalcFockOrbEnergy(i, FDetTemp)
Arr(i, 1) = real(OrbE, dp)
Brr(i) = i
end do
write (stdout, "(A)") "Reordering basis by orbital energies..."
OrbOrder(:, :) = 0
call ORDERBASIS(NBASIS, Arr, Brr, OrbOrder, nBasisMax, G1)
!However, we reorder them here
call writebasis(stdout, G1, nBasis, Arr, Brr)
end if
E0HFDET = ECORE
DO I = 1, NEL
E0HFDET = E0HFDET + ARR(NUHFDET(i), 2)
end do
write (stdout, *) "Fock operator energy:", E0HFDET
! Store the value of Ms for use in other areas
calculated_ms = sum(get_spin_pn(fdet(1:nel)))
End Subroutine DetPreFreezeInit