diff --git a/input/methods b/input/methods index 5b87e02..02ac24b 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) @@ -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 F # G0F2* evGF2* qsGF2* G0F3 evGF3 - F F F F F + 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..f8f9987 100644 --- a/input/options +++ b/input/options @@ -1,18 +1,18 @@ # 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 T 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 + 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 5 T 0.0 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 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 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.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..8d163a5 --- /dev/null +++ b/src/LR/unrestricted_linear_response_B_pp.f90 @@ -0,0 +1,119 @@ +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=a,nBas-nR(1) + ab = ab + 1 + ij = 0 + do i=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 + end do + + ! bbbb block + + ab = 0 + do a=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=i+1,nO(2) + ij = ij + 1 + + B_pp(nPaa+ab,nHaa+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(ab,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..0de2351 --- /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=a,nBas-nR(1) + ab = ab + 1 + cd = 0 + do c=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 + end do + + ! bbbb block + + ab = 0 + do a=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=c,nBas-nR(2) + cd = cd + 1 + + 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 +!----------------------------------------------- + + 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(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 + 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..4ea7fd8 --- /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=i+1,nO(1) + ij = ij + 1 + kl = 0 + do k=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) & + + 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=i+1,nO(2) + ij = ij + 1 + kl = 0 + do k=nC(2)+1,nO(2) + do l=k+1,nO(2) + kl = kl + 1 + + 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 + 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(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 + 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..6c9019a --- /dev/null +++ b/src/LR/unrestricted_linear_response_pp.f90 @@ -0,0 +1,113 @@ +subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt, & +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 + + 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 + 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)) +!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) + 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)) + +!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) + + ! 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/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 new file mode 100644 index 0000000..ec7bb5f --- /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*abs(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*abs(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*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_diag diff --git a/src/MP/MP2.f90 b/src/MP/MP2.f90 index 7dd921f..30e857e 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(-sigm1*Dijab))/Dijab + fs2 = (1d0 - exp(-sigm2*Dijab*Dijab))/Dijab + fk = (1d0 - exp(-kappa*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 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 diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90 new file mode 100644 index 0000000..59513a4 --- /dev/null +++ b/src/RPA/ppURPA.f90 @@ -0,0 +1,146 @@ +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 :: 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_ppURPA(nspin) + double precision :: EcAC(nspin) + +! Hello world + + write(*,*) + write(*,*)'****************************************' + write(*,*)'| particle-particle URPA calculation |' + write(*,*)'****************************************' + write(*,*) + +! Initialization + + Ec_ppURPA(:) = 0d0 + EcAC(:) = 0d0 + +! Useful quantities + +!spin-conserved quantities + + nPaa = nV(1)*(nV(1)-1)/2 + nPbb = nV(2)*(nV(2)-1)/2 + + nP_sc = nPaa + nPbb + + nHaa = nO(1)*(nO(1)-1)/2 + 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(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(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 + + if(spin_conserved) then + + ispin = 1 + +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 + +! Spin-flip manifold + + if(spin_flip) then + + ispin = 2 + +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,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_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(*,*) + +! 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