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