subroutine FillMP2VDM() ! In this routine, the natural orbital 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 SystemData, only: tUEG use constants, only: dp integer :: a, b, c, i, j, a2, b2, c2, i2, j2, x, y, z, w integer :: Startab, Endab, NoOcc, NoOccC, Startc, Endc, Starti, Endi, Startj, Endj real(dp) :: MP2VDMSum character(len=*), parameter :: t_r = 'FillMP2VDM' HElement_t(dp) :: HEl01, HEl02 #ifdef CMPLX_ call stop_all('FillMP2VDM', 'Natural Orbitals not implemented for complex orbitals.') #endif ! Calculating the MP2VDM (D2_ab) matrix whose eigenvectors become the ! transformation matrix. This goes in the natural orbital matrix of ! this module. The eigenvalues are the occupation numbers of the new ! orbitals. These should decrease exponentially so that when we remove ! the orbitals with small occupation numbers we should have little ! effect on the energy. ! For the MP2VDM, we always only rotate the virtual orbitals - ! the denomonator term of the above expression would be 0 if a and b ! were occupied. The orbital labels are ordered occupied then virtual ! if spatial orbitals are being used, otherwise they go occupied beta, ! virtual beta, occupied alpha, virtual alpha. This is so the alpha and ! beta spins can be diagonalised separately and we can keep track of ! which is which when the evectors are reordered and maintain spin ! symmetry. write(stdout, *) 'Filling MP2VDM nat orb matrix' call neci_flush(stdout) FillMP2VDM_Time%timer_name = 'FillMP2VDM' call set_timer(FillMP2VDM_Time, 30) do x = 1, NoSpinCyc if (x == 1) then if (tStoreSpinOrbs) then NoOcc = nOccBeta else NoOcc = NEl / 2 end if Startab = NoOcc + 1 Endab = SpatOrbs else if (x == 2) then NoOcc = nOccAlpha Startab = SpatOrbs + NoOcc + 1 Endab = NoOrbs end if ! a and b must be of the same spin to mix, so only need to run over ! both beta then both alpha. do a2 = Startab, Endab a = SymLabelList2_rot(a2) do b2 = Startab, a2 b = SymLabelList2_rot(b2) MP2VDMSum = 0.0_dp ! When a and b beta, run over both alpha and beta virtual ! for c, then both alpha and beta virtual for both i and j ! etc. do y = 1, NoSpinCyc if (y == 1) then if (tStoreSpinOrbs) then NoOccC = nOccBeta else NoOccC = NEl / 2 end if Startc = NoOccC + 1 Endc = SpatOrbs else if (y == 2) then NoOccC = nOccAlpha Startc = SpatOrbs + NoOccC + 1 Endc = NoOrbs end if do c2 = Startc, Endc c = SymLabelList2_rot(c2) do z = 1, NoSpinCyc if (z == 1) then Starti = 1 if (tStoreSpinOrbs) then Endi = nOccBeta else Endi = NEl / 2 end if else if (z == 2) then Starti = 1 + SpatOrbs Endi = SpatOrbs + nOccAlpha end if do i2 = Starti, Endi i = SymLabelList2_rot(i2) do w = 1, NoSpinCyc if (w == 1) then Startj = 1 if (tStoreSpinOrbs) then Endj = nOccBeta else Endj = NEl / 2 end if else if (w == 2) then Startj = 1 + SpatOrbs Endj = SpatOrbs + nOccAlpha end if do j2 = Startj, Endj j = SymLabelList2_rot(j2) if (tUEG) then HEl01 = get_umat_el(a, c, i, j) HEl02 = get_umat_el(b, c, i, j) MP2VDMSum = MP2VDMSum + & &(((real(HEl01, dp)) * (2.0_dp * (real(HEl02, dp)))) /& &((ARR(2 * i, 2) + ARR(2 * j, 2) - ARR(2 * a, 2) - ARR(2 * c, 2)) & &* (ARR(2 * i, 2) + ARR(2 * j, 2) - ARR(2 * b, 2) - ARR(2 * c, 2)))) HEl02 = get_umat_el(c, b, i, j) MP2VDMSum = MP2VDMSum - & &(((real(HEl01, dp)) * (real(HEl02, dp))) /& &((ARR(2 * i, 2) + ARR(2 * j, 2) - ARR(2 * a, 2) - ARR(2 * c, 2)) * & &(ARR(2 * i, 2) + ARR(2 * j, 2) - ARR(2 * c, 2) - ARR(2 * b, 2)))) else if (tStoreSpinOrbs) then if (abs(ARR(i, 2) + ARR(j, 2) - ARR(a, 2) - ARR(c, 2)) < 1.0e-12_dp) then if (abs(real(UMAT(UMatInd(a, c, i, j)), dp)) > 1.0e-12_dp) then write(stdout, *) i, j, a, c, real(UMAT(UMatInd(a, c, i, j)), dp) call stop_all(t_r, "Dividing a non-zero by zero.") end if end if MP2VDMSum = MP2VDMSum + & (((real(UMAT(UMatInd(a, c, i, j)), dp)) & * (2.0_dp * (real(UMAT(UMatInd(b, c, i, j)), dp)))) / & ((ARR(i, 2) + ARR(j, 2) - ARR(a, 2) - ARR(c, 2)) & * (ARR(i, 2) + ARR(j, 2) - ARR(b, 2) - ARR(c, 2)))) MP2VDMSum = MP2VDMSum - & (((real(UMAT(UMatInd(a, c, i, j)), dp)) & * (real(UMAT(UMatInd(c, b, i, j)), dp))) / & ((ARR(i, 2) + ARR(j, 2) - ARR(a, 2) - ARR(c, 2)) & * (ARR(i, 2) + ARR(j, 2) - ARR(c, 2) - ARR(b, 2)))) else MP2VDMSum = MP2VDMSum + & (((real(UMAT(UMatInd(a, c, i, j)), dp)) & * (2.0_dp * (real(UMAT(UMatInd(b, c, i, j)), dp)))) / & ((ARR(2 * i, 2) + ARR(2 * j, 2) - ARR(2 * a, 2) - ARR(2 * c, 2)) & * (ARR(2 * i, 2) + ARR(2 * j, 2) - ARR(2 * b, 2) - ARR(2 * c, 2)))) MP2VDMSum = MP2VDMSum - & (((real(UMAT(UMatInd(a, c, i, j)), dp)) & * (real(UMAT(UMatInd(c, b, i, j)), dp))) / & ((ARR(2 * i, 2) + ARR(2 * j, 2) - ARR(2 * a, 2) - ARR(2 * c, 2)) & * (ARR(2 * i, 2) + ARR(2 * j, 2) - ARR(2 * c, 2) - ARR(2 * b, 2)))) end if end do end do end do end do end do end do NatOrbMat(a2, b2) = MP2VDMSum NatOrbMat(b2, a2) = MP2VDMSum end do end do end do write(stdout, *) 'Finished filling MP2VDM' call halt_timer(FillMP2VDM_Time) end subroutine FillMP2VDM