call stop_all(t_r,”Calculation of X matrix not the same as via stability matrix”)
call stop_all(t_r,”Calculation of Y matrix not the same as via stability matrix”)
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=dp), | intent(out) | :: | Weight | |||
| real(kind=dp), | intent(out) | :: | Energy |
subroutine RunRPA_QBA(Weight, Energy)
implicit none
integer :: ierr, i, j, m, n, ex(2, maxExcit), ex2(2, maxExcit), mi_ind, nj_ind
integer :: StabilitySize, lWork, info, i_p, m_p, v, mp_ip_ind, ic
integer :: nJ(NEl), exflag, mu, id(2, maxExcit), a, i_ind, iunit, ia_ind
integer(n_int) :: iLutHF(0:NIfTot)
real(dp), intent(out) :: Weight, Energy
real(dp) :: Energy_stab, Temp_real, norm, Energy2, H0tmp, Fii
real(dp) :: X_norm, Y_norm
HElement_t(dp) :: HDiagTemp, hel, hel1, hel2
logical :: tAllExcitsFound, tParity
real(dp), allocatable :: Stability(:, :), temp2(:, :), W2(:)
real(dp), allocatable :: W(:), Work(:), S_half(:, :), temp(:, :)
real(dp), allocatable :: X_stab(:, :), Y_stab(:, :), StabilityCopy(:, :)
real(dp), allocatable :: X_chol(:, :), Y_chol(:, :)
real(dp), allocatable :: AminB(:, :), AplusB(:, :), temp3(:, :)
character(len=*), parameter :: t_r = "RunRPA_QBA"
type(ExcitGenSessionType) :: session
write(stdout, "(A)")
write(stdout, "(A)") "**************************************"
if (tDirectRPA) then
write(stdout, "(A)") "* Entering DIRECT RPA calculation *"
else
write(stdout, "(A)") "* Entering FULL RPA calculation *"
end if
write(stdout, "(A)") "**************************************"
write(stdout, "(A)")
call neci_flush(stdout)
HDiagTemp = get_helement(fDet, fDet, 0)
Energy = real(HDiagTemp, dp)
write(stdout, "(A,G25.10)") "Reference energy is: ", Energy
!Quickly find correlation energy from MP2 for comparison.
Temp_real = 0.0_dp
tAllExcitsFound = .false.
exflag = 3
ex(:, :) = 0
call EncodeBitDet(FDet, iLutHF)
HDiagTemp = GetH0Element3(FDet)
Fii = real(HDiagTemp, dp)
do while (.true.)
if (tReltvy) then
call GenExcitations4(session, FDet, nJ, exFlag, Ex, tParity, tAllExcitsFound, .false.)
else
call GenExcitations3(FDet, iLutHF, nJ, exflag, Ex, tParity, tAllExcitsFound, .false.)
end if
if (tAllExcitsFound) exit !All excits found
if (Ex(1, 2) == 0) then
ic = 1
else
ic = 2
end if
hel = get_helement(FDet, nJ, ic, Ex, tParity)
H0tmp = getH0Element3(nJ)
H0tmp = Fii - H0tmp
Temp_real = Temp_real + (hel**2) / H0tmp
end do
write(stdout, "(A,G25.10)") "For comparison, MP2 correlation energy is: ", Temp_real
Weight = (0.0_dp)
Energy = (0.0_dp)
ov_space = (nBasis - NEl) * NEl
virt_start = NEl + 1
write(stdout, "(A,I8)") "1p-1h space is: ", ov_space
allocate(A_mat(ov_space, ov_space), stat=ierr)
allocate(B_mat(ov_space, ov_space), stat=ierr)
if (ierr /= 0) call stop_all(t_r, "alloc Err")
A_mat(:, :) = 0.0_dp
B_mat(:, :) = 0.0_dp
!First, construct A and B, which correspond to:
! <HF| [a*_i a_m [H, a*_n a_j]] |HF> = A
! -<HF| [a*_i a_m [H, a*_j a_n]] |HF> = B
do j = 1, nel
ex(1, 2) = Brr(j) !Second index in integral
ex2(2, 2) = Brr(j)
do n = virt_start, nBasis
nj_ind = ov_space_ind(n, j)
ex(2, 2) = Brr(n) !fourth index in integral
ex2(1, 2) = Brr(n)
do i = 1, nel
ex(2, 1) = Brr(i) !Third index in integral
ex2(2, 1) = Brr(i)
do m = virt_start, nBasis
mi_ind = ov_space_ind(m, i)
ex(1, 1) = Brr(m) !First index in integral
ex2(1, 1) = Brr(m)
if (tDirectRPA) then
!No exchange interactions
! Obtain spatial rather than spin indices if required
id = gtID(ex)
! Only non-zero contributions if Ms preserved in each term (consider
! physical notation).
if ((G1(ex(1, 1))%Ms == G1(ex(2, 1))%Ms) .and. &
(G1(ex(1, 2))%Ms == G1(ex(2, 2))%Ms)) then
hel1 = get_umat_el(id(1, 1), id(1, 2), id(2, 1), &
id(2, 2))
else
hel1 = (0)
end if
A_mat(mi_ind, nj_ind) = real(hel1, dp)
!Now for B matrix
id = gtID(ex2)
! Only non-zero contributions if Ms preserved in each term (consider
! physical notation).
if ((G1(ex2(1, 1))%Ms == G1(ex2(2, 1))%Ms) .and. &
(G1(ex2(1, 2))%Ms == G1(ex2(2, 2))%Ms)) then
hel1 = get_umat_el(id(1, 1), id(1, 2), id(2, 1), &
id(2, 2))
else
hel1 = (0)
end if
B_mat(mi_ind, nj_ind) = real(hel1, dp)
else
!Full antisymmetrized integrals
HEl1 = sltcnd_excit(nJ, Excite_2_t(ex), .false.)
HEl2 = sltcnd_excit(nJ, Excite_2_t(ex2), .false.)
A_mat(mi_ind, nj_ind) = real(HEl1, dp)
B_mat(mi_ind, nj_ind) = real(HEl2, dp)
end if
end do
end do
end do
end do
!Now add the diagonal part to A
do i = 1, nel
do m = virt_start, nBasis
mi_ind = ov_space_ind(m, i)
A_mat(mi_ind, mi_ind) = A_mat(mi_ind, mi_ind) + (Arr(Brr(m), 2) - Arr(Brr(i), 2))
! write(stdout,*) Arr(Brr(m),2)-Arr(Brr(i),2)
end do
end do
!Check that A is hermitian and B is symmetric
do i = 1, ov_space
do j = 1, ov_space
if (abs(B_mat(i, j) - B_mat(j, i)) > 1.0e-7_dp) then
write(stdout, *) i, j, B_mat(i, j), B_mat(j, i), abs(B_mat(i, j) - B_mat(j, i))
call stop_all(t_r, "B not symmetric")
end if
if (abs(A_mat(i, j) - A_mat(j, i)) > 1.0e-7_dp) then
write(stdout, *) i, j, A_mat(i, j), A_mat(j, i), abs(A_mat(i, j) - A_mat(j, i))
call stop_all(t_r, "A not hermitian")
end if
end do
end do
if (tStabilityAnalysis) then
!Calculate stability analysis of HF solution, and from this, calculate RPA amplitudes
!This should give identical results to other method, will be quite a bit slower, but also
!gives information about the stability of the HF solution in all directions.
write(stdout, "(A)")
write(stdout, "(A)") "Calculating RPA from stability matrix..."
#ifdef CMPLX_
call stop_all(t_r, "Not coded up for complex integrals. Bug ghb24")
#endif
!Stability = ( A B )
! ( B* A* )
StabilitySize = 2 * ov_space
allocate(Stability(StabilitySize, StabilitySize), stat=ierr)
Stability(:, :) = 0.0_dp
Stability(1:ov_space, 1:ov_space) = A_mat(1:ov_space, 1:ov_space)
!Assume all integrals real to start with
! Stability(ov_space+1:StabilitySize,1:ov_space)=conjg(B(1:ov_space,1:ov_space))
Stability(ov_space + 1:StabilitySize, 1:ov_space) = B_mat(1:ov_space, 1:ov_space)
Stability(1:ov_space, ov_space + 1:StabilitySize) = B_mat(1:ov_space, 1:ov_space)
!Assume all integrals real to start with
! Stability(ov_space+1:StabilitySize,ov_space+1:StabilitySize)=conjg(A_mat(1:ov_space,1:ov_space))
Stability(ov_space + 1:StabilitySize, ov_space + 1:StabilitySize) = A_mat(1:ov_space, 1:ov_space)
! call writematrix(Stability,'Stability matrix')
!Now diagonalize
!Find optimal space
allocate(StabilityCopy(StabilitySize, StabilitySize))
StabilityCopy(:, :) = Stability(:, :)
allocate(W(StabilitySize), stat=ierr) !Eigenvalues of stability matrix
allocate(Work(1))
if (ierr /= 0) call stop_all(t_r, "alloc err")
W(:) = 0.0_dp
lWork = -1
info = 0
call dsyev('V', 'U', StabilitySize, Stability, StabilitySize, W, Work, lWork, info)
if (info /= 0) call stop_all(t_r, 'workspace query failed')
lwork = int(work(1)) + 1
deallocate(work)
allocate(work(lwork))
call dsyev('V', 'U', StabilitySize, Stability, StabilitySize, W, Work, lwork, info)
if (info /= 0) call stop_all(t_r, "Diag failed")
deallocate(work)
do i = 1, StabilitySize
if (W(i) < 0.0_dp) then
write(stdout, *) i, W(i)
call stop_all(t_r, "HF solution not stable. Not local minimum. Recompute HF.")
end if
end do
if (.not. tdirectRPA) then
write(stdout, "(A)") "Stability matrix positive definite. HF solution is minimum. RPA stable"
end if
!Now compute S^(1/2), and transform into original basis
allocate(S_half(StabilitySize, StabilitySize))
S_half(:, :) = 0.0_dp
do i = 1, StabilitySize
S_half(i, i) = sqrt(W(i))
end do
allocate(temp(StabilitySize, StabilitySize))
call dgemm('n', 'n', StabilitySize, StabilitySize, StabilitySize, 1.0_dp, Stability, StabilitySize, &
S_half, StabilitySize, 0.0_dp, temp, StabilitySize)
call dgemm('n', 't', StabilitySize, StabilitySize, StabilitySize, 1.0_dp, temp, StabilitySize, Stability, &
StabilitySize, 0.0_dp, S_half, StabilitySize)
!S_half is now S^1/2 in the original basis
!Check this by squaring it.
call dgemm('n', 't', StabilitySize, StabilitySize, StabilitySize, 1.0_dp, S_half, StabilitySize, S_half, &
StabilitySize, 0.0_dp, temp, StabilitySize)
do i = 1, StabilitySize
do j = 1, StabilitySize
if (abs(StabilityCopy(i, j) - temp(i, j)) > 1.0e-7_dp) then
call stop_all(t_r, 'S^1/2 not calculated correctly in original basis')
end if
end do
end do
temp(:, :) = 0.0_dp
do i = 1, ov_space
temp(i, i) = 1.0_dp
end do
do i = ov_space + 1, StabilitySize
temp(i, i) = -1.0_dp
end do
allocate(temp2(StabilitySize, StabilitySize))
call dgemm('n', 'n', StabilitySize, StabilitySize, StabilitySize, 1.0_dp, S_half, StabilitySize, temp, &
StabilitySize, 0.0_dp, temp2, StabilitySize)
call dgemm('n', 'n', StabilitySize, StabilitySize, StabilitySize, 1.0_dp, temp2, StabilitySize, S_half, &
StabilitySize, 0.0_dp, temp, StabilitySize)
!Now diagonalize temp = S^(1/2) (1 0 \\ 0 -1 ) S^(1/2)
lWork = -1
allocate(W2(StabilitySize))
allocate(Work(1))
W2(:) = 0.0_dp
call dsyev('V', 'U', StabilitySize, temp, StabilitySize, W2, Work, lWork, info)
if (info /= 0) call stop_all(t_r, 'workspace query failed')
lwork = int(work(1)) + 1
deallocate(work)
allocate(work(lwork))
call dsyev('V', 'U', StabilitySize, temp, StabilitySize, W2, Work, lwork, info)
if (info /= 0) call stop_all(t_r, "Diag failed")
deallocate(work)
! call writevector(W2,'Excitation energies')
! temp now holds the eigenvectors X~ Y~
! W2 runs over StabilitySize eigenvalues (ov_space*2). Therefor we expect redundant pairs of +-W2, corresponding
! to pairs of eigenvectors (X^v Y^v) and (X^v* Y^v*) (Same in real spaces).
do i = 1, ov_space
!This they are listed in order of increasing eigenvalue, we should be able to easily check that they pair up
if (abs(W2(i) + W2(StabilitySize - i + 1)) > 1.0e-7_dp) then
write(stdout, *) i, StabilitySize - i + 1, W2(i), W2(StabilitySize - i + 1), abs(W2(i) - W2(StabilitySize - i + 1))
call stop_all(t_r, "Excitation energy eigenvalues do not pair")
end if
end do
!We actually have everything we need for the energy already now. However, calculate X and Y too.
!Now construct (X Y) = S^(-1/2) (X~ Y~)
!First get S^(-1/2) in the original basis
S_half(:, :) = 0.0_dp
do i = 1, StabilitySize
S_half(i, i) = -sqrt(W(i))
end do
call dgemm('n', 'n', StabilitySize, StabilitySize, StabilitySize, 1.0_dp, Stability, StabilitySize, S_half, &
StabilitySize, 0.0_dp, temp2, StabilitySize)
call dgemm('n', 't', StabilitySize, StabilitySize, StabilitySize, 1.0_dp, temp2, StabilitySize, Stability, &
StabilitySize, 0.0_dp, S_half, StabilitySize)
!S_half is now S^(-1/2) in the original basis
!Now multiply S^(-1/2) (X~ y~)
call dgemm('n', 'n', StabilitySize, StabilitySize, StabilitySize, 1.0_dp, S_half, StabilitySize, temp, &
StabilitySize, 0.0_dp, temp2, StabilitySize)
!Check that eigenvectors are also paired.
!Rotations among degenerate sets will screw this up though
! do i=1,ov_space
! write(stdout,*) "Eigenvectors: ",i,StabilitySize-i+1,W2(i),W2(StabilitySize-i+1)
! do j=1,StabilitySize
! write(stdout,*) j,temp2(j,i),temp2(j,StabilitySize-i+1)
! end do
! end do
! call writematrix(temp2,'X Y // Y X',.true.)
!temp2 should now be a matrix of (Y X)
! (X Y)
! This is the other way round to normal, but due to the fact that our eigenvalues are ordered -ve -> +ve
! TODO: Are the signs of this matrix correct?
allocate(X_stab(ov_space, ov_space)) !First index is (m,i) compound index. Second is the eigenvector index.
allocate(Y_stab(ov_space, ov_space))
X_stab(:, :) = 0.0_dp
Y_stab(:, :) = 0.0_dp
!Put the eigenvectors corresponding to *positive* eigenvalues into the X_stab and Y_stab arrays.
X_stab(1:ov_space, 1:ov_space) = temp2(1:ov_space, ov_space + 1:StabilitySize)
Y_stab(1:ov_space, 1:ov_space) = -temp2(ov_space + 1:StabilitySize, ov_space + 1:StabilitySize)
deallocate(temp2)
!Normalize the eigenvectors appropriately
do mu = 1, ov_space
norm = 0.0_dp
Y_norm = 0.0_dp
X_norm = 0.0_dp
do i = 1, ov_space
norm = norm + X_stab(i, mu) * X_stab(i, mu) - Y_stab(i, mu) * Y_stab(i, mu)
Y_norm = Y_norm + Y_stab(i, mu) * Y_stab(i, mu)
X_norm = X_norm + X_stab(i, mu) * X_stab(i, mu)
end do
if (norm <= 0.0_dp) then
write(stdout, *) "Norm^2 for vector ", mu, " is: ", norm
call stop_all(t_r, 'norm undefined')
end if
norm = sqrt(norm)
do i = 1, ov_space
X_stab(i, mu) = X_stab(i, mu) / norm
Y_stab(i, mu) = Y_stab(i, mu) / norm
end do
if (Y_norm > X_norm / 2.0_dp) then
write(stdout, *) "Warning: hole amplitudes large for excitation: ", mu, &
" Quasi-boson approximation breaking down."
write(stdout, *) "Norm of X component: ", X_norm
write(stdout, *) "Norm of Y component: ", Y_norm
end if
end do
! call writematrix(X_stab,'X',.true.)
!Now check orthogonality
call Check_XY_orthogonality(X_stab, Y_stab)
! call writevector(W2,'Stab_eigenvalues')
!Now check that we satisfy the original RPA equations
!For the *positive* eigenvalue space (since we have extracted eigenvectors corresponding to this), check that:
! ( A B ) (X) = E_v(X )
! ( B A ) (Y) (-Y)
deallocate(temp)
allocate(temp(StabilitySize, ov_space))
temp(1:ov_space, 1:ov_space) = X_stab(:, :)
temp(ov_space + 1:StabilitySize, 1:ov_space) = Y_stab(:, :)
allocate(temp2(StabilitySize, ov_space))
call dgemm('n', 'n', StabilitySize, ov_space, StabilitySize, 1.0_dp, StabilityCopy, StabilitySize, temp, &
StabilitySize, 0.0_dp, temp2, StabilitySize)
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp2(j, i) - (W2(i + ov_space) * X_stab(j, i))) > 1.0e-6_dp) then
write(stdout, *) i, j, temp2(j, i), (W2(i + ov_space) * X_stab(j, i)), W2(i + ov_space)
call stop_all(t_r, "RPA equations not satisfied for positive frequencies in X matrix")
end if
end do
end do
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp2(j + ov_space, i) - (-W2(i + ov_space) * Y_stab(j, i))) > 1.0e-6_dp) then
write(stdout, *) i, j, temp2(j + ov_space, i), (-W2(i + ov_space) * Y_stab(j, i)), -W2(i + ov_space)
call stop_all(t_r, "RPA equations not satisfied for positive frequencies in Y matrix")
end if
end do
end do
deallocate(temp, temp2)
!Is is also satisfied the other way around?
!Check that we also satisfy (still for the *positive* eigenvalues):
! ( A B ) (Y) = -E_v(Y )
! ( B A ) (X) (-X)
allocate(temp(StabilitySize, ov_space))
temp(1:ov_space, 1:ov_space) = Y_stab(:, :)
temp(ov_space + 1:StabilitySize, 1:ov_space) = X_stab(:, :)
allocate(temp2(StabilitySize, ov_space))
call dgemm('n', 'n', StabilitySize, ov_space, StabilitySize, 1.0_dp, StabilityCopy, StabilitySize, temp, &
StabilitySize, 0.0_dp, temp2, StabilitySize)
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp2(j, i) - (-W2(i + ov_space) * Y_stab(j, i))) > 1.0e-6_dp) then
call stop_all(t_r, "RPA equations not satisfied for negative frequencies in X matrix")
end if
end do
end do
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp2(j + ov_space, i) - (W2(i + ov_space) * X_stab(j, i))) > 1.0e-6_dp) then
call stop_all(t_r, "RPA equations not satisfied for negative frequencies in Y matrix")
end if
end do
end do
deallocate(temp, temp2)
!TODO: Finally, check that we satisfy eq. 1 in the Scuseria paper for X and Y defined for positive eigenvalues...
do i = ov_space + 1, StabilitySize
do j = 1, StabilitySize
StabilityCopy(i, j) = -StabilityCopy(i, j)
end do
end do
!Stability copy is now (A B // -B -A)
allocate(temp(StabilitySize, ov_space))
allocate(temp2(StabilitySize, ov_space))
temp = 0.0_dp
temp(1:ov_space, 1:ov_space) = X_stab(:, :)
temp(ov_space + 1:StabilitySize, 1:ov_space) = Y_stab(:, :)
call dgemm('n', 'n', StabilitySize, ov_space, StabilitySize, 1.0_dp, StabilityCopy, StabilitySize, temp, &
StabilitySize, 0.0_dp, temp2, StabilitySize)
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp2(j, i) - (W2(i + ov_space) * X_stab(j, i))) > 1.0e-7_dp) then
call stop_all(t_r, "RPA equations not satisfied for X")
end if
end do
end do
do i = 1, ov_space
do j = ov_space + 1, StabilitySize
if (abs(temp2(j, i) - (W2(i + ov_space) * Y_stab(j - ov_space, i))) > 1.0e-7_dp) then
call stop_all(t_r, "RPA equations not satisfied for Y")
end if
end do
end do
deallocate(temp, temp2)
!Now calculate energy, in two different ways:
!1. -1/2 Tr[A] + 1/2 sum_v E_v(positive)
Energy_stab = 0.0_dp
do i = 1, ov_space
Energy_stab = Energy_stab + W2(ov_space + i) - A_mat(i, i)
end do
Energy_stab = Energy_stab / 2.0_dp
if (tDirectRPA) then
write(stdout, "(A,G25.10)") "Direct RPA energy from stability analysis (plasmonic RPA-TDA excitation energies): ", &
Energy_stab
else
write(stdout, "(A,G25.10)") "Full RPA energy from stability analysis (plasmonic RPA-TDA excitation energies): ", &
Energy_stab
end if
Energy_stab = 0.0_dp
!E = 0.25 * Tr[BZ] where Z = Y X^-1
allocate(temp2(ov_space, ov_space))
temp2(:, :) = 0.0_dp
!Find X^-1
call d_inv(X_stab, temp2)
! call writematrix(temp2,'X^-1',.true.)
allocate(temp(ov_space, ov_space))
!Find Z (temp)
call dgemm('n', 'n', ov_space, ov_space, ov_space, 1.0_dp, Y_stab, ov_space, temp2, ov_space, 0.0_dp, temp, ov_space)
!Find BZ (temp2)
call dgemm('n', 'n', ov_space, ov_space, ov_space, 1.0_dp, B_mat, ov_space, temp, ov_space, 0.0_dp, temp2, ov_space)
!Take trace of BZ
do i = 1, ov_space
Energy_stab = Energy_stab + temp2(i, i)
end do
Energy_stab = Energy_stab / 2.0_dp
if (tDirectRPA) then
write(stdout, "(A,G25.10)") "Direct RPA energy from stability analysis (Ring-CCD: 1/2 Tr[BZ]): ", Energy_stab
else
write(stdout, "(A,G25.10)") "Full RPA energy from stability analysis (Ring-CCD: 1/2 Tr[BZ]): ", Energy_stab
end if
Energy_stab = 0.0_dp
do i = 1, ov_space
Y_norm = 0.0_dp
do j = 1, ov_space
Y_norm = Y_norm + Y_stab(j, i)**2
end do
Energy_stab = Energy_stab - W2(i + ov_space) * Y_norm
end do
if (tDirectRPA) then
write(stdout, "(A,G25.10)") "Direct RPA energy from stability analysis (Y-matrix): ", Energy_stab
else
write(stdout, "(A,G25.10)") "Full RPA energy from stability analysis (Y-matrix): ", Energy_stab
end if
deallocate(W2, W, temp, temp2, StabilityCopy, Stability)
end if
write(stdout, *)
write(stdout, "(A)") "Calculating RPA via Cholesky decomposition"
!Now, we calculate the RPA amplitudes via consideration of the boson operators.
!This will generally be quicker, and should give the same result.
!Calculate A-B
allocate(AminB(ov_space, ov_space))
do i = 1, ov_space
do j = 1, ov_space
AminB(j, i) = A_mat(j, i) - B_mat(j, i)
end do
end do
if (.true.) then
!Optional
!Check that AminB is positive definite - required for Cholesky decomposition
allocate(temp(ov_space, ov_space))
temp(:, :) = AminB(:, :)
!Also check it is symmetric!
do i = 1, ov_space
do j = 1, ov_space
if (abs(AminB(i, j) - AminB(j, i)) > 1.0e-7_dp) then
call stop_all(t_r, "A - B matrix is not symmetric")
end if
end do
end do
!Now diagonalise
lWork = -1
allocate(W(ov_space))
allocate(Work(1))
W(:) = 0.0_dp
call dsyev('V', 'U', ov_space, temp, ov_space, W, Work, lWork, info)
if (info /= 0) call stop_all(t_r, 'workspace query failed')
lwork = int(work(1)) + 1
deallocate(work)
allocate(work(lwork))
call dsyev('V', 'U', ov_space, temp, ov_space, W, Work, lWork, info)
if (info /= 0) call stop_all(t_r, "Diag failed")
deallocate(work)
do i = 1, ov_space
if (W(i) < 0.0_dp) then
call writevector(W, 'AminB eigenvalues')
call stop_all(t_r, "A-B not positive definite - cannot do cholesky decomposition")
end if
end do
deallocate(temp, W)
end if
!Construct triangular matrices via cholesky decomposition
call dpotrf('U', ov_space, AminB, ov_space, info)
if (info /= 0) call stop_all(t_r, 'Cholesky decomposition failed.')
!Now set lower triangular part to zero
do i = 1, ov_space
do j = 1, i - 1
AminB(i, j) = 0.0_dp
end do
end do
if (.true.) then
!Optional
! call writematrix(AminB,'T')
!Check cholesky decomposition successful
allocate(temp(ov_space, ov_space))
call dgemm('T', 'N', ov_space, ov_space, ov_space, 1.0_dp, AminB, ov_space, AminB, ov_space, 0.0_dp, temp, ov_space)
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp(i, j) - (A_mat(i, j) - B_mat(i, j))) > 1.0e-7_dp) then
call stop_all(t_r, 'Cholesky decomposition not as expected')
end if
end do
end do
deallocate(temp)
end if
allocate(AplusB(ov_space, ov_space), stat=ierr)
do i = 1, ov_space
do j = 1, ov_space
AplusB(j, i) = A_mat(j, i) + B_mat(j, i)
end do
end do
allocate(temp(ov_space, ov_space))
call dgemm('N', 'N', ov_space, ov_space, ov_space, 1.0_dp, AminB, ov_space, AplusB, ov_space, 0.0_dp, temp, ov_space)
call dgemm('N', 'T', ov_space, ov_space, ov_space, 1.0_dp, temp, ov_space, AminB, ov_space, 0.0_dp, AplusB, ov_space)
deallocate(temp)
!Check T (A + B) T^T is actually symmetric
do i = 1, ov_space
do j = 1, ov_space
if (abs(AplusB(i, j) - AplusB(j, i)) > 1.0e-7_dp) then
call stop_all(t_r, 'Not symmetric')
end if
end do
end do
!AplusB is now actually T (A + B) T^T
!Diagonalize this
lWork = -1
allocate(W(ov_space))
allocate(Work(1))
W(:) = 0.0_dp
call dsyev('V', 'U', ov_space, AplusB, ov_space, W, Work, lWork, info)
if (info /= 0) call stop_all(t_r, 'workspace query failed')
lwork = int(work(1)) + 1
deallocate(work)
allocate(work(lwork))
call dsyev('V', 'U', ov_space, AplusB, ov_space, W, Work, lWork, info)
if (info /= 0) call stop_all(t_r, "Diag failed")
deallocate(work)
do i = 1, ov_space
if (W(i) < 0.0_dp) then
call stop_all(t_r, "Excitation energies^2 not positive. Error in diagonalization.")
end if
end do
! call writevector(sqrt(W),'Eigenvalues')
!Excitation energies are now the +- sqrt of eigenvalues
!Construct W^(-1/2) T^T R (temp)
!and W^(1/2) T^-1 R (temp2)
allocate(temp(ov_space, ov_space))
call dgemm('T', 'N', ov_space, ov_space, ov_space, 1.0_dp, AminB, ov_space, AplusB, ov_space, 0.0_dp, temp, ov_space)
do i = 1, ov_space
!Mulitply each column by corresponding eigenvalue^(-1/2). Only deal with positive eigenvalues here
do j = 1, ov_space
temp(j, i) = temp(j, i) * (1.0_dp / (sqrt(sqrt(W(i)))))
end do
end do
allocate(temp2(ov_space, ov_space))
allocate(temp3(ov_space, ov_space)) !Temp3 will hold the inverse of the T matrix temporarily.
call d_inv(AminB, temp3)
call dgemm('N', 'N', ov_space, ov_space, ov_space, 1.0_dp, temp3, ov_space, AplusB, ov_space, 0.0_dp, temp2, ov_space)
deallocate(temp3)
do i = 1, ov_space
!Mulitply each column by corresponding eigenvalue^(1/2). Only deal with positive eigenvalues here
do j = 1, ov_space
temp2(j, i) = temp2(j, i) * (sqrt(sqrt(W(i))))
end do
end do
allocate(X_Chol(ov_space, ov_space))
allocate(Y_Chol(ov_space, ov_space))
do i = 1, ov_space
do j = 1, ov_space
X_Chol(j, i) = 0.5_dp * (temp(j, i) + temp2(j, i))
Y_Chol(j, i) = 0.5_dp * (temp(j, i) - temp2(j, i))
end do
end do
! do i=1,ov_space
! do j=1,ov_space
! if(abs(Y_Chol(i,j)).gt.1.0e-7) then
! write(stdout,*) "Found non-zero Y matrix value..."
! end if
! end do
! end do
!Y_Chol and X_Chol now cover all eigenvectors.
!Eigensystem can be represented as:
! ( X ) val=sqrt(W) and ( Y* ) val=-sqrt(W)
! ( Y ) ( X* )
call Check_XY_orthogonality(X_Chol, Y_Chol)
!Are X and Y vectors the same as from the stability matrix?
!I guess there will be different rotation among degenerate sets, so not a well defined comparison.
! if(tStabilityAnalysis) then
! do i=1,ov_space
! do j=1,ov_space
! if(abs(X_Chol(j,i)-X_Stab(j,i)).gt.1.0e-7) then
! write(stdout,*) j,i,X_Chol(j,i),X_Stab(j,i)
!! call stop_all(t_r,"Calculation of X matrix not the same as via stability matrix")
! end if
! if(abs(Y_Chol(j,i)-Y_Stab(j,i)).gt.1.0e-7) then
! write(stdout,*) j,i,Y_Chol(j,i),Y_Stab(j,i)
!! call stop_all(t_r,"Calculation of Y matrix not the same as via stability matrix")
! end if
! end do
! end do
! end if
!Now check that we satisfy the original RPA equations
!For the *positive* eigenvalue space (since we have extracted eigenvectors corresponding to this), check that:
! ( A B ) (X) = E_v(X )
! ( B A ) (Y) (-Y)
deallocate(temp, temp2)
allocate(Stability(StabilitySize, StabilitySize))
Stability(:, :) = 0.0_dp
Stability(1:ov_space, 1:ov_space) = A_mat(1:ov_space, 1:ov_space)
Stability(ov_space + 1:StabilitySize, 1:ov_space) = B_mat(1:ov_space, 1:ov_space)
Stability(1:ov_space, ov_space + 1:StabilitySize) = B_mat(1:ov_space, 1:ov_space)
Stability(ov_space + 1:StabilitySize, ov_space + 1:StabilitySize) = A_mat(1:ov_space, 1:ov_space)
allocate(temp(StabilitySize, ov_space))
temp(1:ov_space, 1:ov_space) = X_Chol(:, :)
temp(ov_space + 1:StabilitySize, 1:ov_space) = Y_Chol(:, :)
allocate(temp2(StabilitySize, ov_space))
call dgemm('n', 'n', StabilitySize, ov_space, StabilitySize, 1.0_dp, Stability, StabilitySize, temp, &
StabilitySize, 0.0_dp, temp2, StabilitySize)
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp2(j, i) - (sqrt(W(i)) * X_Chol(j, i))) > 1.0e-6_dp) then
write(stdout, *) i, j, temp2(j, i), (sqrt(W(i)) * X_Chol(j, i)), sqrt(W(i))
call stop_all(t_r, "RPA equations not satisfied for positive frequencies in X matrix")
end if
end do
end do
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp2(j + ov_space, i) - (-sqrt(W(i)) * Y_Chol(j, i))) > 1.0e-6_dp) then
write(stdout, *) i, j, temp2(j + ov_space, i), (-sqrt(W(i)) * Y_Chol(j, i)), -sqrt(W(i))
call stop_all(t_r, "RPA equations not satisfied for positive frequencies in Y matrix")
end if
end do
end do
deallocate(temp, temp2)
!Is is also satisfied the other way around?
!Check that we also satisfy (still for the *positive* eigenvalues):
! ( A B ) (Y) = -E_v(Y )
! ( B A ) (X) (-X)
allocate(temp(StabilitySize, ov_space))
temp(1:ov_space, 1:ov_space) = Y_Chol(:, :)
temp(ov_space + 1:StabilitySize, 1:ov_space) = X_Chol(:, :)
allocate(temp2(StabilitySize, ov_space))
call dgemm('n', 'n', StabilitySize, ov_space, StabilitySize, 1.0_dp, Stability, StabilitySize, temp, &
StabilitySize, 0.0_dp, temp2, StabilitySize)
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp2(j, i) - (-sqrt(W(i)) * Y_Chol(j, i))) > 1.0e-6_dp) then
call stop_all(t_r, "RPA equations not satisfied for negative frequencies in X matrix")
end if
end do
end do
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp2(j + ov_space, i) - (sqrt(W(i)) * X_Chol(j, i))) > 1.0e-6_dp) then
call stop_all(t_r, "RPA equations not satisfied for negative frequencies in Y matrix")
end if
end do
end do
deallocate(temp, temp2)
!TODO: Finally, check that we satisfy eq. 1 in the Scuseria paper for X and Y defined for positive eigenvalues...
do i = ov_space + 1, StabilitySize
do j = 1, StabilitySize
Stability(i, j) = -Stability(i, j)
end do
end do
!Stability copy is now (A B // -B -A)
allocate(temp(StabilitySize, ov_space))
allocate(temp2(StabilitySize, ov_space))
temp = 0.0_dp
temp(1:ov_space, 1:ov_space) = X_Chol(:, :)
temp(ov_space + 1:StabilitySize, 1:ov_space) = Y_Chol(:, :)
call dgemm('n', 'n', StabilitySize, ov_space, StabilitySize, 1.0_dp, Stability, StabilitySize, temp, &
StabilitySize, 0.0_dp, temp2, StabilitySize)
do i = 1, ov_space
do j = 1, ov_space
if (abs(temp2(j, i) - (sqrt(W(i)) * X_Chol(j, i))) > 1.0e-7_dp) then
call stop_all(t_r, "RPA equations not satisfied for X")
end if
end do
end do
do i = 1, ov_space
do j = ov_space + 1, StabilitySize
if (abs(temp2(j, i) - (sqrt(W(i)) * Y_Chol(j - ov_space, i))) > 1.0e-7_dp) then
call stop_all(t_r, "RPA equations not satisfied for Y")
end if
end do
end do
deallocate(temp, temp2, Stability)
Energy2 = real(HDiagTemp, dp)
norm = 0.0_dp
do i = 1, ov_space
norm = norm + A_mat(i, i)
end do
Energy2 = Energy2 - norm / 2.0_dp
norm = 0.0_dp
do i = 1, ov_space
!Run over positive eigenvalue pairs
norm = norm + sqrt(W(i))
end do
Energy2 = Energy2 + norm / 2.0_dp
if (tDirectRPA) then
write(stdout, "(A,G25.10)") "Direct RPA energy (plasmonic RPA-TDA excitation energies): ", &
Energy2 - real(HDiagTemp, dp)
else
write(stdout, "(A,G25.10)") "Full RPA energy (plasmonic RPA-TDA excitation energies): ", &
Energy2 - real(HDiagTemp, dp)
end if
Energy = 0.0_dp
allocate(temp(ov_space, ov_space))
allocate(temp2(ov_space, ov_space))
temp(:, :) = 0.0_dp
call d_inv(X_Chol, temp)
call dgemm('n', 'n', ov_space, ov_space, ov_space, 1.0_dp, Y_Chol, ov_space, temp, ov_space, 0.0_dp, temp2, ov_space)
call dgemm('n', 'n', ov_space, ov_space, ov_space, 1.0_dp, B_Mat, ov_space, temp2, ov_space, 0.0_dp, temp, ov_space)
do i = 1, ov_space
Energy = Energy + temp(i, i)
end do
Energy = Energy / 2.0_dp
if (tDirectRPA) then
write(stdout, "(A,G25.10)") "Direct RPA energy (Ring-CCD: 1/2 Tr[BZ]): ", Energy
else
write(stdout, "(A,G25.10)") "Full RPA energy (Ring-CCD: 1/2 Tr[BZ]): ", Energy
end if
deallocate(temp, temp2)
HDiagTemp = get_helement(fDet, fDet, 0)
Energy = real(HDiagTemp, dp)
do v = 1, ov_space
norm = 0.0_dp
do i = 1, ov_space
norm = norm + (abs(Y_Chol(i, v)))**2.0_dp
end do
Energy = Energy - norm * sqrt(W(v))
end do
if (tDirectRPA) then
write(stdout, "(A,G25.10)") "Direct RPA energy (Y matrix): ", Energy - real(HDiagTemp, dp)
else
write(stdout, "(A,G25.10)") "Full RPA energy (Y matrix): ", Energy - real(HDiagTemp, dp)
end if
write(stdout, "(A)")
write(stdout, "(A)") "RPA calculation completed successfully"
write(stdout, "(A)")
write(stdout, "(A)") "Calculating 1RDM..."
allocate(OccNumbers(nBasis)) !This is the GS 1RDM, which is diagonal in RPA
OccNumbers(:) = 0.0_dp
!For occupied orbital: Gamma(i,i) = 1 - 1/2 \sum_{mu,a} |Y^{mu}_ia|^2
!For virtual orbital: Gamma(a,a) = 1/2 \sum_{mu,i} |Y^{mu}_ia|^2
do i = 1, nel
do mu = 1, ov_space
do a = virt_start, nBasis
ia_ind = ov_space_ind(a, i)
OccNumbers(i) = OccNumbers(i) + 0.5_dp * abs(Y_Chol(ia_ind, mu))**2.0_dp
end do
end do
OccNumbers(i) = 1.0_dp - OccNumbers(i)
end do
do a = virt_start, nBasis
do mu = 1, ov_space
do i = 1, nel
ia_ind = ov_space_ind(a, i)
OccNumbers(a) = OccNumbers(a) + 0.5_dp * abs(Y_Chol(ia_ind, mu))**2.0_dp
end do
end do
end do
write(stdout, "(A)") "Writing RPA occupation numbers to file: RPA_OccNumbers"
write(stdout, "(A)") ""
iunit = get_free_unit()
open(iunit, file='RPA_OccNumbers', status='unknown')
write(iunit, "(A)") "# Orbital(Energy_order) Orbital Fock eigenvalue RPA_Occupation"
do i = 1, nBasis
i_ind = Brr(i)
write(iunit, "(2I9,2G25.10)") i, i_ind, Arr(i_ind, 2), OccNumbers(i)
end do
close(iunit)
deallocate(A_mat, B_mat, X_Chol, Y_Chol, OccNumbers)
end subroutine RunRPA_QBA