From 5d5824ce14efe0cabd25d0f2aa449f0c480fd66f Mon Sep 17 00:00:00 2001 From: EnzoMonino Date: Wed, 15 Dec 2021 11:41:06 +0100 Subject: [PATCH] ppURPA --- input/methods | 4 +- input/options | 4 +- mol/h2.xyz | 2 +- src/LR/unrestricted_linear_response.f90 | 2 +- src/LR/unrestricted_linear_response_B_pp.f90 | 116 +++++++++++++++++++ src/LR/unrestricted_linear_response_C_pp.f90 | 115 ++++++++++++++++++ src/LR/unrestricted_linear_response_D_pp.f90 | 115 ++++++++++++++++++ src/LR/unrestricted_linear_response_pp.f90 | 102 ++++++++++++++++ src/RPA/ppURPA.f90 | 36 +++--- 9 files changed, 473 insertions(+), 23 deletions(-) create mode 100644 src/LR/unrestricted_linear_response_B_pp.f90 create mode 100644 src/LR/unrestricted_linear_response_C_pp.f90 create mode 100644 src/LR/unrestricted_linear_response_D_pp.f90 create mode 100644 src/LR/unrestricted_linear_response_pp.f90 diff --git a/input/methods b/input/methods index 0444112..77c21ce 100644 --- a/input/methods +++ b/input/methods @@ -13,9 +13,9 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - F F F F F + F F T F F # G0T0 evGT qsGT - T F F + F F F # MCMP2 F # * unrestricted version available diff --git a/input/options b/input/options index 1e505e3..37b1199 100644 --- a/input/options +++ b/input/options @@ -11,8 +11,8 @@ # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 256 0.00001 T 15 T 0.00367493 F F F F F # ACFDT: AC Kx XBS - F T F + F F F # BSE: BSE dBSE dTDA evDyn - T T T F + F F F F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/mol/h2.xyz b/mol/h2.xyz index 21fde66..e329359 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 0.5 +H 0. 0. 1.2 diff --git a/src/LR/unrestricted_linear_response.f90 b/src/LR/unrestricted_linear_response.f90 index a23e073..b940f45 100644 --- a/src/LR/unrestricted_linear_response.f90 +++ b/src/LR/unrestricted_linear_response.f90 @@ -106,7 +106,7 @@ subroutine unrestricted_linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR, if(minval(Omega) < 0d0) & call print_warning('You may have instabilities in linear response: negative excitations!!') - + ! do ia=1,nSt ! if(Omega(ia) < 0d0) Omega(ia) = 0d0 ! end do diff --git a/src/LR/unrestricted_linear_response_B_pp.f90 b/src/LR/unrestricted_linear_response_B_pp.f90 new file mode 100644 index 0000000..5f055ad --- /dev/null +++ b/src/LR/unrestricted_linear_response_B_pp.f90 @@ -0,0 +1,116 @@ +subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa, & + nHab,nHbb,nHt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,B_pp) + +! Compute linear response + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nPaa + integer,intent(in) :: nPab + integer,intent(in) :: nPbb + integer,intent(in) :: nPt + integer,intent(in) :: nHaa + integer,intent(in) :: nHab + integer,intent(in) :: nHbb + integer,intent(in) :: nHt + double precision,intent(in) :: lambda + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: eF + double precision,external :: Kronecker_delta + + integer :: i,j,a,b,c,d,ij,ab + +! Output variables + + double precision,intent(out) :: B_pp(nPt,nHt) + + + eF = 0d0 + +!----------------------------------------------- +! Build B matrix for spin-conserving transitions +!----------------------------------------------- + + if(ispin == 1) then + + ! aaaa block + + ab = 0 + do a=nO(1)+1,nBas-nR(1) + do b=nO(1)+1,nBas-nR(1) + ab = ab + 1 + ij = 0 + do i=nC(1)+1,nO(1) + do j=nC(1)+1,nO(1) + ij = ij + 1 + + B_pp(ab,ij) = lambda*(ERI_aaaa(a,b,i,j) - ERI_aaaa(a,b,j,i)) + + end do + end do + end do + end do + + ! bbbb block + + ab = 0 + do a=nO(2)+1,nBas-nR(2) + do b=nO(2)+1,nBas-nR(2) + ab = ab + 1 + ij = 0 + do i=nC(2)+1,nO(2) + do j=nC(2)+1,nO(2) + ij = ij + 1 + + B_pp(nPaa+nPab+ab,nHaa+nHab+ij) = lambda*(ERI_bbbb(a,b,i,j) - ERI_bbbb(a,b,j,i)) + + end do + end do + end do + end do + + end if + +!----------------------------------------------- +! Build B matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + + B_pp(:,:) = 0d0 + + ! abab block + + ab = 0 + do a=nO(1)+1,nBas-nR(1) + do b=nO(2)+1,nBas-nR(2) + ab = ab + 1 + ij = 0 + do i=nC(1)+1,nO(1) + do j=nC(2)+1,nO(2) + ij = ij + 1 + B_pp(nPaa+ab,nHaa+ij) = lambda*ERI_aabb(a,b,i,j) + end do + end do + end do + end do + + end if + + +end subroutine unrestricted_linear_response_B_pp diff --git a/src/LR/unrestricted_linear_response_C_pp.f90 b/src/LR/unrestricted_linear_response_C_pp.f90 new file mode 100644 index 0000000..fbcfa6f --- /dev/null +++ b/src/LR/unrestricted_linear_response_C_pp.f90 @@ -0,0 +1,115 @@ +subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,lambda,& + e,ERI_aaaa,ERI_aabb,ERI_bbbb,C_pp) + +! Compute linear response + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nPaa + integer,intent(in) :: nPab + integer,intent(in) :: nPbb + integer,intent(in) :: nPt + double precision,intent(in) :: lambda + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: eF + double precision,external :: Kronecker_delta + + integer :: i,j,a,b,c,d,ab,cd + +! Output variables + + double precision,intent(out) :: C_pp(nPt,nPt) + + + eF = 0d0 + +!----------------------------------------------- +! Build C matrix for spin-conserving transitions +!----------------------------------------------- + + if(ispin == 1) then + + ! aaaa block + + ab = 0 + do a=nO(1)+1,nBas-nR(1) + do b=nO(1)+1,nBas-nR(1) + ab = ab + 1 + cd = 0 + do c=nO(1)+1,nBas-nR(1) + do d=nO(1)+1,nBas-nR(1) + cd = cd + 1 + + C_pp(ab,cd) = (e(a,1) + e(b,1) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + lambda*(ERI_aaaa(a,b,c,d) - ERI_aaaa(a,b,d,c)) + + end do + end do + end do + end do + + ! bbbb block + + ab = 0 + do a=nO(2)+1,nBas-nR(2) + do b=nO(2)+1,nBas-nR(2) + ab = ab + 1 + cd = 0 + do c=nO(2)+1,nBas-nR(2) + do d=nO(2)+1,nBas-nR(2) + cd = cd + 1 + + C_pp(nPaa+nPab+ab,nPaa+nPab+cd) = (e(a,2) + e(b,2) - eF)*Kronecker_delta(a,c) & + *Kronecker_delta(b,d) + lambda*(ERI_bbbb(a,b,c,d) - ERI_bbbb(a,b,d,c)) + + end do + end do + end do + end do + + end if + +!----------------------------------------------- +! Build C matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + + C_pp(:,:) = 0d0 + + ! abab block + + ab = 0 + do a=nO(1)+1,nBas-nR(1) + do b=nO(2)+1,nBas-nR(2) + ab = ab + 1 + cd = 0 + do c=nO(1)+1,nBas-nR(1) + do d=nO(2)+1,nBas-nR(2) + cd = cd + 1 + C_pp(nPaa+ab,nPaa+cd) = (e(a,1) + e(b,2))*Kronecker_delta(a,c) & +*Kronecker_delta(b,c) + lambda*ERI_aabb(a,b,c,d) + end do + end do + end do + end do + + end if + + +end subroutine unrestricted_linear_response_C_pp diff --git a/src/LR/unrestricted_linear_response_D_pp.f90 b/src/LR/unrestricted_linear_response_D_pp.f90 new file mode 100644 index 0000000..0845eb5 --- /dev/null +++ b/src/LR/unrestricted_linear_response_D_pp.f90 @@ -0,0 +1,115 @@ +subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,lambda,& + e,ERI_aaaa,ERI_aabb,ERI_bbbb,D_pp) + +! Compute linear response + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nHaa + integer,intent(in) :: nHab + integer,intent(in) :: nHbb + integer,intent(in) :: nHt + double precision,intent(in) :: lambda + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: eF + double precision,external :: Kronecker_delta + + integer :: i,j,k,l,ij,kl + +! Output variables + + double precision,intent(out) :: D_pp(nHt,nHt) + + + eF = 0d0 + +!----------------------------------------------- +! Build D matrix for spin-conserving transitions +!----------------------------------------------- + + if(ispin == 1) then + + ! aaaa block + + ij = 0 + do i=nC(1)+1,nO(1) + do j=nC(1)+1,nO(1) + ij = ij + 1 + kl = 0 + do k=nC(1)+1,nO(1) + do l=nC(1)+1,nO(1) + kl = kl + 1 + + D_pp(ij,kl) = -(e(i,1) + e(j,1) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + + lambda*(ERI_aaaa(i,j,k,l) - ERI_aaaa(i,j,l,k)) + + end do + end do + end do + end do + + ! bbbb block + + ij = 0 + do i=nC(2)+1,nO(2) + do j=nC(2)+1,nO(2) + ij = ij + 1 + kl = 0 + do k=nC(2)+1,nO(2) + do l=nC(2)+1,nO(2) + kl = kl + 1 + + D_pp(nHaa+nHab+ij,nHaa+nHab+kl) = -(e(i,2) + e(j,2) - eF)*Kronecker_delta(i,k) & + *Kronecker_delta(j,l) + lambda*(ERI_bbbb(i,j,k,l) - ERI_bbbb(i,j,l,k)) + + end do + end do + end do + end do + + end if + +!----------------------------------------------- +! Build D matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + + D_pp(:,:) = 0d0 + + ! abab block + + ij = 0 + do i=nC(1)+1,nO(1) + do j=nC(2)+1,nO(2) + ij = ij + 1 + kl = 0 + do k=nC(1)+1,nO(1) + do l=nC(2)+1,nO(2) + kl = kl + 1 + D_pp(nHaa+ij,nHaa+kl) = -(e(i,1) + e(j,2))*Kronecker_delta(i,k)& + *Kronecker_delta(j,l) +lambda*ERI_aabb(i,j,k,l) + end do + end do + end do + end do + + end if + + +end subroutine unrestricted_linear_response_D_pp diff --git a/src/LR/unrestricted_linear_response_pp.f90 b/src/LR/unrestricted_linear_response_pp.f90 new file mode 100644 index 0000000..a4ac53b --- /dev/null +++ b/src/LR/unrestricted_linear_response_pp.f90 @@ -0,0 +1,102 @@ +subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt, & +nHaa,nHab,nHbb,nHt,nS_sc,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1,X1,Y1,Omega2,X2,Y2,& +EcRPA) + +! Compute linear response for unrestricted formalism + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + logical,intent(in) :: TDA + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nPaa + integer,intent(in) :: nPab + integer,intent(in) :: nPbb + integer,intent(in) :: nPt + integer,intent(in) :: nHaa + integer,intent(in) :: nHab + integer,intent(in) :: nHbb + integer,intent(in) :: nHt + integer,intent(in) :: nS_sc + double precision,intent(in) :: lambda + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + +! Local variables + + double precision,external :: trace_matrix + double precision :: EcRPA1 + double precision :: EcRPA2 + double precision,allocatable :: C(:,:) + double precision,allocatable :: B(:,:) + double precision,allocatable :: D(:,:) + double precision,allocatable :: M(:,:) + double precision,allocatable :: Z(:,:) + double precision,allocatable :: Omega(:) + +! Output variables + + double precision,intent(out) :: Omega1(nPt) + double precision,intent(out) :: X1(nPt,nPt) + double precision,intent(out) :: Y1(nHt,nPt) + double precision,intent(out) :: Omega2(nHt) + double precision,intent(out) :: X2(nPt,nHt) + double precision,intent(out) :: Y2(nHt,nHt) + double precision,intent(out) :: EcRPA + + +! Memory allocation + + allocate(C(nPt,nPt),B(nPt,nHt),D(nHt,nHt),M(nPt+nHt,nPt+nHt),Z(nPt+nHt,nPt+nHt),& +Omega(nPt+nHt)) + +! Build C, B and D matrices for the pp channel + + call unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,lambda,& +e,ERI_aaaa,ERI_aabb,ERI_bbbb,C) + + call unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,& +nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,B) + +call unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,lambda,& +ERI_aaaa,ERI_aabb,ERI_bbbb,D) + +! Diagonal blocks + + M( 1:nPt , 1:nPt) = + C(1:nPt,1:nPt) + M(nPt+1:nPt+nHt,nPt+1:nPt+nHt) = - D(1:nHt,1:nHt) + + ! Off-diagonal blocks + + M( 1:nPt ,nPt+1:nHt+nPt) = - B(1:nPt,1:nHt) + M(nPt+1:nHt+nPt, 1:nPt) = + transpose(B(1:nPt,1:nHt)) + +! Diagonalize the p-h matrix + + if(nHt+nPt > 0) call diagonalize_general_matrix(nHt+nPt,M,Omega,Z) + + ! Split the various quantities in p-p and h-h parts + + call sort_ppRPA(nHt,nPt,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),& +Y2(:,:)) + + ! end if Pourquoi ne faut-il pas de end if ici ? + +! Compute the RPA correlation energy + + EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nPt,C(:,:)) - trace_matrix(nHt,D(:,:)) ) + EcRPA1 = +sum(Omega1(:)) - trace_matrix(nPt,C(:,:)) + EcRPA2 = -sum(Omega2(:)) - trace_matrix(nHt,D(:,:)) + if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & + print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' + +end subroutine unrestricted_linear_response_pp diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90 index ac09ec1..54fe6e8 100644 --- a/src/RPA/ppURPA.f90 +++ b/src/RPA/ppURPA.f90 @@ -27,14 +27,14 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH integer :: ispin integer :: nS - integer :: nOOs,nOOt - integer :: nVVs,nVVt - double precision,allocatable :: Omega1s(:),Omega1t(:) - double precision,allocatable :: X1s(:,:),X1t(:,:) - double precision,allocatable :: Y1s(:,:),Y1t(:,:) - double precision,allocatable :: Omega2s(:),Omega2t(:) - double precision,allocatable :: X2s(:,:),X2t(:,:) - double precision,allocatable :: Y2s(:,:),Y2t(:,:) + integer :: nPaa,nPbb,nPab,nPt + integer :: nHaa,nHbb,nHab,nHt + double precision,allocatable :: Omega1(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + double precision,allocatable :: Omega2(:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) double precision :: Ec_ppRPA(nspin) double precision :: EcAC(nspin) @@ -54,18 +54,20 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ! Useful quantities -! nS = nO*nV - -! nOOs = nO*(nO+1)/2 -! nVVs = nV*(nV+1)/2 - -! nOOt = nO*(nO-1)/2 -! nVVt = nV*(nV-1)/2 + nPaa = nV(1)*(nV(1)-1)/2 + nPab = nV(1)*nV(2) + nPbb = nV(2)*nV(2) + nPt = nPaa + nPab + nPbb + + nHaa = nO(1)*(nO(1)-1)/2 + nHab = nO(1)*nO(2) + nHbb = nO(2)*nO(2) + nHt = nHaa + nHab + nHbb ! Memory allocation -! allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & -! Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs)) + allocate(Omega1(nPt),X1(nPt,nPt),Y1(nHt,nPt), & + Omega2(nHt),X2(nPt,nHt),Y2(nHt,nHt)) ! allocate(Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & ! Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt))