diff --git a/input/methods b/input/methods index f51432f..e2bef5b 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,13 +9,13 @@ # 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 - T F F F F + F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - T F F F F + F F F 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 5b13b5b..f9f075c 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # 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 @@ -9,12 +9,12 @@ # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 3 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg - 256 0.00001 T 5 T 0.0 F F F F F T + 256 0.00001 T 5 T 0.0 F F F F F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg 256 0.00001 T 5 T 0.0 F F # ACFDT: AC Kx XBS 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 3c8e04d..bb00204 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 1.5 +H 0. 0. 1.0 diff --git a/src/LR/unrestricted_linear_response_B_pp.f90 b/src/LR/unrestricted_linear_response_B_pp.f90 index 8d163a5..6f5db5b 100644 --- a/src/LR/unrestricted_linear_response_B_pp.f90 +++ b/src/LR/unrestricted_linear_response_B_pp.f90 @@ -1,5 +1,5 @@ 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) + nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,B_pp) ! Compute linear response @@ -22,8 +22,7 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP 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) :: lambda 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) @@ -43,56 +42,10 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP eF = 0d0 !----------------------------------------------- -! Build B matrix for spin-conserving transitions +! Build B matrix for spin-conserved 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 + if(ispin == 1) then ! abab block @@ -107,6 +60,52 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP B_pp(ab,ij) = lambda*ERI_aabb(a,b,i,j) + end do + end do + end do + end do + + end if + +!----------------------------------------------- +! Build B matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + + ! aaaa block + + ab = 0 + do a=nO(1)+1,nBas-nR(1) + do b=a+1,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 + end if + + if (ispin == 3) then + + ! 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(ab,ij) = lambda*(ERI_bbbb(a,b,i,j) - ERI_bbbb(a,b,j,i)) end do end do @@ -115,5 +114,4 @@ subroutine unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP 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 index 0de2351..45e8643 100644 --- a/src/LR/unrestricted_linear_response_C_pp.f90 +++ b/src/LR/unrestricted_linear_response_C_pp.f90 @@ -37,61 +37,12 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP eF = 0d0 - !----------------------------------------------- -! Build C matrix for spin-conserving transitions +! Build C matrix for spin-conserved 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 @@ -102,8 +53,8 @@ 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(ab,cd) = (e(a,1) + e(b,2))*Kronecker_delta(a,c) & -*Kronecker_delta(b,c) + lambda*ERI_aabb(a,b,c,d) + C_pp(ab,cd) = (e(a,1) + e(b,2))*Kronecker_delta(a,c) & +*Kronecker_delta(b,d) + lambda*ERI_aabb(a,b,c,d) end do end do end do @@ -111,5 +62,54 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP end if +!----------------------------------------------- +! Build C matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + + ! aaaa block + + ab = 0 + do a=nO(1)+1,nBas-nR(1) + do b=a+1,nBas-nR(1) + ab = ab + 1 + cd = 0 + do c=nO(1)+1,nBas-nR(1) + do d=c+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)) +!write(*,*) C_pp(ab,cd) + end do + end do + end do + end do + end if + + + if (ispin == 3) then + + ! bbbb block + + ab = 0 + do a=nO(2)+1,nBas-nR(2) + do b=a+1,nBas-nR(2) + ab = ab + 1 + cd = 0 + do c=nO(2)+1,nBas-nR(2) + do d=c+1,nBas-nR(2) + cd = cd + 1 + + C_pp(ab,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 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 index 4ea7fd8..02ab1d0 100644 --- a/src/LR/unrestricted_linear_response_D_pp.f90 +++ b/src/LR/unrestricted_linear_response_D_pp.f90 @@ -39,11 +39,37 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH eF = 0d0 !----------------------------------------------- -! Build D matrix for spin-conserving transitions +! Build D matrix for spin-conserved transitions !----------------------------------------------- if(ispin == 1) then + ! 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 + + +!----------------------------------------------- +! Build D matrix for spin-flip transitions +!----------------------------------------------- + + if(ispin == 2) then + ! aaaa block ij = 0 @@ -62,6 +88,9 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH end do end do end do + end if + + if (ispin == 3) then ! bbbb block @@ -74,7 +103,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH 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) & + D_pp(ij,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 @@ -84,32 +113,4 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH 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 index 6c9019a..3111e89 100644 --- a/src/LR/unrestricted_linear_response_pp.f90 +++ b/src/LR/unrestricted_linear_response_pp.f90 @@ -55,25 +55,22 @@ 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 + + allocate(C(nPt,nPt),B(nPt,nHt),D(nHt,nHt),M(nPt+nHt,nPt+nHt),Z(nPt+nHt,nPt+nHt)& +,Omega(nPt+nHt)) +!write(*,*) 'Hello' ! 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) +e,ERI_aaaa,ERI_aabb,ERI_bbbb,D) !call matout(nHt,nHt,D) !write(*,*) 'Hello' @@ -91,22 +88,20 @@ ERI_aaaa,ERI_aabb,ERI_bbbb,D) ! 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(:,:)) - -! end if Pourquoi ne faut-il pas de end if ici ? + call sort_ppRPA(nHt,nPt,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),& +Y2(:,:)) ! 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 !!!' diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90 index 59513a4..395e078 100644 --- a/src/RPA/ppURPA.f90 +++ b/src/RPA/ppURPA.f90 @@ -51,46 +51,29 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH 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 +!spin-conserved quantities + + nPab = nV(1)*nV(2) + nHab = nO(1)*nO(2) + + nP_sc = nPab + nH_sc = nHab + +! 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)) + 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) @@ -102,6 +85,33 @@ write(*,*) 'Hello!' ispin = 2 +!spin-flip quantities + + nPaa = nV(1)*(nV(1)-1)/2 + nPbb = nV(2)*(nV(2)-1)/2 + + nP_sf = nPaa + + nHaa = nO(1)*(nO(1)-1)/2 + nHbb = nO(2)*(nO(2)-1)/2 + + nH_sf = nHaa + +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)) + +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)) + + ispin = 3 + + nP_sf = nPbb + nH_sf = nHbb + +!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)) + 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))