From 1dfee1e9fa8719e37ff805394fb52063a6f38037 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 13 Dec 2021 11:40:16 +0100 Subject: [PATCH 1/8] set Fermi energy to zero in pp --- input/methods | 4 ++-- input/options | 6 +++--- mol/h2.xyz | 2 +- src/LR/linear_response_C_pp.f90 | 4 ++-- src/LR/linear_response_D_pp.f90 | 4 ++-- 5 files changed, 10 insertions(+), 10 deletions(-) diff --git a/input/methods b/input/methods index 5b87e02..77c21ce 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - F F T F + T F F F # MP2* MP3 MP2-F12 F F F # CCD pCCD DCD CCSD CCSD(T) @@ -13,7 +13,7 @@ # 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 F F F # MCMP2 diff --git a/input/options b/input/options index 2592c87..1e505e3 100644 --- a/input/options +++ b/input/options @@ -1,15 +1,15 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability - 1024 0.0000001 T 5 1 1 F F + 1024 0.00001 T 5 1 1 F F # MP: # CC: maxSCF thresh DIIS n_diis - 64 0.0000001 T 5 + 64 0.00001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.00367493 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 - 256 0.00001 T 5 T 0.00 F F F F F + 256 0.00001 T 15 T 0.00367493 F F F F F # ACFDT: AC Kx XBS F T F # BSE: BSE dBSE dTDA evDyn diff --git a/mol/h2.xyz b/mol/h2.xyz index 83cc4a4..21fde66 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 3.2 +H 0. 0. 0.5 diff --git a/src/LR/linear_response_C_pp.f90 b/src/LR/linear_response_C_pp.f90 index 724dc59..b9315a7 100644 --- a/src/LR/linear_response_C_pp.f90 +++ b/src/LR/linear_response_C_pp.f90 @@ -25,8 +25,8 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C_pp ! Define the chemical potential - eF = e(nO) + e(nO+1) -! eF = 0d0 +! eF = e(nO) + e(nO+1) + eF = 0d0 ! Build C matrix for the singlet manifold diff --git a/src/LR/linear_response_D_pp.f90 b/src/LR/linear_response_D_pp.f90 index 398b96a..e05ed6c 100644 --- a/src/LR/linear_response_D_pp.f90 +++ b/src/LR/linear_response_D_pp.f90 @@ -25,8 +25,8 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D_pp ! Define the chemical potential - eF = e(nO) + e(nO+1) -! eF = 0d0 +! eF = e(nO) + e(nO+1) + eF = 0d0 ! Build the D matrix for the singlet manifold From 3d3806da641e6a9ae9f092c5b471692f6dd5ef98 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 13 Dec 2021 14:28:05 +0100 Subject: [PATCH 2/8] ppURPA --- input/methods | 4 +- src/RPA/ppURPA.f90 | 133 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 135 insertions(+), 2 deletions(-) create mode 100644 src/RPA/ppURPA.f90 diff --git a/input/methods b/input/methods index 77c21ce..0444112 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 T F F + F F F F F # G0T0 evGT qsGT - F F F + T F F # MCMP2 F # * unrestricted version available diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90 new file mode 100644 index 0000000..ac09ec1 --- /dev/null +++ b/src/RPA/ppURPA.f90 @@ -0,0 +1,133 @@ +subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUHF,ERI_aaaa,ERI_aabb,ERI_bbbb,e) + +! Perform unrestricted pp-RPA calculations + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA + logical,intent(in) :: doACFDT + logical,intent(in) :: spin_conserved + logical,intent(in) :: spin_flip + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + double precision,intent(in) :: ENuc + double precision,intent(in) :: EUHF + 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 + + 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(:,:) + + double precision :: Ec_ppRPA(nspin) + double precision :: EcAC(nspin) + +! Hello world + + write(*,*) + write(*,*)'****************************************' + write(*,*)'| particle-particle URPA calculation |' + write(*,*)'****************************************' + write(*,*) + +! Initialization + + Ec_ppRPA(:) = 0d0 + EcAC(:) = 0d0 + +! 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 + + ! Memory allocation + +! allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & +! Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs)) + +! allocate(Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & +! Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt)) + +! Spin-conserved manifold + + if(spin_conserved) then + + ispin = 1 + +! call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,e,ERI, & +! Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,Ec_ppRPA(ispin)) + +! call print_excitation('pp-RPA (N+2)',5,nVVs,Omega1s) +! call print_excitation('pp-RPA (N-2)',5,nOOs,Omega2s) + + endif + +! Spin-flip manifold + + if(spin_flip) then + + ispin = 2 + +! call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,e,ERI, & +! Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,Ec_ppRPA(ispin)) + +! call print_excitation('pp-RPA (N+2)',6,nVVt,Omega1t) +! call print_excitation('pp-RPA (N-2)',6,nOOt,Omega2t) + + endif + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-conserved) =',Ec_ppRPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-flip) =',3d0*Ec_ppRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + EUHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + +! Compute the correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '---------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of pp-RPA correlation energy' +! write(*,*) '---------------------------------------------------------' +! write(*,*) + +! call ACFDT_pp(TDA,singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI,e,EcAC) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (singlet) =',EcAC(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy (triplet) =',EcAC(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA correlation energy =',EcAC(1) + EcAC(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@ppRPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + + +end subroutine ppURPA From fcf46af6002af165d50e0f8a75a224c5bda1c2d8 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 13 Dec 2021 14:30:49 +0100 Subject: [PATCH 3/8] modif QuAcK for ppURPA --- src/QuAcK/QuAcK.f90 | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index edaa4ab..c6b0444 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -824,7 +824,15 @@ program QuAcK if(doppRPA) then call cpu_time(start_RPA) - call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) + if(unrestricted) then + + call ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUHF,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,eHF) + + else + + call ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) + + end if call cpu_time(end_RPA) t_RPA = end_RPA - start_RPA From 5d5824ce14efe0cabd25d0f2aa449f0c480fd66f Mon Sep 17 00:00:00 2001 From: EnzoMonino Date: Wed, 15 Dec 2021 11:41:06 +0100 Subject: [PATCH 4/8] 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)) From 55187b39c46ffd620cf633ee13d8bf61fa744b5b Mon Sep 17 00:00:00 2001 From: EnzoMonino Date: Thu, 16 Dec 2021 12:28:28 +0100 Subject: [PATCH 5/8] ppURPA updates --- input/methods | 8 +- input/options | 4 +- mol/h2.xyz | 2 +- src/LR/linear_response_pp.f90 | 10 ++- src/LR/unrestricted_linear_response_B_pp.f90 | 19 +++-- src/LR/unrestricted_linear_response_C_pp.f90 | 22 +++--- src/LR/unrestricted_linear_response_D_pp.f90 | 12 +-- src/LR/unrestricted_linear_response_pp.f90 | 39 ++++++---- src/RPA/ppURPA.f90 | 81 +++++++++++--------- 9 files changed, 112 insertions(+), 85 deletions(-) diff --git a/input/methods b/input/methods index 77c21ce..4836a6a 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - T F F F + F T F F # MP2* MP3 MP2-F12 F F F # CCD pCCD DCD CCSD CCSD(T) @@ -9,11 +9,11 @@ # CIS* CIS(D) CID CISD FCI F F F F F # RPA* RPAx* crRPA ppRPA - F F F F + F F F T # G0F2* evGF2* qsGF2* G0F3 evGF3 - F F F F F + F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - F F T F F + F F F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 37b1199..d6745cf 100644 --- a/input/options +++ b/input/options @@ -1,11 +1,11 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability - 1024 0.00001 T 5 1 1 F F + 1024 0.00001 T 5 1 1 T F # MP: # CC: maxSCF thresh DIIS n_diis 64 0.00001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T T T T + F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.00367493 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 diff --git a/mol/h2.xyz b/mol/h2.xyz index e329359..21fde66 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 1.2 +H 0. 0. 0.5 diff --git a/src/LR/linear_response_pp.f90 b/src/LR/linear_response_pp.f90 index a1decea..894d111 100644 --- a/src/LR/linear_response_pp.f90 +++ b/src/LR/linear_response_pp.f90 @@ -47,7 +47,8 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om ! Memory allocation allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Omega(nOO+nVV)) - +write(*,*) 'nOO', nOO +write(*,*) 'nVV', nVV !-------------------------------------------------! ! Solve the p-p eigenproblem ! !-------------------------------------------------! @@ -62,7 +63,7 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om call linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,C) call linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,D) - +!call matout(nVV,nVV,C) if(TDA) then X1(:,:) = +C(:,:) @@ -76,7 +77,8 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om else call linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,B) - +!call matout(nVV,nOO,B) +!call matout(nOO,nOO,D) ! Diagonal blocks M( 1:nVV , 1:nVV) = + C(1:nVV,1:nVV) @@ -86,7 +88,7 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om M( 1:nVV ,nVV+1:nOO+nVV) = - B(1:nVV,1:nOO) M(nVV+1:nOO+nVV, 1:nVV) = + transpose(B(1:nVV,1:nOO)) - +call matout(nOO+nVV,nOO+nVV,M) ! Diagonalize the p-h matrix if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z) diff --git a/src/LR/unrestricted_linear_response_B_pp.f90 b/src/LR/unrestricted_linear_response_B_pp.f90 index 5f055ad..8d163a5 100644 --- a/src/LR/unrestricted_linear_response_B_pp.f90 +++ b/src/LR/unrestricted_linear_response_B_pp.f90 @@ -46,21 +46,21 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP ! Build B matrix for spin-conserving transitions !----------------------------------------------- - if(ispin == 1) then + if(ispin == 1) then ! aaaa block ab = 0 do a=nO(1)+1,nBas-nR(1) - do b=nO(1)+1,nBas-nR(1) + do b=a,nBas-nR(1) ab = ab + 1 ij = 0 do i=nC(1)+1,nO(1) - do j=nC(1)+1,nO(1) + do j=i+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 @@ -70,14 +70,14 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP ab = 0 do a=nO(2)+1,nBas-nR(2) - do b=nO(2)+1,nBas-nR(2) + do b=a+1,nBas-nR(2) ab = ab + 1 ij = 0 do i=nC(2)+1,nO(2) - do j=nC(2)+1,nO(2) + do j=i+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)) + B_pp(nPaa+ab,nHaa+ij) = lambda*(ERI_bbbb(a,b,i,j) - ERI_bbbb(a,b,j,i)) end do end do @@ -104,7 +104,10 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP 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) + + B_pp(ab,ij) = lambda*ERI_aabb(a,b,i,j) + + end do end do end do diff --git a/src/LR/unrestricted_linear_response_C_pp.f90 b/src/LR/unrestricted_linear_response_C_pp.f90 index fbcfa6f..0de2351 100644 --- a/src/LR/unrestricted_linear_response_C_pp.f90 +++ b/src/LR/unrestricted_linear_response_C_pp.f90 @@ -48,16 +48,16 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP ab = 0 do a=nO(1)+1,nBas-nR(1) - do b=nO(1)+1,nBas-nR(1) + do b=a,nBas-nR(1) ab = ab + 1 cd = 0 do c=nO(1)+1,nBas-nR(1) - do d=nO(1)+1,nBas-nR(1) + do d=c,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)) - +!write(*,*) C_pp(ab,cd) end do end do end do @@ -67,23 +67,23 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP ab = 0 do a=nO(2)+1,nBas-nR(2) - do b=nO(2)+1,nBas-nR(2) + do b=a,nBas-nR(2) ab = ab + 1 cd = 0 do c=nO(2)+1,nBas-nR(2) - do d=nO(2)+1,nBas-nR(2) + do d=c,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 + C_pp(nPaa+ab,nPaa+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)) +!write(*,*) 'nPaa+ab',nPaa+ab + end do end do end do end do end if - +! !----------------------------------------------- ! Build C matrix for spin-flip transitions !----------------------------------------------- @@ -102,7 +102,7 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP 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) & + C_pp(ab,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 diff --git a/src/LR/unrestricted_linear_response_D_pp.f90 b/src/LR/unrestricted_linear_response_D_pp.f90 index 0845eb5..4ea7fd8 100644 --- a/src/LR/unrestricted_linear_response_D_pp.f90 +++ b/src/LR/unrestricted_linear_response_D_pp.f90 @@ -48,11 +48,11 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH ij = 0 do i=nC(1)+1,nO(1) - do j=nC(1)+1,nO(1) + do j=i+1,nO(1) ij = ij + 1 kl = 0 do k=nC(1)+1,nO(1) - do l=nC(1)+1,nO(1) + do l=k+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) & @@ -67,14 +67,14 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH ij = 0 do i=nC(2)+1,nO(2) - do j=nC(2)+1,nO(2) + do j=i+1,nO(2) ij = ij + 1 kl = 0 do k=nC(2)+1,nO(2) - do l=nC(2)+1,nO(2) + do l=k+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) & + D_pp(nHaa+ij,nHaa+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 @@ -102,7 +102,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH 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)& + D_pp(ij,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 diff --git a/src/LR/unrestricted_linear_response_pp.f90 b/src/LR/unrestricted_linear_response_pp.f90 index a4ac53b..6c9019a 100644 --- a/src/LR/unrestricted_linear_response_pp.f90 +++ b/src/LR/unrestricted_linear_response_pp.f90 @@ -1,5 +1,5 @@ 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,& +nHaa,nHab,nHbb,nHt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1,X1,Y1,Omega2,X2,Y2,& EcRPA) ! Compute linear response for unrestricted formalism @@ -23,8 +23,7 @@ EcRPA) integer,intent(in) :: nHaa integer,intent(in) :: nHab integer,intent(in) :: nHbb - integer,intent(in) :: nHt - integer,intent(in) :: nS_sc + 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) @@ -56,20 +55,28 @@ EcRPA) ! Memory allocation - allocate(C(nPt,nPt),B(nPt,nHt),D(nHt,nHt),M(nPt+nHt,nPt+nHt),Z(nPt+nHt,nPt+nHt),& + allocate(C(nPt,nPt),B(nPt,nHt),D(nHt,nHt),M(nPt+nHt,nPt+nHt),Z(nPt+nHt,nPt+nHt),& Omega(nPt+nHt)) - +!write(*,*) 'ispin', ispin +!write(*,*) 'nPt', nPt +!write(*,*) 'nHt', 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 matout(nPt,nPt,C) +!write(*,*) 'Hello' 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 matout(nPt,nHt,B) +!write(*,*) 'Hello' call unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,lambda,& ERI_aaaa,ERI_aabb,ERI_bbbb,D) +!call matout(nHt,nHt,D) +!write(*,*) 'Hello' ! Diagonal blocks M( 1:nPt , 1:nPt) = + C(1:nPt,1:nPt) @@ -80,23 +87,27 @@ ERI_aaaa,ERI_aabb,ERI_bbbb,D) 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)) +!call matout(nPt+nHt,nPt+nHt,M) + ! Diagonalize the p-h matrix - if(nHt+nPt > 0) call diagonalize_general_matrix(nHt+nPt,M,Omega,Z) +! 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(:,:)) +! call sort_ppRPA(nHt,nPt,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),& +!Y2(:,:)) - ! end if Pourquoi ne faut-il pas de end if ici ? +! 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 !!!' +! 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 54fe6e8..59513a4 100644 --- a/src/RPA/ppURPA.f90 +++ b/src/RPA/ppURPA.f90 @@ -25,18 +25,17 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ! Local variables - integer :: ispin - integer :: nS - 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(:,:) + integer :: ispin + integer :: nPaa,nPbb,nPab,nP_sc,nP_sf + integer :: nHaa,nHbb,nHab,nH_sc,nH_sf + double precision,allocatable :: Omega1sc(:),Omega1sf(:) + double precision,allocatable :: X1sc(:,:),X1sf(:,:) + double precision,allocatable :: Y1sc(:,:),Y1sf(:,:) + double precision,allocatable :: Omega2sc(:),Omega2sf(:) + double precision,allocatable :: X2sc(:,:),X2sf(:,:) + double precision,allocatable :: Y2sc(:,:),Y2sf(:,:) - double precision :: Ec_ppRPA(nspin) + double precision :: Ec_ppURPA(nspin) double precision :: EcAC(nspin) ! Hello world @@ -49,28 +48,38 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ! Initialization - Ec_ppRPA(:) = 0d0 + Ec_ppURPA(:) = 0d0 EcAC(:) = 0d0 ! Useful quantities +!spin-conserved quantities + nPaa = nV(1)*(nV(1)-1)/2 - nPab = nV(1)*nV(2) - nPbb = nV(2)*nV(2) - nPt = nPaa + nPab + nPbb + nPbb = nV(2)*(nV(2)-1)/2 + + nP_sc = nPaa + nPbb nHaa = nO(1)*(nO(1)-1)/2 - nHab = nO(1)*nO(2) - nHbb = nO(2)*nO(2) - nHt = nHaa + nHab + nHbb + nHbb = nO(2)*(nO(2)-1)/2 + nH_sc = nHaa + nHbb + +!spin-flip quantities + + nPab = nV(1)*nV(2) + nHab = nO(1)*nO(2) + + nP_sf = nPab + nH_sf = nPab + ! Memory allocation - allocate(Omega1(nPt),X1(nPt,nPt),Y1(nHt,nPt), & - Omega2(nHt),X2(nPt,nHt),Y2(nHt,nHt)) + allocate(Omega1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), & + Omega2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc)) -! allocate(Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & -! Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt)) + allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & + Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) ! Spin-conserved manifold @@ -78,11 +87,12 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ispin = 1 -! call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,e,ERI, & -! Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,Ec_ppRPA(ispin)) - -! call print_excitation('pp-RPA (N+2)',5,nVVs,Omega1s) -! call print_excitation('pp-RPA (N-2)',5,nOOs,Omega2s) +call unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nP_sc, & +nHaa,nHab,nHbb,nH_sc,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1sc,X1sc,Y1sc,Omega2sc,X2sc,Y2sc,& +Ec_ppURPA(ispin)) +write(*,*) 'Hello!' + call print_excitation('pp-RPA (N+2)',5,nP_sc,Omega1sc) + call print_excitation('pp-RPA (N-2)',5,nH_sc,Omega2sc) endif @@ -92,20 +102,21 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ispin = 2 -! call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,e,ERI, & -! Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,Ec_ppRPA(ispin)) +call unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nP_sf, & +nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf,Omega2sf,X2sf,Y2sf,& +Ec_ppURPA(ispin)) -! call print_excitation('pp-RPA (N+2)',6,nVVt,Omega1t) -! call print_excitation('pp-RPA (N-2)',6,nOOt,Omega2t) + call print_excitation('pp-RPA (N+2)',6,nP_sf,Omega1sf) + call print_excitation('pp-RPA (N-2)',6,nH_sf,Omega2sf) endif write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-conserved) =',Ec_ppRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-flip) =',3d0*Ec_ppRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + EUHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-conserved) =',Ec_ppURPA(1) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy (spin-flip) =',3d0*Ec_ppURPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA correlation energy =',Ec_ppURPA(1) + 3d0*Ec_ppURPA(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppRPA total energy =',ENuc + EUHF + Ec_ppURPA(1) + 3d0*Ec_ppURPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) From e006fb5a01e4fb3850997af242e1268df06f6bb9 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 16 Dec 2021 14:24:02 +0100 Subject: [PATCH 6/8] regularized MP2 --- input/methods | 4 +- src/MP/MP2.f90 | 130 ++++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 114 insertions(+), 20 deletions(-) diff --git a/input/methods b/input/methods index 77c21ce..a1580bd 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF KS MOM T F F F # MP2* MP3 MP2-F12 - F F F + T F F # CCD pCCD DCD CCSD CCSD(T) F F F F F # drCCD rCCD crCCD lCCD @@ -13,7 +13,7 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - F F T F F + F F F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/src/MP/MP2.f90 b/src/MP/MP2.f90 index 7dd921f..cfca705 100644 --- a/src/MP/MP2.f90 +++ b/src/MP/MP2.f90 @@ -1,6 +1,6 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) -! Perform second-order Moller-Plesset calculation +! Perform second-order Moller-Plesset calculation with and without regularizers implicit none @@ -19,7 +19,15 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) ! Local variables integer :: i,j,a,b - double precision :: eps,E2a,E2b + + double precision :: kappa,sigm1,sigm2 + double precision :: Dijab + double precision :: fs,fs2,fk + + double precision :: E2d,E2ds,E2ds2,E2dk + double precision :: E2x,E2xs,E2xs2,E2xk + + double precision :: EcsMP2,Ecs2MP2,EckMP2 ! Output variables @@ -33,43 +41,129 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) write(*,*)'************************************************' write(*,*) -! Compute MP2 energy +!---------------------------------------------! +! Parameters for regularized MP2 calculations ! +!---------------------------------------------! + + kappa = 1.1d0 + sigm1 = 0.7d0 + sigm2 = 0.4d0 + +!--------------------------------------------------! +! Compute conventinal and regularized MP2 energies ! +!--------------------------------------------------! + + E2d = 0d0 + E2ds = 0d0 + E2ds2 = 0d0 + E2dk = 0d0 + + E2x = 0d0 + E2xs = 0d0 + E2xs2 = 0d0 + E2xk = 0d0 - E2a = 0d0 - E2b = 0d0 do i=nC+1,nO - do j=nC+1,nO - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do j=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR - eps = e(i) + e(j) - e(a) - e(b) + Dijab = e(a) + e(b) - e(i) - e(j) -! Second-order ring diagram + ! Second-order ring diagram - E2a = E2a + ERI(i,j,a,b)*ERI(i,j,a,b)/eps + fs = (1d0 - exp(-kappa*Dijab))/Dijab + fs2 = (1d0 - exp(-sigm1*Dijab*Dijab))/Dijab + fk = (1d0 - exp(-sigm2*Dijab))**2/Dijab -! Second-order exchange diagram + E2d = E2d - ERI(i,j,a,b)*ERI(i,j,a,b)/Dijab + E2ds = E2ds - ERI(i,j,a,b)*ERI(i,j,a,b)*fs + E2ds2 = E2ds2 - ERI(i,j,a,b)*ERI(i,j,a,b)*fs2 + E2dk = E2dk - ERI(i,j,a,b)*ERI(i,j,a,b)*fk - E2b = E2b + ERI(i,j,a,b)*ERI(i,j,b,a)/eps + ! Second-order exchange diagram - enddo + E2x = E2x - ERI(i,j,a,b)*ERI(i,j,b,a)/Dijab + E2xs = E2xs - ERI(i,j,a,b)*ERI(i,j,b,a)*fs + E2xs2 = E2xs2 - ERI(i,j,a,b)*ERI(i,j,b,a)*fs2 + E2xk = E2xk - ERI(i,j,a,b)*ERI(i,j,b,a)*fk + + enddo + enddo enddo enddo - enddo - EcMP2 = 2d0*E2a - E2b + EcMP2 = 2d0*E2d - E2x + EcsMP2 = 2d0*E2ds - E2xs + Ecs2MP2 = 2d0*E2ds2 - E2xs2 + EckMP2 = 2d0*E2dk - E2xk + +!------------! +! MP2 energy ! +!------------! write(*,*) write(*,'(A32)') '--------------------------' write(*,'(A32)') ' MP2 calculation ' write(*,'(A32)') '--------------------------' write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcMP2 - write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2a - write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2b + write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2d + write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2x write(*,'(A32)') '--------------------------' write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EcMP2 write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EcMP2 write(*,'(A32)') '--------------------------' write(*,*) +!-------------------! +! sigma1-MP2 energy ! +!-------------------! + + write(*,*) + write(*,'(A32)') '--------------------------' + write(*,'(A32)') ' sigma-MP2 calculation ' + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcsMP2 + write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2ds + write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EcsMP2 + write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EcsMP2 + write(*,'(A32)') '--------------------------' + write(*,*) + +!--------------------! +! sigma^2-MP2 energy ! +!--------------------! + + write(*,*) + write(*,'(A32)') '--------------------------' + write(*,'(A32)') ' sigma^2-MP2 calculation ' + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',Ecs2MP2 + write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2ds2 + write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs2 + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + Ecs2MP2 + write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + Ecs2MP2 + write(*,'(A32)') '--------------------------' + write(*,*) + +!------------------! +! kappa-MP2 energy ! +!------------------! + + write(*,*) + write(*,'(A32)') '--------------------------' + write(*,'(A32)') ' kappa-MP2 calculation ' + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EckMP2 + write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2dk + write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xk + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EckMP2 + write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EckMP2 + write(*,'(A32)') '--------------------------' + write(*,*) + end subroutine MP2 From a64bdce3c463ff5a2ce6163c32755e787e0d6288 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 16 Dec 2021 14:47:13 +0100 Subject: [PATCH 7/8] regularized GW --- ...gularized_self_energy_correlation_diag.f90 | 124 ++++++++++++++++++ src/MP/MP2.f90 | 6 +- 2 files changed, 127 insertions(+), 3 deletions(-) create mode 100644 src/MBPT/regularized_self_energy_correlation_diag.f90 diff --git a/src/MBPT/regularized_self_energy_correlation_diag.f90 b/src/MBPT/regularized_self_energy_correlation_diag.f90 new file mode 100644 index 0000000..d0dba7d --- /dev/null +++ b/src/MBPT/regularized_self_energy_correlation_diag.f90 @@ -0,0 +1,124 @@ +subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) + +! Compute diagonal of the correlation part of the regularized self-energy + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: COHSEX + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,p,q,jb + double precision :: eps + double precision,external :: SigC_dcgw + + double precision :: kappa + double precision :: fk + +! Output variables + + double precision,intent(out) :: SigC(nBas) + double precision,intent(out) :: EcGM + +! Initialize + + SigC(:) = 0d0 + +!---------------------------------------------! +! Parameters for regularized MP2 calculations ! +!---------------------------------------------! + + kappa = 1.1d0 + +!----------------------------- +! COHSEX static self-energy +!----------------------------- + + if(COHSEX) then + + ! COHSEX: SEX part of the COHSEX correlation self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do jb=1,nS + SigC(p) = SigC(p) + 4d0*rho(p,i,jb)**2/Omega(jb) + end do + end do + end do + + ! COHSEX: COH part of the COHSEX correlation self-energy + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do jb=1,nS + SigC(p) = SigC(p) - 2d0*rho(p,q,jb)**2/Omega(jb) + end do + end do + end do + + ! GM correlation energy + + EcGM = 0d0 + do i=nC+1,nO + EcGM = EcGM - SigC(i) + end do + +!----------------------------- +! GW self-energy +!----------------------------- + + else + + ! Occupied part of the correlation self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do jb=1,nS + eps = e(p) - e(i) + Omega(jb) + fk = (1d0 - exp(-kappa*eps))**2/eps + SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*fk + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do jb=1,nS + eps = e(p) - e(a) - Omega(jb) + fk = (1d0 - exp(-kappa*eps))**2/eps + SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*fk + end do + end do + end do + + ! GM correlation energy + + EcGM = 0d0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do jb=1,nS + eps = e(a) - e(i) + Omega(jb) + fk = (1d0 - exp(-kappa*eps))**2/eps + EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk + end do + end do + end do + + end if + +end subroutine regularized_self_energy_correlation_diag diff --git a/src/MP/MP2.f90 b/src/MP/MP2.f90 index cfca705..30e857e 100644 --- a/src/MP/MP2.f90 +++ b/src/MP/MP2.f90 @@ -72,9 +72,9 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) ! Second-order ring diagram - fs = (1d0 - exp(-kappa*Dijab))/Dijab - fs2 = (1d0 - exp(-sigm1*Dijab*Dijab))/Dijab - fk = (1d0 - exp(-sigm2*Dijab))**2/Dijab + fs = (1d0 - exp(-sigm1*Dijab))/Dijab + fs2 = (1d0 - exp(-sigm2*Dijab*Dijab))/Dijab + fk = (1d0 - exp(-kappa*Dijab))**2/Dijab E2d = E2d - ERI(i,j,a,b)*ERI(i,j,a,b)/Dijab E2ds = E2ds - ERI(i,j,a,b)*ERI(i,j,a,b)*fs From d24b729fc32d70e15263761b326e0f019dc29dab Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Thu, 16 Dec 2021 15:19:44 +0100 Subject: [PATCH 8/8] regularized GW --- input/methods | 8 +- input/options | 2 +- src/MBPT/evGW.f90 | 2 + src/MBPT/qsGW.f90 | 2 + .../regularized_self_energy_correlation.f90 | 126 ++++++++++++++++++ ...gularized_self_energy_correlation_diag.f90 | 6 +- 6 files changed, 138 insertions(+), 8 deletions(-) create mode 100644 src/MBPT/regularized_self_energy_correlation.f90 diff --git a/input/methods b/input/methods index c1c696e..02ac24b 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF KS MOM - F T F F + T F F F # MP2* MP3 MP2-F12 - T F F + F F F # CCD pCCD DCD CCSD CCSD(T) F F F F F # drCCD rCCD crCCD lCCD @@ -9,11 +9,11 @@ # CIS* CIS(D) CID CISD FCI F F F F F # RPA* RPAx* crRPA ppRPA - F F F T + F F F F # 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 F F F # MCMP2 diff --git a/input/options b/input/options index d6745cf..f8f9987 100644 --- a/input/options +++ b/input/options @@ -9,7 +9,7 @@ # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.00367493 3 # 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 + 256 0.00001 T 5 T 0.0 F F F F F # ACFDT: AC Kx XBS F F F # BSE: BSE dBSE dTDA evDyn diff --git a/src/MBPT/evGW.f90 b/src/MBPT/evGW.f90 index 9a74734..c25f9ab 100644 --- a/src/MBPT/evGW.f90 +++ b/src/MBPT/evGW.f90 @@ -161,11 +161,13 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(G0W) then +! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) else +! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) diff --git a/src/MBPT/qsGW.f90 b/src/MBPT/qsGW.f90 index d6bc2bb..2b7b289 100644 --- a/src/MBPT/qsGW.f90 +++ b/src/MBPT/qsGW.f90 @@ -192,11 +192,13 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(G0W) then +! call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) else +! call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) diff --git a/src/MBPT/regularized_self_energy_correlation.f90 b/src/MBPT/regularized_self_energy_correlation.f90 new file mode 100644 index 0000000..1a56dcd --- /dev/null +++ b/src/MBPT/regularized_self_energy_correlation.f90 @@ -0,0 +1,126 @@ +subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) + +! Compute correlation part of the regularized self-energy + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: COHSEX + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,j,a,b + integer :: p,q,r + integer :: jb + double precision :: eps + + double precision :: kappa + double precision :: fk + +! Output variables + + double precision,intent(out) :: SigC(nBas,nBas) + double precision,intent(out) :: EcGM + +! Initialize + + SigC(:,:) = 0d0 + +!---------------------------------------------! +! Parameters for regularized MP2 calculations ! +!---------------------------------------------! + + kappa = 1.1d0 + +!-----------------------------! +! COHSEX static approximation ! +!-----------------------------! + + if(COHSEX) then + + ! COHSEX: SEX of the COHSEX correlation self-energy + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do i=nC+1,nO + do jb=1,nS + SigC(p,q) = SigC(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb) + end do + end do + end do + end do + + ! COHSEX: COH part of the COHSEX correlation self-energy + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do r=nC+1,nBas-nR + do jb=1,nS + SigC(p,q) = SigC(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb) + end do + end do + end do + end do + + EcGM = 0d0 + do i=nC+1,nO + EcGM = EcGM + 0.5d0*SigC(i,i) + end do + + else + +!----------------! +! GW self-energy ! +!----------------! + + ! Occupied part of the correlation self-energy + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do i=nC+1,nO + do jb=1,nS + eps = e(p) - e(i) + Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*fk + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do a=nO+1,nBas-nR + do jb=1,nS + eps = e(p) - e(a) - Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*fk + end do + end do + end do + end do + + ! GM correlation energy + + EcGM = 0d0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do jb=1,nS + eps = e(a) - e(i) + Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk + end do + end do + end do + + end if + +end subroutine regularized_self_energy_correlation diff --git a/src/MBPT/regularized_self_energy_correlation_diag.f90 b/src/MBPT/regularized_self_energy_correlation_diag.f90 index d0dba7d..ec7bb5f 100644 --- a/src/MBPT/regularized_self_energy_correlation_diag.f90 +++ b/src/MBPT/regularized_self_energy_correlation_diag.f90 @@ -88,7 +88,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, do i=nC+1,nO do jb=1,nS eps = e(p) - e(i) + Omega(jb) - fk = (1d0 - exp(-kappa*eps))**2/eps + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*fk end do end do @@ -100,7 +100,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, do a=nO+1,nBas-nR do jb=1,nS eps = e(p) - e(a) - Omega(jb) - fk = (1d0 - exp(-kappa*eps))**2/eps + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*fk end do end do @@ -113,7 +113,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, do a=nO+1,nBas-nR do jb=1,nS eps = e(a) - e(i) + Omega(jb) - fk = (1d0 - exp(-kappa*eps))**2/eps + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk end do end do