From b3f376a81b94b9e35f09dc4c1bfc07d5c36f1043 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 7 Oct 2019 22:31:45 +0200 Subject: [PATCH] debug ppRPA --- input/options | 2 +- src/QuAcK/linear_response_B_pp.f90 | 70 ++++++++++++++++------------ src/QuAcK/linear_response_C_pp.f90 | 75 ++++++++++++++++++------------ src/QuAcK/linear_response_D_pp.f90 | 75 ++++++++++++++++++------------ src/QuAcK/linear_response_pp.f90 | 11 ++--- src/QuAcK/ppRPA.f90 | 59 ++++++++++++----------- 6 files changed, 169 insertions(+), 123 deletions(-) diff --git a/input/options b/input/options index da3080c..721bbaa 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/src/QuAcK/linear_response_B_pp.f90 b/src/QuAcK/linear_response_B_pp.f90 index f3f5dd6..0dc6cbd 100644 --- a/src/QuAcK/linear_response_B_pp.f90 +++ b/src/QuAcK/linear_response_B_pp.f90 @@ -1,4 +1,4 @@ -subroutine linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) +subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) ! Compute the B matrix of the pp channel @@ -8,14 +8,11 @@ subroutine linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) ! Input variables integer,intent(in) :: ispin - logical,intent(in) :: dRPA integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) ! Local variables - double precision :: delta_spin - double precision :: delta_dRPA double precision,external :: Kronecker_delta integer :: a,b,i,j,ab,ij @@ -24,33 +21,48 @@ subroutine linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) double precision,intent(out) :: B_pp(nVV,nOO) -! Singlet or triplet manifold? - - delta_spin = 0d0 - if(ispin == 1) delta_spin = +1d0 - if(ispin == 2) delta_spin = -1d0 - -! Direct RPA - - delta_dRPA = 0d0 - if(dRPA) delta_dRPA = 1d0 - -! Build A matrix +! Build B matrix for the singlet manifold - ab = 0 - do a=nO+1,nBas-nR - do b=a+1,nBas-nR - ab = ab + 1 - ij = 0 - do i=nC+1,nO - do j=i+1,nO - ij = ij + 1 + if(ispin == 1) then - B_pp(ab,ij) = (1d0 + delta_spin)*ERI(a,b,i,j) - (1d0 - delta_dRPA)*ERI(a,b,j,i) + ab = 0 + do a=nO+1,nBas-nR + do b=nO+1,a + ab = ab + 1 + ij = 0 + do i=nC+1,nO + do j=nC+1,i + ij = ij + 1 + + B_pp(ab,ij) = (ERI(a,b,i,j) + ERI(a,b,j,i))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) + + end do + end do + end do + end do - enddo - enddo - enddo - enddo + end if + +! Build the B matrix for the triplet manifold + + if(ispin == 2) then + + ab = 0 + do a=nO+1,nBas-nR + do b=nO+1,a-1 + ab = ab + 1 + ij = 0 + do i=nC+1,nO + do j=nC+1,i-1 + ij = ij + 1 + + B_pp(ab,ij) = ERI(a,b,i,j) - ERI(a,b,j,i) + + end do + end do + end do + end do + + end if end subroutine linear_response_B_pp diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/QuAcK/linear_response_C_pp.f90 index 6f5599c..c68b257 100644 --- a/src/QuAcK/linear_response_C_pp.f90 +++ b/src/QuAcK/linear_response_C_pp.f90 @@ -1,4 +1,4 @@ -subroutine linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) +subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) ! Compute the C matrix of the pp channel @@ -8,14 +8,11 @@ subroutine linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) ! Input variables integer,intent(in) :: ispin - logical,intent(in) :: dRPA integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) ! Local variables - double precision :: delta_spin - double precision :: delta_dRPA double precision :: eF double precision,external :: Kronecker_delta @@ -25,36 +22,54 @@ subroutine linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) double precision,intent(out) :: C_pp(nVV,nVV) -! Singlet or triplet manifold? - - delta_spin = 0d0 - if(ispin == 1) delta_spin = +1d0 - if(ispin == 2) delta_spin = -1d0 - -! Direct RPA - - delta_dRPA = 0d0 - if(dRPA) delta_dRPA = 1d0 - -! Build C matrix +! Define the chemical potential eF = e(nO) + e(nO+1) - ab = 0 - do a=nO+1,nBas-nR - do b=a+1,nBas-nR - ab = ab + 1 - cd = 0 - do c=nO+1,nBas-nR - do d=c+1,nBas-nR - cd = cd + 1 +! Build C matrix for the singlet manifold - 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) + if(ispin == 1) then - enddo - enddo - enddo - enddo + ab = 0 + do a=nO+1,nBas-nR + do b=nO+1,a + ab = ab + 1 + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,c + cd = cd + 1 + + C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + (ERI(a,b,c,d) + ERI(a,b,d,c))/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + + end do + end do + end do + end do + + end if + +! Build C matrix for the triplet manifold + + if(ispin == 2) then + + ab = 0 + do a=nO+1,nBas-nR + do b=nO+1,a-1 + ab = ab + 1 + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,c-1 + cd = cd + 1 + + C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + ERI(a,b,c,d) - ERI(a,b,d,c) + + end do + end do + end do + end do + + end if end subroutine linear_response_C_pp diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/QuAcK/linear_response_D_pp.f90 index 4daeab5..6220297 100644 --- a/src/QuAcK/linear_response_D_pp.f90 +++ b/src/QuAcK/linear_response_D_pp.f90 @@ -1,4 +1,4 @@ -subroutine linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) +subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) ! Compute the D matrix of the pp channel @@ -8,15 +8,12 @@ subroutine linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) ! Input variables integer,intent(in) :: ispin - logical,intent(in) :: dRPA integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) ! Local variables double precision :: eF - double precision :: delta_spin - double precision :: delta_dRPA double precision,external :: Kronecker_delta integer :: i,j,k,l,ij,kl @@ -25,36 +22,54 @@ subroutine linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) double precision,intent(out) :: D_pp(nOO,nOO) -! Singlet or triplet manifold? - - delta_spin = 0d0 - if(ispin == 1) delta_spin = +1d0 - if(ispin == 2) delta_spin = -1d0 - -! Direct RPA - - delta_dRPA = 0d0 - if(dRPA) delta_dRPA = 1d0 - -! Build D matrix +! Define the chemical potential eF = e(nO) + e(nO+1) - ij = 0 - do i=nC+1,nO - do j=i+1,nO - ij = ij + 1 - kl = 0 - do k=nC+1,nO - do l=k+1,nO - kl = kl + 1 +! Build the D matrix for the singlet manifold - 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) + if(ispin == 1) then - enddo - enddo - enddo - enddo + ij = 0 + do i=nC+1,nO + do j=nC+1,i + ij = ij + 1 + kl = 0 + do k=nC+1,nO + do l=nC+1,k + kl = kl + 1 + + D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + + ERI(k,l,i,j) - ERI(k,l,j,i) + + end do + end do + end do + end do + + end if + +! Build the D matrix for the triplet manifold + + if(ispin == 2) then + + ij = 0 + do i=nC+1,nO + do j=nC+1,i-1 + ij = ij + 1 + kl = 0 + do k=nC+1,nO + do l=nC+1,k-1 + kl = kl + 1 + + D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + + (ERI(k,l,i,j) + ERI(k,l,j,i))/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + + end do + end do + end do + end do + + end if end subroutine linear_response_D_pp diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index 044c8dc..3a546fa 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -1,5 +1,4 @@ -subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & - Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA) +subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA) ! Compute the p-p channel of the linear response: see Scueria et al. JCP 139, 104113 (2013) @@ -8,8 +7,6 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, ! Input variables - logical,intent(in) :: dRPA - logical,intent(in) :: TDA logical,intent(in) :: BSE integer,intent(in) :: ispin,nBas,nC,nO,nV,nR integer,intent(in) :: nOO @@ -43,9 +40,9 @@ subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, ! Build B, C and D matrices for the pp channel - call linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B) - call linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C) - call linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D) + call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B) + call linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C) + call linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D) !------------------------------------------------------------------------ ! Solve the p-p eigenproblem diff --git a/src/QuAcK/ppRPA.f90 b/src/QuAcK/ppRPA.f90 index c5447df..ae29426 100644 --- a/src/QuAcK/ppRPA.f90 +++ b/src/QuAcK/ppRPA.f90 @@ -21,8 +21,6 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER ! Local variables - logical :: dRPA - logical :: TDA logical :: BSE integer :: ispin integer :: nOO @@ -44,45 +42,40 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER write(*,*)'****************************************' write(*,*) -! Useful quantities - - nOO = nO*(nO-1)/2 - nVV = nV*(nV-1)/2 - ! Initialization Ec_ppRPA(:) = 0d0 -! Switch on exchange for TDHF - - dRPA = .false. - -! Switch off Tamm-Dancoff approximation for TDHF - - TDA = .false. - ! Switch off Bethe-Salpeter equation for TDHF BSE = .false. -! Memory allocation - - allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & - Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin)) - ! Singlet manifold if(singlet_manifold) then ispin = 1 - call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & + ! Useful quantities + + nOO = nO*(nO+1)/2 + nVV = nV*(nV+1)/2 + + ! Memory allocation + + allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & + Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin)) + + call linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & + Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & + Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & Ec_ppRPA(ispin)) + call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin)) + deallocate(Omega1,X1,Y1,Omega2,X2,Y2) + endif ! Triplet manifold @@ -91,13 +84,27 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER ispin = 2 - call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & - Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & - Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & + ! Useful quantities + + nOO = nO*(nO-1)/2 + nVV = nV*(nV-1)/2 + + ! Memory allocation + + allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), & + Omega2(nOO,nspin),X2(nVV,nOO,nspin),Y2(nOO,nOO,nspin)) + + + call linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI, & + Omega1(:,ispin),X1(:,:,ispin),Y1(:,:,ispin), & + Omega2(:,ispin),X2(:,:,ispin),Y2(:,:,ispin), & Ec_ppRPA(ispin)) + call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1(:,ispin)) call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2(:,ispin)) + deallocate(Omega1,X1,Y1,Omega2,X2,Y2) + endif write(*,*)