diff --git a/input/dft b/input/dft index aaee5fd..ef6a4c2 100644 --- a/input/dft +++ b/input/dft @@ -20,9 +20,9 @@ 4 # occupation numbers 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 @@ -31,9 +31,9 @@ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.0 0.0 0.0 + 1 0.0 0.0 # Ncentered ? -F +T # Parameters for CC weight-dependent exchange functional 4 -0.718713,-0.133321,0.226288,-0.250718 diff --git a/input/methods b/input/methods index 2021856..07cc57d 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - F T F F + T F F F # MP2* MP3 MP2-F12 F F F # CCD pCCD DCD CCSD CCSD(T) diff --git a/input/options b/input/options index 32acad5..90a6cbd 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability - 256 0.0000001 T 5 2 1 T 0.0 F + 256 0.0000001 T 5 2 1 F 0.0 F # MP: # CC: maxSCF thresh DIIS n_diis @@ -15,6 +15,6 @@ # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn - F F T F + T T T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/GT/Bethe_Salpeter_Tmatrix.f90 b/src/GT/Bethe_Salpeter_Tmatrix.f90 index ecda52d..99c7a4e 100644 --- a/src/GT/Bethe_Salpeter_Tmatrix.f90 +++ b/src/GT/Bethe_Salpeter_Tmatrix.f90 @@ -83,6 +83,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eT,ERI, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) +! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,Omega1s,rho1s,Omega2s,rho2s,TAs) if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,Omega1s,rho1s,Omega2s,rho2s,TBs) @@ -96,6 +98,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eT,ERI, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) +! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,Omega1t,rho1t,Omega2t,rho2t,TAt) if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,Omega1t,rho1t,Omega2t,rho2t,TBt) @@ -110,8 +114,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! Compute BSE singlet excitation energies - call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt+TAs,TBt+TBs, & - EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAs+TAt,TBs+TBt, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & diff --git a/src/GT/G0T0.f90 b/src/GT/G0T0.f90 index e75ba9d..5194ed4 100644 --- a/src/GT/G0T0.f90 +++ b/src/GT/G0T0.f90 @@ -106,8 +106,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eHF,ERI_MO, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) -! EcRPA(ispin) = 1d0*EcRPA(ispin) - call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:)) call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:)) @@ -123,9 +121,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eHF,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) -! EcRPA(ispin) = 2d0*EcRPA(ispin) -! EcRPA(ispin) = 3d0*EcRPA(ispin) - call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:)) call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:)) diff --git a/src/GT/excitation_density_Tmatrix.f90 b/src/GT/excitation_density_Tmatrix.f90 index fa58d75..6ee28f6 100644 --- a/src/GT/excitation_density_Tmatrix.f90 +++ b/src/GT/excitation_density_Tmatrix.f90 @@ -47,71 +47,49 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do p=nC+1,nBas-nR do q=nC+1,nBas-nR -! do ab=1,nVV - ab = 0 - do a=nO+1,nBas-nR - do b=a,nBas-nR - ab = ab + 1 + do ab=1,nVV cd = 0 do c=nO+1,nBas-nR -! do d=nO+1,c do d=c,nBas-nR cd = cd + 1 rho1(p,q,ab) = rho1(p,q,ab) & - + (1d0*ERI(p,q,c,d) + 0d0*ERI(p,q,d,c))*X1(cd,ab) -! + ERI(p,q,c,d)*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(c,d))) -! + (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) + + ERI(p,q,c,d)*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(c,d))) end do end do kl = 0 do k=nC+1,nO -! do l=nC+1,k do l=k,nO kl = kl + 1 rho1(p,q,ab) = rho1(p,q,ab) & - + (1d0*ERI(p,q,k,l) + 0d0*ERI(p,q,l,k))*Y1(kl,ab) -! + ERI(p,q,k,l)*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(k,l))) -! + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(k,l))) + + ERI(p,q,k,l)*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(k,l))) end do end do end do - end do -! do ij=1,nOO - ij = 0 - do i=nC+1,nO - do j=i,nO - ij = ij + 1 + do ij=1,nOO cd = 0 do c=nO+1,nBas-nR -! do d=nO+1,c do d=c,nBas-nR cd = cd + 1 rho2(p,q,ij) = rho2(p,q,ij) & - + (1d0*ERI(p,q,c,d) + 0d0*ERI(p,q,d,c))*X2(cd,ij) -! + ERI(p,q,c,d)*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(c,d))) -! + (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(c,d))) + + ERI(p,q,c,d)*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(c,d))) end do end do kl = 0 do k=nC+1,nO -! do l=nC+1,k do l=k,nO kl = kl + 1 rho2(p,q,ij) = rho2(p,q,ij) & - + (1d0*ERI(p,q,k,l) + 0d0*ERI(p,q,l,k))*Y2(kl,ij) -! + ERI(p,q,k,l)*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(k,l))) -! + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) + + ERI(p,q,k,l)*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(k,l))) end do end do end do - end do end do end do @@ -127,7 +105,11 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do p=nC+1,nBas-nR do q=nC+1,nBas-nR - do ab=1,nVV +! do ab=1,nVV + ab = 0 + do a=nO+1,nBas-nR + do b=a+1,nBas-nR + ab = ab + 1 cd = 0 do c=nO+1,nBas-nR @@ -148,8 +130,13 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r end do end do + end do - do ij=1,nOO +! do ij=1,nOO + ij = 0 + do i=nC+1,nO + do j=i+1,nO + ij = ij + 1 cd = 0 do c=nO+1,nBas-nR @@ -170,6 +157,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r end do end do + end do end do end do @@ -185,7 +173,11 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do p=nC+1,nBas-nR do q=nC+1,nBas-nR - do ab=1,nVV +! do ab=1,nVV + ab = 0 + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + ab = ab + 1 cd = 0 do c=nO+1,nBas-nR @@ -204,8 +196,13 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r end do end do + end do - do ij=1,nOO +! do ij=1,nOO + ij = 0 + do i=nC+1,nO + do j=nC+1,nO + ij = ij + 1 cd = 0 do c=nO+1,nBas-nR @@ -224,6 +221,7 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r end do end do + end do end do end do diff --git a/src/GT/static_Tmatrix_A.f90 b/src/GT/static_Tmatrix_A.f90 index ee8158d..b9fc05c 100644 --- a/src/GT/static_Tmatrix_A.f90 +++ b/src/GT/static_Tmatrix_A.f90 @@ -26,7 +26,7 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O double precision :: chi double precision :: eps - integer :: i,j,a,b,ia,jb,kl,cd + integer :: i,j,a,b,ia,jb,kl,cd,c,d ! Output variables @@ -47,13 +47,14 @@ subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,O do cd=1,nVV eps = + Omega1(cd) -! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2) +! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2) chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2) + enddo do kl=1,nOO eps = - Omega2(kl) -! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) +! chi = chi + lambda*rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2) enddo