From 6ca04d4b44626701fd9d8af4af9fc35fd835eadf Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 6 Oct 2019 20:08:38 +0200 Subject: [PATCH] pp-RPA --- examples/molecule.H2 | 2 +- input/basis | 21 ++++++++----- input/options | 2 +- input/weight | 21 ++++++++----- src/QuAcK/linear_response_B_pp.f90 | 6 ++-- src/QuAcK/linear_response_C_pp.f90 | 13 +++++---- src/QuAcK/linear_response_D_pp.f90 | 11 ++++--- src/QuAcK/linear_response_pp.f90 | 47 +++++++++++++++++++++--------- 8 files changed, 81 insertions(+), 42 deletions(-) diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 32aa31b..81c624a 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -2,4 +2,4 @@ 2 1 1 0 0 # Znuc x y z H 0. 0. 0. - H 0. 0. 0.6 + H 0. 0. 1.4 diff --git a/input/basis b/input/basis index b9ca7b5..4c61da7 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,16 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 +1 6 +S 4 1.00 + 234.0000000 0.0025870 + 35.1600000 0.0195330 + 7.9890000 0.0909980 + 2.2120000 0.2720500 S 1 1.00 - 0.2976000 1.0000000 + 0.6669000 1.0000000 +S 1 1.00 + 0.2089000 1.0000000 P 1 1.00 - 1.2750000 1.0000000 + 3.0440000 1.0000000 +P 1 1.00 + 0.7580000 1.0000000 +D 1 1.00 + 1.9650000 1.0000000 diff --git a/input/options b/input/options index 1b1eb4a..d1a5689 100644 --- a/input/options +++ b/input/options @@ -5,7 +5,7 @@ # CC: maxSCF thresh DIIS n_diis 64 0.00001 F 1 # CIS/TDHF/BSE: singlet triplet - T F + T T # GF: maxSCF thresh DIIS n_diis renormalization 64 0.00001 T 10 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta diff --git a/input/weight b/input/weight index b9ca7b5..4c61da7 100644 --- a/input/weight +++ b/input/weight @@ -1,9 +1,16 @@ -1 3 -S 3 1.00 - 38.3600000 0.0238090 - 5.7700000 0.1548910 - 1.2400000 0.4699870 +1 6 +S 4 1.00 + 234.0000000 0.0025870 + 35.1600000 0.0195330 + 7.9890000 0.0909980 + 2.2120000 0.2720500 S 1 1.00 - 0.2976000 1.0000000 + 0.6669000 1.0000000 +S 1 1.00 + 0.2089000 1.0000000 P 1 1.00 - 1.2750000 1.0000000 + 3.0440000 1.0000000 +P 1 1.00 + 0.7580000 1.0000000 +D 1 1.00 + 1.9650000 1.0000000 diff --git a/src/QuAcK/linear_response_B_pp.f90 b/src/QuAcK/linear_response_B_pp.f90 index 8ab931a..f3f5dd6 100644 --- a/src/QuAcK/linear_response_B_pp.f90 +++ b/src/QuAcK/linear_response_B_pp.f90 @@ -39,14 +39,14 @@ subroutine linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) ab = 0 do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do b=a+1,nBas-nR ab = ab + 1 ij = 0 do i=nC+1,nO - do j=nC+1,nO + do j=i+1,nO ij = ij + 1 - B_pp(ab,ij) = (1d0 + delta_spin)*ERI(a,b,i,j) - (1d0 - delta_dRPA)*ERI(a,b,j,j) + B_pp(ab,ij) = (1d0 + delta_spin)*ERI(a,b,i,j) - (1d0 - delta_dRPA)*ERI(a,b,j,i) enddo enddo diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/QuAcK/linear_response_C_pp.f90 index ef3b3ee..6f5599c 100644 --- a/src/QuAcK/linear_response_C_pp.f90 +++ b/src/QuAcK/linear_response_C_pp.f90 @@ -16,6 +16,7 @@ subroutine linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) double precision :: delta_spin double precision :: delta_dRPA + double precision :: eF double precision,external :: Kronecker_delta integer :: a,b,c,d,ab,cd @@ -35,19 +36,21 @@ subroutine linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) delta_dRPA = 0d0 if(dRPA) delta_dRPA = 1d0 -! Build A matrix +! Build C matrix + eF = e(nO) + e(nO+1) + ab = 0 do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do b=a+1,nBas-nR ab = ab + 1 cd = 0 do c=nO+1,nBas-nR - do d=nO+1,nBas-nR + do d=c+1,nBas-nR cd = cd + 1 - C_pp(ab,cd) = (e(a) + e(b))*Kronecker_delta(a,c)*Kronecker_delta(b,d) & - + (1d0 + delta_spin)*ERI(a,b,c,d) - (1d0 - delta_dRPA)*ERI(a,b,d,c) + C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + (1d0 + delta_spin)*ERI(a,b,c,d) - (1d0 - delta_dRPA)*ERI(a,b,d,c) enddo enddo diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/QuAcK/linear_response_D_pp.f90 index 28455b6..4daeab5 100644 --- a/src/QuAcK/linear_response_D_pp.f90 +++ b/src/QuAcK/linear_response_D_pp.f90 @@ -14,6 +14,7 @@ subroutine linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) ! Local variables + double precision :: eF double precision :: delta_spin double precision :: delta_dRPA double precision,external :: Kronecker_delta @@ -35,18 +36,20 @@ subroutine linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) delta_dRPA = 0d0 if(dRPA) delta_dRPA = 1d0 -! Build A matrix +! Build D matrix + + eF = e(nO) + e(nO+1) ij = 0 do i=nC+1,nO - do j=nC+1,nO + do j=i+1,nO ij = ij + 1 kl = 0 do k=nC+1,nO - do l=nC+1,nO + do l=k+1,nO kl = kl + 1 - D_pp(ij,kl) = - (e(i) + e(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + (1d0 + delta_spin)*ERI(k,l,i,j) - (1d0 - delta_dRPA)*ERI(k,l,j,i) enddo diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index 3b9874a..6b7339e 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -1,6 +1,6 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,Ec_ppRPA) -! Compute the p-p channel of the linear response +! Compute the p-p channel of the linear response: see Scueria et al. JCP 139, 104113 (2013) implicit none include 'parameters.h' @@ -24,7 +24,14 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,E double precision,allocatable :: C(:,:) double precision,allocatable :: D(:,:) double precision,allocatable :: M(:,:) + double precision,allocatable :: Z(:,:) double precision,allocatable :: w(:) + double precision,allocatable :: w1(:) + double precision,allocatable :: w2(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) ! Output variables @@ -32,12 +39,13 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,E ! Useful quantities - nOO = nO*nO - nVV = nV*nV + nOO = nO*(nO-1)/2 + nVV = nV*(nV-1)/2 ! Memory allocation - allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),w(nOO+nVV)) + allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),w(nOO+nVV), & + w1(nVV),w2(nOO),X1(nVV,nVV),Y1(nVV,nOO),X2(nOO,nVV),Y2(nOO,nOO)) ! Build B, C and D matrices for the pp channel @@ -56,8 +64,8 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,E ! Diagonal blocks - M( 1:nVV , 1:nVV) = +C(1:nVV,1:nVV) - M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = -D(1:nOO,1:nOO) + M( 1:nVV , 1:nVV) = + C(1:nVV,1:nVV) + M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - D(1:nOO,1:nOO) ! Off-diagonal blocks @@ -69,22 +77,33 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,E ! Diagonalize the p-h matrix - call diagonalize_matrix(nOO+nVV,M(:,:),w(:)) + Z(:,:) = M(:,:) + call diagonalize_matrix(nOO+nVV,Z(:,:),w(:)) -! Build X+Y + write(*,*) 'pp-RPA excitation energies' + call matout(nOO+nVV,1,w(:)) + write(*,*) -! XpY(1:nS,1:nS) = M(nS+1:2*nS,1:nS) + M(nS+1:2*nS,nS+1:2*nS) +! Split the various quantities in p-p and h-h parts -! print*,'X+Y' -! call matout(nS,nS,XpY) + w1(:) = w(nOO+1:nOO+nVV) + w2(:) = w(1:nOO) - print*,'pp-RPA excitation energies' - call matout(nOO+nVV,1,w) + X1(:,:) = Z(nOO+1:nOO+nVV, 1:nVV ) + Y1(:,:) = Z(nOO+1:nOO+nVV,nVV+1:nOO+nVV) + X2(:,:) = Z( 1:nOO , 1:nVV ) + Y2(:,:) = Z( 1:nOO ,nVV+1:nOO+nVV) + + if(minval(w1(:)) < 0d0) call print_warning('You may have instabilities in pp-RPA!!') + if(maxval(w2(:)) > 0d0) call print_warning('You may have instabilities in pp-RPA!!') ! Compute the RPA correlation energy - Ec_ppRPA = 0.5d0*(sum(abs(w(:))) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:))) + Ec_ppRPA = 0.5d0*( sum(w1(:)) - sum(w2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) print*,'Ec(pp-RPA) = ',Ec_ppRPA + print*,'Ec(pp-RPA) = ',0.5d0*( sum(abs(w(:))) - trace_matrix(nVV*nOO,M(:,:))) + print*,'Ec(pp-RPA) = ',+sum(w1(:)) - trace_matrix(nVV,C(:,:)) + print*,'Ec(pp-RPA) = ',-sum(w2(:)) - trace_matrix(nOO,D(:,:)) end subroutine linear_response_pp