subroutine FillOneRDM()
! Det is the number of determinants in FCIDets.
! FCIDets contains the list of all determinants in the system in bit
! string representation, FCIDets(0:NIfTot,1:Det)
! ICILevel is the max excitation level of the calculation - as in
! EXCITE ICILevel.
! FCIDetIndex(1:NEl) contains the index of FCIDets where each
! excitation level starts.
! As in FCIDetIndex(1) = 2 always I think - Excitation level 1 starts
! at the second determinant (after HF).
! Pretty sure FCIDetIndex always goes from 1:NEl even from truncated
! excite calculations.
! The elements of AllHistogram correspond to the rows of FCIDets - i.e
! to each determinant in the system. AllHistogram contains the final
! (normalised) amplitude of the determinant - with sign.
use DetCalcData, only: Det, FCIDets, FCIDetIndex, ICILevel
use DetBitOps, only: FindBitExcitLevel
use bit_reps, only: decode_bit_det
use hist_data, only: AllHistogram
integer :: excit, i, j
integer :: Starti, Endi, Startj, Endj, ExcitLevel, Ex(2, 1)
integer :: Orbi, Orbj, nJ(NEl), Orbk, k, nI(NEl), MaxExcit
integer :: Spins
logical :: tSign
real(dp) :: SignDet
#ifdef DEBUG_
character(*), parameter :: this_routine = "FillOneRDM"
#endif
! Density matrix = D_pq = < Psi | a_q+ a_p | Psi >
! = sum_ij [ c_i* c_j < D_i | a_q+ a_p | D_j > ]
! Where a_p is the annihilation, and a_q+ the creation operators.
! In other words, < D_i | a_q+ a_p | D_j > will only be non-zero if
! D_i and D_j are connected by an annihilation in p and a creation in q.
! Want to run through all determinants D_i in the final wavefunction.
! For each, find all determinants, D_j that are connected to D_i by a
! single excitation - i.e. those which differ by just one orbital. Only
! need to consider those of the same excitation level or one above or
! one below. Find the differing orbitals - these will be p and q. Sum
! in the occupations of D_i and D_j (c_i x c_j) to the matrix element
! D_pq. Take, for instance, p always =< q.
! Will get the orbitals in the original labelling system - convert it
! to this system.
! Get a list of the wavefunction with amplitudes in order of excitation.
! Depending on the type of reduced density matrix want to:
! Run through the determinants with excitation level one less, the same
! and one more.
FillOneRDM_Time%timer_name = 'FillOneRDM'
call set_timer(FillOneRDM_Time, 30)
write(stdout, *) '*** The weight of the HF determinant is : ', AllHistogram(1, 1)
write(stdout, *) 'Beginning to fill the one-electron reduced density matrix.'
if (ICILevel == 0) then
MaxExcit = NEl
else
MaxExcit = ICILevel
end if
! Run through all determinants D_i, in the final wavefunction, Psi.
! If this is done by excitation block, we then don't have to check
! the excitation level of the determinant each time.
do excit = 0, MaxExcit
! The HF only involves 'occupied' orbitals - these are not required
! if only rotating virt.
if (tRotateVirtOnly .and. tSeparateOccVirt .and. (excit == 0)) cycle
! This next bit is a bit messy because there is no row in
! FCIDetIndex for the HF - there is probably a tidier way to
! achieve the same thing, but it does the trick for now.
if (excit == 0) then ! i is the HF det.
Starti = 1
Endi = 1
Startj = 1
Endj = min((FCIDetIndex(2) - 1), Det) ! If i is the HF det, just run over singly excited j.
else if (excit == MaxExcit) then
Starti = FCIDetIndex(excit)
Endi = Det
Startj = FCIDetIndex(excit - 1)
Endj = Det
else
Starti = FCIDetIndex(excit)
Endi = FCIDetIndex(excit + 1) - 1
if (excit == 1) then
Startj = 1
if (NEl < 3) then
Endj = Det
else
Endj = FCIDetIndex(3) - 1
end if
else if (excit == (MaxExcit - 1)) then
Startj = FCIDetIndex(excit - 1)
Endj = Det
else
Startj = FCIDetIndex(excit - 1)
Endj = FCIDetIndex(excit + 2) - 1
end if
end if
! Then run through the determinants in that excitation level.
do i = Starti, Endi
! Run through all determinants D_j, with the potential to be
! connected to i by a single excitation, i.e from one
! excitation lower to one excitation higher.
do j = Startj, i
if ((i > Det) .or. (j > Det)) then
call stop_all('FillOneRDM', 'Running through i or j larger than the number of determinants.')
end if
ExcitLevel = FindBitExcitLevel(FCIDets(:, i), FCIDets(:, j), 2)
! Need to find the excitation level between D_i and D_j.
! If this is 1 - go on to add their contributions to the
! OneRDM.
if (ExcitLevel == 1) then
Ex(:, :) = 0
Ex(1, 1) = ExcitLevel
ASSERT(.not. t_3_body_excits)
call GetBitExcitation(FCIDets(:, i), FCIDets(:, j), Ex, tSign)
! Gives the orbitals involved in the excitation Ex(1,1)
! in i -> Ex(2,1) in j (in spin orbitals).
if (tStoreSpinOrbs) then
! OneRDM will be in spin orbitals - simply add the
! orbitals involved.
Orbi = SymLabelListInv_rot(Ex(1, 1))
Orbj = SymLabelListInv_rot(Ex(2, 1))
Spins = 1
else
Orbi = SymLabelListInv_rot(ceiling(real(Ex(1, 1), dp) / 2.0_dp))
Orbj = SymLabelListInv_rot(ceiling(real(Ex(2, 1), dp) / 2.0_dp))
Spins = 2
end if
if (tSign) then
SignDet = -1.0_dp
else
SignDet = 1.0_dp
end if
NatOrbMat(Orbi, Orbj) = NatOrbMat(Orbi, Orbj) + (SignDet * AllHistogram(1, i) * AllHistogram(1, j))
NatOrbMat(Orbj, Orbi) = NatOrbMat(Orbj, Orbi) + (SignDet * AllHistogram(1, i) * AllHistogram(1, j))
if ((abs(AllHistogram(1, i) * AllHistogram(1, j)) > 1.0e-12_dp) .and. &
(int(G1(SymLabelList2_rot(Orbi) * Spins)%sym%S, 4) /= &
int(G1(SymLabelList2_rot(Orbj) * Spins)%sym%S, 4))) then
write(stdout, *) 'ERROR in symmetries'
write(stdout, *) 'Ex,', Ex(1, 1), Ex(2, 1)
write(stdout, *) ceiling(real(Ex(1, 1) / 2.0_dp, dp)), ceiling(real(Ex(2, 1) / 2.0_dp, dp))
write(stdout, *) 'Orbi,', Orbi, 'Orbj,', Orbj
write(stdout, *) 'Sym(Orbi)', int(G1(SymLabelList2_rot(Orbi) * Spins)%sym%S, 4), 'Sym(Orbj)', &
int(G1(SymLabelList2_rot(Orbj) * Spins)%sym%S, 4)
call decode_bit_det(nI, FCIDets(0:NIfTot, i))
write(stdout, *) 'i', nI
call decode_bit_det(nJ, FCIDets(0:NIfTot, j))
write(stdout, *) 'j', nJ
write(stdout, *) 'AllHistogram(1,i)', AllHistogram(1, i)
write(stdout, *) 'AllHistogram(1,j)', AllHistogram(1, j)
call neci_flush(stdout)
call stop_all('FillOneRDM', 'Non-zero element between different symmetries.')
end if
else if (ExcitLevel == 0) then
call Decode_Bit_Det(nJ, FCIDets(0:NIfTot, j))
do k = 1, NEl
if (tStoreSpinOrbs) then
Orbk = SymLabelListInv_rot(nJ(k))
else
Orbk = SymLabelListInv_rot(ceiling(real(nJ(k), dp) / 2.0_dp))
end if
NatOrbMat(Orbk, Orbk) = NatOrbMat(Orbk, Orbk) + (AllHistogram(1, j)**2)
! 0.5 x because this will be added twice since we
! are not currently restricting i<k or anything.
end do
end if
end do
end do
end do
write(stdout, *) 'Filled OneRDM'
call halt_timer(FillOneRDM_Time)
end subroutine FillOneRDM