diff --git a/input/methods b/input/methods index a1580bd..c1c696e 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 T F F # CCD pCCD DCD CCSD CCSD(T) @@ -9,9 +9,9 @@ # 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 F F F # G0T0 evGT qsGT 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(*,*)