subroutine FindNatOrbitals()
! This routine simply takes a transformation matrix and rotates the integrals to produce a new FCIDUMP file.
! In one case the transformation matrix is read in from a file TRANSFORMMAT.
! In the other, the transformation matrix is calculated from the MP2 variational density matrix.
! MP2VDM = D2_ab = sum_ijc [ t_ij^ac ( 2 t_ij^bc - t_ji^bc ) ]
! Where : t_ij^ac = - < ab | ij > / ( E_a - E_i + E_b - Ej )
! Ref : J. Chem. Phys. 131, 034113 (2009) - note: in Eqn 1, the cb indices are the wrong way round (should be bc).
use NatOrbsMod, only: SetUpNatOrbLabels, FindNatOrbs, FillCoeffT1, DeallocateNatOrbs, PrintOccTable
integer :: i, a, ierr, MinReadIn, MaxReadIn, iunit
character(len=*), parameter :: this_routine = 'FindNatOrbitals'
if (tUseMP2VarDenMat) write(stdout, *) '*** Transforming the HF orbitals into the MP2 approximate natural orbitals. ***'
if (tFindCINatOrbs) then
write(stdout, *) '*** Transforming the HF orbitals into approximate natural orbitals'
write(stdout, *) 'based on the one-electron density matrix found from the wavefunction calculated above. ***'
end if
if (tSpinOrbs) then
if (.not. tStoreSpinOrbs) then
write(stdout, *) "We want to use spin orbitals - turning on tStoreSpinOrbs."
tStoreSpinOrbs = .true.
end if
end if
if (tROHF .and. tStoreSpinOrbs) call Stop_All(this_routine, "Cannot compress open shell systems into spatial " &
& //"orbitals when rotating, turn off ROHF.")
if (tTruncRODump .and. (.not. tTruncDumpbyVal)) then
NoFrozenVirt = NoTruncOrbs(1)
else if (tTruncRODump) then
! If the 'number of frozen orbitals' is given as a cutoff - take NoFrozenVirt to be 0
!for all the allocation purposes - will set this later when
! we have the eigenvalues and know how many orbitals lie below it.
NoFrozenVirt = 0
TruncEval = TruncEvalues(1)
else
NoFrozenVirt = 0
end if
SpatOrbs = nBasis / 2
if (tStoreSpinOrbs) then
NoOrbs = nBasis
NoOcc = NEl
MinReadIn = 1
MaxReadIn = nBasis
if (tRotateVirtOnly) MinReadIn = NEl + 1
if (tRotateOccOnly) MaxReadIn = NEl
! If tStoreSpinOrbs ARR(:,2) is not filled, but we want to use it later, so just fill it here.
do i = 1, NoOrbs
ARR(BRR(i), 2) = ARR(i, 1)
end do
allocate(SymLabelCounts2_rot(2, 32), stat=ierr)
call LogMemAlloc('SymLabelCounts2_rot', 2 * 32, 4, this_routine, SymLabelCounts2_rotTag, ierr)
SymLabelCounts2_rot(:, :) = 0
! first 8 refer to the occupied, and the second to the virtual beta spin.
! third and fourth to the occupied and virtual alpha spin.
else
NoOrbs = SpatOrbs
NoOcc = NEl / 2
MinReadIn = 1
MaxReadIn = SpatOrbs
if (tRotateVirtOnly) MinReadIn = (NEl / 2) + 1
if (tRotateOccOnly) MaxReadIn = NEl / 2
allocate(SymLabelCounts2_rot(2, 16), stat=ierr)
call LogMemAlloc('SymLabelCounts2_rot', 2 * 16, 4, this_routine, SymLabelCounts2_rotTag, ierr)
SymLabelCounts2_rot(:, :) = 0
! first 8 refer to the occupied, and the second to the virtual.
end if
NoRotOrbs = NoOrbs
call ApproxMemReq()
allocate(SymLabelList2_rot(NoOrbs), stat=ierr)
call LogMemAlloc('SymLabelList2_rot', NoOrbs, 4, this_routine, SymLabelList2_rotTag, ierr)
SymLabelList2_rot(:) = 0
allocate(SymLabelList3_rot(NoOrbs), stat=ierr)
call LogMemAlloc('SymLabelList3_rot', NoOrbs, 4, this_routine, SymLabelList3_rotTag, ierr)
SymLabelList3_rot(:) = 0
allocate(SymLabelListInv_rot(NoOrbs), stat=ierr)
call LogMemAlloc('SymLabelListInv_rot', NoOrbs, 4, this_routine, SymLabelListInv_rotTag, ierr)
SymLabelListInv_rot(:) = 0
if (tReadInCoeff .or. tUseHFOrbs) then
! No symmetry, so no reordering of the orbitals - symlabellist just goes from 1-NoOrbs.
! When we are just reading in the coefficients and transforming, it does not matter about the ordering of the orbitals.
do i = 1, NoOrbs
SymLabelList2_rot(i) = i
SymLabelListInv_rot(i) = i
end do
else if (tFindCINatOrbs .or. tUseMP2VarDenMat) then
call SetupNatOrbLabels()
end if
! Yet another labelling system, SymLabelList3_rot is created here.
! This indicates the label of the transformed orbital.
! In the case where we are truncating the space, the transformed orbitals are ordered according
!to the size of the eigenvalues of the MP2VDM
! matrix when it is diagonalised. We wish to keep them in this order when transforming
!the integrals etc, so that when we truncate the last
! NoFrozenVirt orbitals, we are removing those with the smallest MP2VDM eigenvalues (occupation numbers).
! In the case where no truncation is made however, SymLabelList3_rot is the same as SymLabelList2_rot,
!so that the indexes remain the same as previously.
! This allows for the option of going straight into a spawning calc from the rotation, which is not
!possible when a truncation is performed
! because of the messed up indices.
if (tTruncRODump) then
if (MOD(NoFrozenVirt, 2) /= 0) call Stop_All(this_routine, "Must freeze virtual spin orbitals in pairs of 2.")
if (tStoreSpinOrbs) then
NoRotOrbs = NoOrbs - NoFrozenVirt
else
NoFrozenVirt = NoFrozenVirt / 2
NoRotOrbs = NoOrbs - NoFrozenVirt
end if
do i = 1, NoOrbs
SymLabelList3_rot(i) = i
end do
else
do i = 1, NoOrbs
SymLabelList3_rot(i) = SymLabelList2_rot(i)
end do
end if
allocate(CoeffT1(NoOrbs, NoRotOrbs), stat=ierr)
call LogMemAlloc(this_routine, NoRotOrbs * NoOrbs, 8, this_routine, CoeffT1Tag, ierr)
CoeffT1(:, :) = 0.0_dp
if (tSeparateOccVirt) then
do i = 1, NoRotOrbs
CoeffT1(i, i) = 1.0_dp
end do
end if
if (tUEG) then
call FindNatOrbs()
call FillCoeffT1()
else
if (tReadInCoeff) then
write(stdout, '(A)') " Reading in the transformation matrix from TRANSFORMMAT, and using this to rotate the HF orbitals."
iunit = get_free_unit()
open(iunit, file='TRANSFORMMAT', status='old')
do i = 1, NoOrbs
do a = 1, NoOrbs
read(iunit, *) CoeffT1(a, i)
end do
end do
close(iunit)
else if (tFindCINatOrbs .or. tUseMP2VarDenMat .or. tUseHFOrbs) then
if (.not. tUseHFOrbs) call FindNatOrbs()
if (tUseHFOrbs) then
call PrintOccTable()
else
call FillCoeffT1()
end if
end if
if (tPrintRODump) then
allocate(FourIndInts(NoOrbs, NoOrbs, NoOrbs, NoOrbs), stat=ierr)
call LogMemAlloc('FourIndInts', (NoOrbs**4), 8, this_routine, FourIndIntsTag, ierr)
! Then, transform2ElInts
write(stdout, *) 'Transforming the four index integrals'
call Transform2ElIntsMemSave()
write(stdout, *) 'Re-calculating the fock matrix'
call CalcFOCKMatrix()
write(stdout, *) 'Refilling the UMAT and TMAT2D'
! The ROFCIDUMP is also printed out in here.
call RefillUMATandTMAT2D()
call neci_flush(stdout)
if ((tFindCINatOrbs .or. tUseMP2VarDenMat) .and. (NoDumpTruncs > 1)) call ReTruncROFciDump()
if ((.not. tUseHFOrbs) .and. (.not. tReadInCoeff)) call DeallocateNatOrbs()
end if
if (tWriteTransMat) call WriteTransformMat()
! If a truncation is being made, the new basis will not be in the correct energetic ordering - this does not matter, as we
! never go straight into a spawning and they will be reordered when the ROFCIDUMP file is read in again.
call writebasis(stdout, G1, nBasis, ARR, BRR)
deallocate(CoeffT1)
call LogMemDeAlloc(this_routine, CoeffT1Tag)
deallocate(SymLabelList2_rot)
call LogMemDeAlloc(this_routine, SymLabelList2_rotTag)
deallocate(SymLabelListInv_rot)
call LogMemDeAlloc(this_routine, SymLabelListInv_rotTag)
if (tPrintRODump) then
deallocate(FourIndInts)
call LogMemDeAlloc(this_routine, FourIndIntsTag)
end if
end if
end subroutine FindNatOrbitals