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