diff --git a/input/methods b/input/methods index 07cc57d..df59e1c 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 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW 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 bac7e06..66bba38 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 diff --git a/mol/h2.xyz b/mol/h2.xyz index 3c8e04d..21fde66 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. 0.5 diff --git a/src/LR/unrestricted_linear_response_C_pp.f90 b/src/LR/unrestricted_linear_response_C_pp.f90 index 45e8643..7e37f78 100644 --- a/src/LR/unrestricted_linear_response_C_pp.f90 +++ b/src/LR/unrestricted_linear_response_C_pp.f90 @@ -1,5 +1,5 @@ -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) +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 @@ -54,7 +54,7 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP 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,d) + lambda*ERI_aabb(a,b,c,d) + *Kronecker_delta(b,d) + lambda*ERI_aabb(a,b,c,d) end do end do end do @@ -79,9 +79,10 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP 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) + 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 @@ -103,12 +104,13 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP 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)) + *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 do + end do end if diff --git a/src/LR/unrestricted_linear_response_D_pp.f90 b/src/LR/unrestricted_linear_response_D_pp.f90 index 02ab1d0..3d7c10e 100644 --- a/src/LR/unrestricted_linear_response_D_pp.f90 +++ b/src/LR/unrestricted_linear_response_D_pp.f90 @@ -1,5 +1,5 @@ -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) +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 @@ -55,7 +55,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH 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) + *Kronecker_delta(j,l) +lambda*ERI_aabb(i,j,k,l) end do end do end do @@ -81,14 +81,15 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH 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)) + 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 - end if + end if if (ispin == 3) then @@ -104,7 +105,8 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH kl = kl + 1 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)) + *Kronecker_delta(j,l) + lambda*(ERI_bbbb(i,j,k,l) & + - ERI_bbbb(i,j,l,k)) end do end do diff --git a/src/LR/unrestricted_linear_response_pp.f90 b/src/LR/unrestricted_linear_response_pp.f90 index 3111e89..89acb7f 100644 --- a/src/LR/unrestricted_linear_response_pp.f90 +++ b/src/LR/unrestricted_linear_response_pp.f90 @@ -1,6 +1,7 @@ -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) +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 @@ -56,53 +57,47 @@ 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(*,*) 'Hello' + 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_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_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,& -e,ERI_aaaa,ERI_aabb,ERI_bbbb,D) + call unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,& + lambda,e,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) + 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 +! 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) + 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 +! 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(:,:)) ! Compute the RPA correlation energy - EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nPt,C(:,:)) - trace_matrix(nHt,D(:,:)) ) + 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 73f960f..b6de2f7 100644 --- a/src/RPA/ppURPA.f90 +++ b/src/RPA/ppURPA.f90 @@ -57,24 +57,25 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ispin = 1 -!spin-conserved quantities +!Spin-conserved quantities - nPab = nV(1)*nV(2) - nHab = nO(1)*nO(2) + nPab = nV(1)*nV(2) + nHab = nO(1)*nO(2) - nP_sc = nPab - nH_sc = nHab + 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)) + 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)) + 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)) - call print_excitation('pp-RPA (N+2)',5,nP_sc,Omega1sc) - call print_excitation('pp-RPA (N-2)',5,nH_sc,Omega2sc) + call print_excitation('pp-RPA (N+2)',5,nP_sc,Omega1sc) + call print_excitation('pp-RPA (N-2)',5,nH_sc,Omega2sc) endif @@ -84,39 +85,40 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ispin = 2 -!spin-flip quantities +!Spin-flip quantities - nPaa = nV(1)*(nV(1)-1)/2 - nPbb = nV(2)*(nV(2)-1)/2 + nPaa = nV(1)*(nV(1)-1)/2 + nPbb = nV(2)*(nV(2)-1)/2 - nP_sf = nPaa + nP_sf = nPaa - nHaa = nO(1)*(nO(1)-1)/2 - nHbb = nO(2)*(nO(2)-1)/2 + nHaa = nO(1)*(nO(1)-1)/2 + nHbb = nO(2)*(nO(2)-1)/2 - nH_sf = nHaa + 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)) + 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)) + 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 + ispin = 3 - nP_sf = nPbb - nH_sf = nHbb + 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)) + 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) + call print_excitation('pp-RPA (N+2)',6,nP_sf,Omega1sf) + call print_excitation('pp-RPA (N-2)',6,nH_sf,Omega2sf) endif