From 196e4c917afe5ad765546e013d8b261af53f3246 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 30 Nov 2022 17:18:51 +0100 Subject: [PATCH] introduce SRG regularizer in GF, GW, and GT and add option for reg MP2 --- input/methods | 2 +- input/options | 4 +- src/GF/regularized_self_energy_GF2.f90 | 12 ++- src/GF/regularized_self_energy_GF2_diag.f90 | 12 ++- ...larized_renormalization_factor_Tmatrix.f90 | 13 ++- src/GT/regularized_self_energy_Tmatrix.f90 | 10 +-- .../regularized_self_energy_Tmatrix_diag.f90 | 6 +- src/GW/regularized_renormalization_factor.f90 | 12 ++- .../regularized_self_energy_correlation.f90 | 8 +- ...gularized_self_energy_correlation_diag.f90 | 8 +- src/MP/MP2.f90 | 79 ++++++++++--------- src/QuAcK/QuAcK.f90 | 5 +- src/QuAcK/read_options.f90 | 8 +- 13 files changed, 92 insertions(+), 87 deletions(-) diff --git a/input/methods b/input/methods index 2aa967f..c780ffc 100644 --- a/input/methods +++ b/input/methods @@ -1,7 +1,7 @@ # RHF UHF KS MOM T F F F # MP2* MP3 MP2-F12 - F F F + T F F # CCD pCCD DCD CCSD CCSD(T) F F F F F # drCCD rCCD crCCD lCCD diff --git a/input/options b/input/options index dcd942a..d18fff5 100644 --- a/input/options +++ b/input/options @@ -1,7 +1,7 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability 512 0.0000001 T 5 1 1 F 0.0 F -# MP: - +# MP: reg + F # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip diff --git a/src/GF/regularized_self_energy_GF2.f90 b/src/GF/regularized_self_energy_GF2.f90 index e4b0475..53facf8 100644 --- a/src/GF/regularized_self_energy_GF2.f90 +++ b/src/GF/regularized_self_energy_GF2.f90 @@ -37,7 +37,7 @@ subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC ! Parameters for regularized calculations ! !-----------------------------------------! - kappa = 1.1d0 + kappa = 1d0 !----------------------------------------------------! ! Compute GF2 self-energy and renormalization factor ! @@ -52,9 +52,8 @@ subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j) num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps - dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) - dfk = dfk*fk + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) SigC(p,q) = SigC(p,q) + num*fk if(p == q) Z(p) = Z(p) - num*dfk @@ -74,9 +73,8 @@ subroutine regularized_self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b) num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps - dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) - dfk = dfk*fk + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) SigC(p,q) = SigC(p,q) + num*fk if(p == q) Z(p) = Z(p) - num*dfk diff --git a/src/GF/regularized_self_energy_GF2_diag.f90 b/src/GF/regularized_self_energy_GF2_diag.f90 index 80693e8..634eb9d 100644 --- a/src/GF/regularized_self_energy_GF2_diag.f90 +++ b/src/GF/regularized_self_energy_GF2_diag.f90 @@ -37,7 +37,7 @@ subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI ! Parameters for regularized calculations ! !-----------------------------------------! - kappa = 1.1d0 + kappa = 1d0 !----------------------------------------------------! ! Compute GF2 self-energy and renormalization factor ! @@ -51,9 +51,8 @@ subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j) num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps - dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) - dfk = dfk*fk + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) SigC(p) = SigC(p) + num*fk Z(p) = Z(p) - num*dfk @@ -71,9 +70,8 @@ subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b) num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps - dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) - dfk = dfk*fk + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) SigC(p) = SigC(p) + num*fk Z(p) = Z(p) - num*dfk diff --git a/src/GT/regularized_renormalization_factor_Tmatrix.f90 b/src/GT/regularized_renormalization_factor_Tmatrix.f90 index 457b93d..68e311f 100644 --- a/src/GT/regularized_renormalization_factor_Tmatrix.f90 +++ b/src/GT/regularized_renormalization_factor_Tmatrix.f90 @@ -33,7 +33,7 @@ subroutine regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,n ! Parameters for regularized calculations ! !-----------------------------------------! - kappa = 1.1d0 + kappa = 1d0 ! Occupied part of the T-matrix self-energy @@ -43,10 +43,8 @@ subroutine regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,n eps = e(p) + e(i) - Omega1(cd) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps - dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) - dfk = dfk*fk - + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) Z(p) = Z(p) - rho1(p,i,cd)**2*dfk enddo @@ -61,9 +59,8 @@ subroutine regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,n eps = e(p) + e(nO+a) - Omega2(kl) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps - dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) - dfk = dfk*fk + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) Z(p) = Z(p) - rho2(p,nO+a,kl)**2*dfk diff --git a/src/GT/regularized_self_energy_Tmatrix.f90 b/src/GT/regularized_self_energy_Tmatrix.f90 index b9f3cd0..ee807d9 100644 --- a/src/GT/regularized_self_energy_Tmatrix.f90 +++ b/src/GT/regularized_self_energy_Tmatrix.f90 @@ -37,7 +37,7 @@ subroutine regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1 ! Parameters for regularized calculations ! !-----------------------------------------! - kappa = 1.1d0 + kappa = 1d0 !---------------------------------------------- ! Occupied part of the T-matrix self-energy @@ -49,9 +49,9 @@ subroutine regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1 do cd=1,nVV eps = e(p) + e(i) - Omega1(cd) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps - SigT(p,q) = SigT(p,q) + rho1(p,i,cd)*rho1(q,i,cd)*eps/(eps**2 + eta**2) + SigT(p,q) = SigT(p,q) + rho1(p,i,cd)*rho1(q,i,cd)*fk enddo enddo @@ -68,9 +68,9 @@ subroutine regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1 do kl=1,nOO eps = e(p) + e(a) - Omega2(kl) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps - SigT(p,q) = SigT(p,q) + rho2(p,a,kl)*rho2(q,a,kl)*eps/(eps**2 + eta**2) + SigT(p,q) = SigT(p,q) + rho2(p,a,kl)*rho2(q,a,kl)*fk enddo enddo diff --git a/src/GT/regularized_self_energy_Tmatrix_diag.f90 b/src/GT/regularized_self_energy_Tmatrix_diag.f90 index 90b5f89..978bd2a 100644 --- a/src/GT/regularized_self_energy_Tmatrix_diag.f90 +++ b/src/GT/regularized_self_energy_Tmatrix_diag.f90 @@ -37,7 +37,7 @@ subroutine regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,O ! Parameters for regularized calculations ! !-----------------------------------------! - kappa = 1.1d0 + kappa = 1d0 !---------------------------------------------- ! Occupied part of the T-matrix self-energy @@ -48,7 +48,7 @@ subroutine regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,O do cd=1,nVV eps = e(p) + e(i) - Omega1(cd) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps SigT(p) = SigT(p) + rho1(p,i,cd)**2*fk @@ -65,7 +65,7 @@ subroutine regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,O do kl=1,nOO eps = e(p) + e(a) - Omega2(kl) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps SigT(p) = SigT(p) + rho2(p,a,kl)**2*fk diff --git a/src/GW/regularized_renormalization_factor.f90 b/src/GW/regularized_renormalization_factor.f90 index dea118f..f1ddc05 100644 --- a/src/GW/regularized_renormalization_factor.f90 +++ b/src/GW/regularized_renormalization_factor.f90 @@ -39,7 +39,7 @@ subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,O ! Parameters for regularized calculations ! !-----------------------------------------! - kappa = 1.1d0 + kappa = 1d0 ! static COHSEX approximation @@ -56,9 +56,8 @@ subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,O do i=nC+1,nO do jb=1,nS eps = e(p) - e(i) + Omega(jb) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps - dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) - dfk = dfk*fk + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*dfk end do end do @@ -70,9 +69,8 @@ subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,O do a=nO+1,nBas-nR do jb=1,nS eps = e(p) - e(a) - Omega(jb) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps - dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) - dfk = dfk*fk + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*dfk end do end do diff --git a/src/GW/regularized_self_energy_correlation.f90 b/src/GW/regularized_self_energy_correlation.f90 index dcb6ffc..16808a9 100644 --- a/src/GW/regularized_self_energy_correlation.f90 +++ b/src/GW/regularized_self_energy_correlation.f90 @@ -42,7 +42,7 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e, ! Parameters for regularized MP2 calculations ! !---------------------------------------------! - kappa = 1.1d0 + kappa = 1d0 !-----------------------------! ! COHSEX static approximation ! @@ -92,7 +92,7 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e, do i=nC+1,nO do jb=1,nS eps = e(p) - e(i) + Omega(jb) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*fk end do end do @@ -106,7 +106,7 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e, do a=nO+1,nBas-nR do jb=1,nS eps = e(p) - e(a) - Omega(jb) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*fk end do end do @@ -120,7 +120,7 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e, do a=nO+1,nBas-nR do jb=1,nS eps = e(a) - e(i) + Omega(jb) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk end do end do diff --git a/src/GW/regularized_self_energy_correlation_diag.f90 b/src/GW/regularized_self_energy_correlation_diag.f90 index 2856e84..7f4ea9e 100644 --- a/src/GW/regularized_self_energy_correlation_diag.f90 +++ b/src/GW/regularized_self_energy_correlation_diag.f90 @@ -40,7 +40,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, ! Parameters for regularized calculations ! !-----------------------------------------! - kappa = 1.1d0 + kappa = 1d0 !----------------------------- ! COHSEX static self-energy @@ -87,7 +87,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, do i=nC+1,nO do jb=1,nS eps = e(p) - e(i) + Omega(jb) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*fk end do end do @@ -99,7 +99,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, do a=nO+1,nBas-nR do jb=1,nS eps = e(p) - e(a) - Omega(jb) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*fk end do end do @@ -112,7 +112,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, do a=nO+1,nBas-nR do jb=1,nS eps = e(a) - e(i) + Omega(jb) - fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk end do end do diff --git a/src/MP/MP2.f90 b/src/MP/MP2.f90 index 30e857e..dd28df8 100644 --- a/src/MP/MP2.f90 +++ b/src/MP/MP2.f90 @@ -1,4 +1,4 @@ -subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) +subroutine MP2(regularize,nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) ! Perform second-order Moller-Plesset calculation with and without regularizers @@ -6,6 +6,7 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) ! Input variables + logical,intent(in) :: regularize integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -115,55 +116,59 @@ subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2) write(*,'(A32)') '--------------------------' write(*,*) + if(regularize) then + !-------------------! ! sigma1-MP2 energy ! !-------------------! - write(*,*) - write(*,'(A32)') '--------------------------' - write(*,'(A32)') ' sigma-MP2 calculation ' - write(*,'(A32)') '--------------------------' - write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcsMP2 - write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2ds - write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs - write(*,'(A32)') '--------------------------' - write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EcsMP2 - write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EcsMP2 - write(*,'(A32)') '--------------------------' - write(*,*) + write(*,*) + write(*,'(A32)') '--------------------------' + write(*,'(A32)') ' sigma-MP2 calculation ' + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EcsMP2 + write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2ds + write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EcsMP2 + write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EcsMP2 + write(*,'(A32)') '--------------------------' + write(*,*) !--------------------! ! sigma^2-MP2 energy ! !--------------------! - write(*,*) - write(*,'(A32)') '--------------------------' - write(*,'(A32)') ' sigma^2-MP2 calculation ' - write(*,'(A32)') '--------------------------' - write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',Ecs2MP2 - write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2ds2 - write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs2 - write(*,'(A32)') '--------------------------' - write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + Ecs2MP2 - write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + Ecs2MP2 - write(*,'(A32)') '--------------------------' - write(*,*) + write(*,*) + write(*,'(A32)') '--------------------------' + write(*,'(A32)') ' sigma^2-MP2 calculation ' + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',Ecs2MP2 + write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2ds2 + write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xs2 + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + Ecs2MP2 + write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + Ecs2MP2 + write(*,'(A32)') '--------------------------' + write(*,*) !------------------! ! kappa-MP2 energy ! !------------------! - write(*,*) - write(*,'(A32)') '--------------------------' - write(*,'(A32)') ' kappa-MP2 calculation ' - write(*,'(A32)') '--------------------------' - write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EckMP2 - write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2dk - write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xk - write(*,'(A32)') '--------------------------' - write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EckMP2 - write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EckMP2 - write(*,'(A32)') '--------------------------' - write(*,*) + write(*,*) + write(*,'(A32)') '--------------------------' + write(*,'(A32)') ' kappa-MP2 calculation ' + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',EckMP2 + write(*,'(A32,1X,F16.10)') ' Direct part = ',2d0*E2dk + write(*,'(A32,1X,F16.10)') ' Exchange part = ',-E2xk + write(*,'(A32)') '--------------------------' + write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ',EHF + EckMP2 + write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + EckMP2 + write(*,'(A32)') '--------------------------' + write(*,*) + + end if end subroutine MP2 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index b1ae456..ccda5a9 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -108,6 +108,8 @@ program QuAcK double precision :: thresh_HF,level_shift logical :: DIIS_HF,guess_type,ortho_type,mix + logical :: regMP + integer :: maxSCF_CC,n_diis_CC double precision :: thresh_CC logical :: DIIS_CC @@ -174,6 +176,7 @@ program QuAcK ! Read options for methods call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, & + regMP, & maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & TDA,singlet,triplet,spin_conserved,spin_flip, & maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, & @@ -479,7 +482,7 @@ program QuAcK else - call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,EcMP2) + call MP2(regMP,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,EcMP2) end if diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 2138b05..5c93fe0 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -1,4 +1,5 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, & + regMP, & maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & TDA,singlet,triplet,spin_conserved,spin_flip, & maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, & @@ -24,6 +25,8 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t double precision,intent(out) :: level_shift logical,intent(out) :: dostab + logical,intent(out) :: regMP + integer,intent(out) :: maxSCF_CC double precision,intent(out) :: thresh_CC logical,intent(out) :: DIIS_CC @@ -106,8 +109,11 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t ! Read MPn options - read(1,*) + regMP = .false. read(1,*) + read(1,*) answer1 + + if(answer1 == 'T') regMP = .true. ! Read CC options