diff --git a/input/methods b/input/methods index 688d56a..9b32026 100644 --- a/input/methods +++ b/input/methods @@ -3,7 +3,7 @@ # MP2* MP3 MP2-F12 F F F # CCD pCCD DCD CCSD CCSD(T) - F F F F F + F F F T F # drCCD rCCD crCCD lCCD F F F F # CIS* CIS(D) CID CISD FCI diff --git a/input/options b/input/options index 6d53fb5..2dc2665 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 T F + 256 0.0000001 T 5 2 1 T 1.0 F # MP: # CC: maxSCF thresh DIIS n_diis @@ -13,8 +13,8 @@ # 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 + F T T # BSE: BSE dBSE dTDA evDyn - T T T F + T F T 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 a4e936a..3a1f2fe 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 2.000000 +H 0. 0. 0.741 diff --git a/src/GT/G0T0.f90 b/src/GT/G0T0.f90 index 11f4415..e75ba9d 100644 --- a/src/GT/G0T0.f90 +++ b/src/GT/G0T0.f90 @@ -238,11 +238,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing write(*,*)'-------------------------------------------------------------------------------' write(*,*) -! Free memory - - deallocate(Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s, & - Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t) - ! Compute the BSE correlation energy via the adiabatic connection if(doACFDT) then @@ -260,7 +255,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing end if call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - ERI_MO,eHF,eG0T0,EcAC) + nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, & + Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eHF,eG0T0,EcAC) if(exchange_kernel) then diff --git a/src/GT/evGT.f90 b/src/GT/evGT.f90 index 1c5e518..2af3bdd 100644 --- a/src/GT/evGT.f90 +++ b/src/GT/evGT.f90 @@ -265,7 +265,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & end if call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - ERI_MO,eGT,eGT,EcAC) + nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, & + Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC) if(exchange_kernel) then diff --git a/src/GT/qsGT.f90 b/src/GT/qsGT.f90 index 1e4c349..4344878 100644 --- a/src/GT/qsGT.f90 +++ b/src/GT/qsGT.f90 @@ -389,7 +389,8 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T end if call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - ERI_MO,eGT,eGT,EcAC) + nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, & + Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 6cba349..0e3989e 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -12,7 +12,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc integer,intent(in) :: max_diis integer,intent(in) :: guess_type double precision,intent(in) :: thresh - logical,intent(in) :: level_shift + double precision,intent(in) :: level_shift integer,intent(in) :: nBas integer,intent(in) :: nO @@ -135,7 +135,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc ! Level-shifting - if(level_shift .and. Conv > thresh) call level_shifting(nBas,nO,S,c,F) + if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,c,F) ! Diagonalize Fock matrix diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 7917d00..39bbaf6 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -12,7 +12,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, integer,intent(in) :: max_diis integer,intent(in) :: guess_type logical,intent(in) :: mix - logical,intent(in) :: level_shift + double precision,intent(in) :: level_shift double precision,intent(in) :: thresh integer,intent(in) :: nBas @@ -187,10 +187,10 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, ! Level-shifting - if(level_shift .and. Conv > thresh) then + if(level_shift > 0d0 .and. Conv > thresh) then do ispin=1,nspin - call level_shifting(nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin)) + call level_shifting(level_shift,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin)) end do end if diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 890b48f..40d1415 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -111,8 +111,8 @@ program QuAcK double precision :: start_Bas ,end_Bas ,t_Bas integer :: maxSCF_HF,n_diis_HF - double precision :: thresh_HF - logical :: DIIS_HF,guess_type,ortho_type,mix,level_shift + double precision :: thresh_HF,level_shift + logical :: DIIS_HF,guess_type,ortho_type,mix integer :: maxSCF_CC,n_diis_CC double precision :: thresh_CC diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 6393922..8a18c34 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -22,7 +22,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t integer,intent(out) :: guess_type integer,intent(out) :: ortho_type logical,intent(out) :: mix - logical,intent(out) :: level_shift + double precision,intent(out) :: level_shift logical,intent(out) :: dostab integer,intent(out) :: maxSCF_CC @@ -101,16 +101,15 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t guess_type = 1 ortho_type = 1 mix = .false. - level_shift = .false. + level_shift = 0d0 dostab = .false. read(1,*) - read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,answer3,answer4 + read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,level_shift,answer3 if(answer1 == 'T') DIIS_HF = .true. if(answer2 == 'T') mix = .true. - if(answer3 == 'T') level_shift = .true. - if(answer4 == 'T') dostab = .true. + if(answer3 == 'T') dostab = .true. if(.not.DIIS_HF) n_diis_HF = 1 diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90 index 9f48261..051deab 100644 --- a/src/RPA/ACFDT_Tmatrix.f90 +++ b/src/RPA/ACFDT_Tmatrix.f90 @@ -1,5 +1,6 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - ERI,eT,eGT,EcAC) + nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, & + Omega2t,X2t,Y2t,rho1t,rho2t,ERI,eT,eGT,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix @@ -26,6 +27,28 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple integer,intent(in) :: nR integer,intent(in) :: nS + integer,intent(in) :: nOOs + integer,intent(in) :: nOOt + integer,intent(in) :: nVVs + integer,intent(in) :: nVVt + + double precision,intent(in) :: Omega1s(nVVs) + double precision,intent(in) :: X1s(nVVs,nVVs) + double precision,intent(in) :: Y1s(nOOs,nVVs) + double precision,intent(in) :: Omega2s(nOOs) + double precision,intent(in) :: X2s(nVVs,nOOs) + double precision,intent(in) :: Y2s(nOOs,nOOs) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs) + double precision,intent(in) :: Omega1t(nVVt) + double precision,intent(in) :: X1t(nVVt,nVVt) + double precision,intent(in) :: Y1t(nOOt,nVVt) + double precision,intent(in) :: Omega2t(nOOt) + double precision,intent(in) :: X2t(nVVt,nOOt) + double precision,intent(in) :: Y2t(nOOt,nOOt) + double precision,intent(in) :: rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: rho2t(nBas,nBas,nOOt) + double precision,intent(in) :: eT(nBas) double precision,intent(in) :: eGT(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -39,46 +62,23 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple double precision :: lambda double precision,allocatable :: Ec(:,:) - integer :: nOOs,nOOt - integer :: nVVs,nVVt - double precision :: EcRPA(nspin) - double precision,allocatable :: TA(:,:) - double precision,allocatable :: TB(:,:) + double precision,allocatable :: TAs(:,:) + double precision,allocatable :: TBs(:,:) + double precision,allocatable :: TAt(:,:) + double precision,allocatable :: TBt(:,:) double precision,allocatable :: Omega(:,:) double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) - double precision,allocatable :: Omega1s(:),Omega1t(:) - double precision,allocatable :: X1s(:,:),X1t(:,:) - double precision,allocatable :: Y1s(:,:),Y1t(:,:) - double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:) - double precision,allocatable :: Omega2s(:),Omega2t(:) - double precision,allocatable :: X2s(:,:),X2t(:,:) - double precision,allocatable :: Y2s(:,:),Y2t(:,:) - double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) - ! Output variables double precision,intent(out) :: EcAC(nspin) -! Useful quantities - - nOOs = nO*nO - nVVs = nV*nV - - nOOt = nO*(nO-1)/2 - nVVt = nV*(nV-1)/2 - ! Memory allocation - allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & - Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & - rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & - Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & - Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & - rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) - allocate(TA(nS,nS),TB(nS,nS),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + allocate(TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS), & + Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) allocate(Ec(nAC,nspin)) ! Antisymmetrized kernel version @@ -113,11 +113,6 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple lambda = rAC(iAC) - ! Initialize T matrix - - TA(:,:) = 0d0 - TB(:,:) = 0d0 - if(doXBS) then isp_T = 1 @@ -128,8 +123,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple 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,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TAs) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TBs) isp_T = 2 iblock = 4 @@ -139,12 +134,12 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple 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,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TAt) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TBt) end if - call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & + call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt+TAs,TBt+TBs, & EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) @@ -183,11 +178,6 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple lambda = rAC(iAC) - ! Initialize T matrix - - TA(:,:) = 0d0 - TB(:,:) = 0d0 - if(doXBS) then isp_T = 1 @@ -198,8 +188,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple 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,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TAs) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TBs) isp_T = 2 iblock = 4 @@ -209,12 +199,12 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple 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,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TAt) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TBt) end if - call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & + call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt-TAs,TBt-TBs, & EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) diff --git a/src/eDFT/UKS.f90 b/src/eDFT/UKS.f90 index 3d55bf1..152b35f 100644 --- a/src/eDFT/UKS.f90 +++ b/src/eDFT/UKS.f90 @@ -20,7 +20,7 @@ subroutine UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,t integer,intent(in) :: max_diis integer,intent(in) :: guess_type logical,intent(in) :: mix - logical,intent(in) :: level_shift + double precision,intent(in) :: level_shift double precision,intent(in) :: thresh integer,intent(in) :: nBas double precision,intent(in) :: AO(nBas,nGrid) @@ -275,10 +275,10 @@ subroutine UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,t ! Level-shifting - if(level_shift .and. Conv > thresh) then + if(level_shift > 0d0 .and. Conv > thresh) then do ispin=1,nspin - call level_shifting(nBas,maxval(nO(ispin,:)),S,c,F(:,:,ispin)) + call level_shifting(level_shift,nBas,maxval(nO(ispin,:)),S,c,F(:,:,ispin)) end do end if diff --git a/src/utils/level_shifting.f90 b/src/utils/level_shifting.f90 index 022b524..9e8e23c 100644 --- a/src/utils/level_shifting.f90 +++ b/src/utils/level_shifting.f90 @@ -1,4 +1,4 @@ -subroutine level_shifting(nBas,nO,S,c,F) +subroutine level_shifting(level_shift,nBas,nO,S,c,F) ! Perform level-shifting on the Fock matrix @@ -6,6 +6,7 @@ subroutine level_shifting(nBas,nO,S,c,F) ! Input variables + double precision,intent(in) :: level_shift integer,intent(in) :: nBas integer,intent(in) :: nO double precision,intent(in) :: S(nBas,nBas) @@ -15,7 +16,6 @@ subroutine level_shifting(nBas,nO,S,c,F) double precision,allocatable :: F_MO(:,:) double precision,allocatable :: Sc(:,:) - double precision,parameter :: LS = 0.1d0 integer :: a @@ -28,7 +28,7 @@ subroutine level_shifting(nBas,nO,S,c,F) F_MO(:,:) = matmul(transpose(c),matmul(F,c)) do a=nO+1,nBas - F_MO(a,a) = F_MO(a,a) + LS + F_MO(a,a) = F_MO(a,a) + level_shift end do Sc(:,:) = matmul(S,c)