From 203ce9541cd05b80740028555af5728c0faf3c83 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 17 Dec 2021 11:41:40 +0100 Subject: [PATCH] regularized GW --- input/methods | 2 +- input/options | 14 +- mol/h2.xyz | 2 +- src/GF/G0F2.f90 | 4 +- src/GF/UG0F2.f90 | 3 +- src/GF/evGF2.f90 | 4 +- src/GF/evUGF2.f90 | 3 +- src/GF/qsGF2.f90 | 3 +- src/GF/qsUGF2.f90 | 3 +- src/{MBPT => GT}/Bethe_Salpeter_Tmatrix.f90 | 0 ..._Salpeter_Tmatrix_dynamic_perturbation.f90 | 0 src/{MBPT => GT}/G0T0.f90 | 3 +- src/{MBPT => GT}/dynamic_Tmatrix_A.f90 | 0 src/{MBPT => GT}/evGT.f90 | 3 +- .../excitation_density_Tmatrix.f90 | 0 src/{MBPT => GT}/print_G0T0.f90 | 0 src/{MBPT => GT}/print_evGT.f90 | 0 src/{MBPT => GT}/print_qsGT.f90 | 0 src/{MBPT => GT}/qsGT.f90 | 3 +- .../renormalization_factor_Tmatrix.f90 | 0 src/{MBPT => GT}/self_energy_Tmatrix.f90 | 0 src/{MBPT => GT}/self_energy_Tmatrix_diag.f90 | 0 src/{MBPT => GT}/static_Tmatrix_TA.f90 | 2 +- src/{MBPT => GT}/static_Tmatrix_TB.f90 | 2 +- src/{MBPT => GW}/Bethe_Salpeter.f90 | 0 .../Bethe_Salpeter_AB_matrix_dynamic.f90 | 0 src/{MBPT => GW}/Bethe_Salpeter_A_matrix.f90 | 0 .../Bethe_Salpeter_A_matrix_dynamic.f90 | 0 src/{MBPT => GW}/Bethe_Salpeter_B_matrix.f90 | 0 .../Bethe_Salpeter_ZAB_matrix_dynamic.f90 | 0 .../Bethe_Salpeter_ZA_matrix_dynamic.f90 | 0 .../Bethe_Salpeter_dynamic_perturbation.f90 | 0 ...alpeter_dynamic_perturbation_iterative.f90 | 0 src/{MBPT => GW}/G0W0.f90 | 18 ++- src/{MBPT => GW}/G0W0_SOSEX.f90 | 0 src/{MBPT => GW}/QP_graph.f90 | 0 src/{MBPT => GW}/QP_roots.f90 | 0 .../Sangalli_dynamic_perturbation.f90 | 0 src/{MBPT => GW}/SigmaC.f90 | 0 src/{MBPT => GW}/UG0W0.f90 | 3 +- src/{MBPT => GW}/USigmaC.f90 | 0 src/{MBPT => GW}/dSigmaC.f90 | 0 src/{MBPT => GW}/dUSigmaC.f90 | 0 src/{MBPT => GW}/evGW.f90 | 46 ++++-- src/{MBPT => GW}/evUGW.f90 | 3 +- src/{MBPT => GW}/exchange_matrix_MO_basis.f90 | 0 src/{MBPT => GW}/excitation_density.f90 | 0 src/{MBPT => GW}/excitation_density_SOSEX.f90 | 0 src/{MBPT => GW}/plot_GW.f90 | 0 src/{MBPT => GW}/print_G0W0.f90 | 0 src/{MBPT => GW}/print_SOSEX.f90 | 0 src/{MBPT => GW}/print_UG0W0.f90 | 0 src/{MBPT => GW}/print_evGW.f90 | 0 src/{MBPT => GW}/print_evUGW.f90 | 0 src/{MBPT => GW}/print_qsGW.f90 | 0 src/{MBPT => GW}/print_qsUGW.f90 | 0 src/{MBPT => GW}/qsGW.f90 | 31 +++- src/{MBPT => GW}/qsGW_PT.f90 | 0 src/{MBPT => GW}/qsUGW.f90 | 3 +- src/GW/regularized_renormalization_factor.f90 | 87 ++++++++++++ .../regularized_self_energy_correlation.f90 | 7 +- ...gularized_self_energy_correlation_diag.f90 | 7 +- src/{MBPT => GW}/renormalization_factor.f90 | 33 ++--- .../renormalization_factor_SOSEX.f90 | 0 src/{MBPT => GW}/self_energy_correlation.f90 | 7 +- .../self_energy_correlation_SOSEX_diag.f90 | 0 .../self_energy_correlation_diag.f90 | 1 - src/{MBPT => GW}/self_energy_exchange.f90 | 0 .../self_energy_exchange_diag.f90 | 0 src/{MBPT => GW}/static_screening_WA.f90 | 0 src/{MBPT => GW}/static_screening_WB.f90 | 0 src/{MBPT => GW}/ufBSE.f90 | 3 +- src/{MBPT => GW}/ufG0W0.f90 | 3 +- src/{MBPT => GW}/ufGW.f90 | 3 +- .../unrestricted_Bethe_Salpeter.f90 | 0 .../unrestricted_Bethe_Salpeter_A_matrix.f90 | 0 ...ricted_Bethe_Salpeter_A_matrix_dynamic.f90 | 0 .../unrestricted_Bethe_Salpeter_B_matrix.f90 | 0 ...icted_Bethe_Salpeter_ZA_matrix_dynamic.f90 | 0 ...ed_Bethe_Salpeter_dynamic_perturbation.f90 | 0 src/{MBPT => GW}/unrestricted_QP_graph.f90 | 0 .../unrestricted_excitation_density.f90 | 0 ...ted_regularized_renormalization_factor.f90 | 111 +++++++++++++++ ...ed_regularized_self_energy_correlation.f90 | 133 ++++++++++++++++++ ...gularized_self_energy_correlation_diag.f90 | 126 +++++++++++++++++ .../unrestricted_renormalization_factor.f90 | 2 +- .../unrestricted_self_energy_correlation.f90 | 2 +- ...estricted_self_energy_correlation_diag.f90 | 2 +- src/MBPT/Makefile | 10 -- src/MBPT/Sangalli_A_matrix_dynamic.f90.x | 78 ---------- src/MBPT/excitation_density_RI.f90 | 65 --------- src/MBPT/excitation_density_SOSEX_RI.f90 | 65 --------- src/QuAcK/QuAcK.f90 | 66 +++++---- src/QuAcK/read_options.f90 | 49 ++++++- 94 files changed, 679 insertions(+), 336 deletions(-) rename src/{MBPT => GT}/Bethe_Salpeter_Tmatrix.f90 (100%) rename src/{MBPT => GT}/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 (100%) rename src/{MBPT => GT}/G0T0.f90 (98%) rename src/{MBPT => GT}/dynamic_Tmatrix_A.f90 (100%) rename src/{MBPT => GT}/evGT.f90 (99%) rename src/{MBPT => GT}/excitation_density_Tmatrix.f90 (100%) rename src/{MBPT => GT}/print_G0T0.f90 (100%) rename src/{MBPT => GT}/print_evGT.f90 (100%) rename src/{MBPT => GT}/print_qsGT.f90 (100%) rename src/{MBPT => GT}/qsGT.f90 (98%) rename src/{MBPT => GT}/renormalization_factor_Tmatrix.f90 (100%) rename src/{MBPT => GT}/self_energy_Tmatrix.f90 (100%) rename src/{MBPT => GT}/self_energy_Tmatrix_diag.f90 (100%) rename src/{MBPT => GT}/static_Tmatrix_TA.f90 (96%) rename src/{MBPT => GT}/static_Tmatrix_TB.f90 (96%) rename src/{MBPT => GW}/Bethe_Salpeter.f90 (100%) rename src/{MBPT => GW}/Bethe_Salpeter_AB_matrix_dynamic.f90 (100%) rename src/{MBPT => GW}/Bethe_Salpeter_A_matrix.f90 (100%) rename src/{MBPT => GW}/Bethe_Salpeter_A_matrix_dynamic.f90 (100%) rename src/{MBPT => GW}/Bethe_Salpeter_B_matrix.f90 (100%) rename src/{MBPT => GW}/Bethe_Salpeter_ZAB_matrix_dynamic.f90 (100%) rename src/{MBPT => GW}/Bethe_Salpeter_ZA_matrix_dynamic.f90 (100%) rename src/{MBPT => GW}/Bethe_Salpeter_dynamic_perturbation.f90 (100%) rename src/{MBPT => GW}/Bethe_Salpeter_dynamic_perturbation_iterative.f90 (100%) rename src/{MBPT => GW}/G0W0.f90 (92%) rename src/{MBPT => GW}/G0W0_SOSEX.f90 (100%) rename src/{MBPT => GW}/QP_graph.f90 (100%) rename src/{MBPT => GW}/QP_roots.f90 (100%) rename src/{MBPT => GW}/Sangalli_dynamic_perturbation.f90 (100%) rename src/{MBPT => GW}/SigmaC.f90 (100%) rename src/{MBPT => GW}/UG0W0.f90 (98%) rename src/{MBPT => GW}/USigmaC.f90 (100%) rename src/{MBPT => GW}/dSigmaC.f90 (100%) rename src/{MBPT => GW}/dUSigmaC.f90 (100%) rename src/{MBPT => GW}/evGW.f90 (87%) rename src/{MBPT => GW}/evUGW.f90 (99%) rename src/{MBPT => GW}/exchange_matrix_MO_basis.f90 (100%) rename src/{MBPT => GW}/excitation_density.f90 (100%) rename src/{MBPT => GW}/excitation_density_SOSEX.f90 (100%) rename src/{MBPT => GW}/plot_GW.f90 (100%) rename src/{MBPT => GW}/print_G0W0.f90 (100%) rename src/{MBPT => GW}/print_SOSEX.f90 (100%) rename src/{MBPT => GW}/print_UG0W0.f90 (100%) rename src/{MBPT => GW}/print_evGW.f90 (100%) rename src/{MBPT => GW}/print_evUGW.f90 (100%) rename src/{MBPT => GW}/print_qsGW.f90 (100%) rename src/{MBPT => GW}/print_qsUGW.f90 (100%) rename src/{MBPT => GW}/qsGW.f90 (91%) rename src/{MBPT => GW}/qsGW_PT.f90 (100%) rename src/{MBPT => GW}/qsUGW.f90 (99%) create mode 100644 src/GW/regularized_renormalization_factor.f90 rename src/{MBPT => GW}/regularized_self_energy_correlation.f90 (92%) rename src/{MBPT => GW}/regularized_self_energy_correlation_diag.f90 (93%) rename src/{MBPT => GW}/renormalization_factor.f90 (63%) rename src/{MBPT => GW}/renormalization_factor_SOSEX.f90 (100%) rename src/{MBPT => GW}/self_energy_correlation.f90 (91%) rename src/{MBPT => GW}/self_energy_correlation_SOSEX_diag.f90 (100%) rename src/{MBPT => GW}/self_energy_correlation_diag.f90 (98%) rename src/{MBPT => GW}/self_energy_exchange.f90 (100%) rename src/{MBPT => GW}/self_energy_exchange_diag.f90 (100%) rename src/{MBPT => GW}/static_screening_WA.f90 (100%) rename src/{MBPT => GW}/static_screening_WB.f90 (100%) rename src/{MBPT => GW}/ufBSE.f90 (97%) rename src/{MBPT => GW}/ufG0W0.f90 (97%) rename src/{MBPT => GW}/ufGW.f90 (97%) rename src/{MBPT => GW}/unrestricted_Bethe_Salpeter.f90 (100%) rename src/{MBPT => GW}/unrestricted_Bethe_Salpeter_A_matrix.f90 (100%) rename src/{MBPT => GW}/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 (100%) rename src/{MBPT => GW}/unrestricted_Bethe_Salpeter_B_matrix.f90 (100%) rename src/{MBPT => GW}/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 (100%) rename src/{MBPT => GW}/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 (100%) rename src/{MBPT => GW}/unrestricted_QP_graph.f90 (100%) rename src/{MBPT => GW}/unrestricted_excitation_density.f90 (100%) create mode 100644 src/GW/unrestricted_regularized_renormalization_factor.f90 create mode 100644 src/GW/unrestricted_regularized_self_energy_correlation.f90 create mode 100644 src/GW/unrestricted_regularized_self_energy_correlation_diag.f90 rename src/{MBPT => GW}/unrestricted_renormalization_factor.f90 (97%) rename src/{MBPT => GW}/unrestricted_self_energy_correlation.f90 (98%) rename src/{MBPT => GW}/unrestricted_self_energy_correlation_diag.f90 (98%) delete mode 100644 src/MBPT/Makefile delete mode 100644 src/MBPT/Sangalli_A_matrix_dynamic.f90.x delete mode 100644 src/MBPT/excitation_density_RI.f90 delete mode 100644 src/MBPT/excitation_density_SOSEX_RI.f90 diff --git a/input/methods b/input/methods index 02ac24b..df2ad59 100644 --- a/input/methods +++ b/input/methods @@ -13,7 +13,7 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - F F T F F + F T F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index f8f9987..e8d2643 100644 --- a/input/options +++ b/input/options @@ -1,15 +1,17 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability - 1024 0.00001 T 5 1 1 T F + 1024 0.00001 T 5 1 1 T F # MP: # CC: maxSCF thresh DIIS n_diis 64 0.00001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T T T T -# GF: maxSCF thresh DIIS n_diis lin eta renorm - 256 0.00001 T 5 T 0.00367493 3 -# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 - 256 0.00001 T 5 T 0.0 F F F F F + F T T T T +# 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 +# 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 diff --git a/mol/h2.xyz b/mol/h2.xyz index 21fde66..d955cc4 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 0.5 +H 0. 0. 0.7 diff --git a/src/GF/G0F2.f90 b/src/GF/G0F2.f90 index cfc7750..2e98066 100644 --- a/src/GF/G0F2.f90 +++ b/src/GF/G0F2.f90 @@ -1,4 +1,5 @@ -subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) +subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform a one-shot second-order Green function calculation @@ -16,6 +17,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nBas integer,intent(in) :: nO integer,intent(in) :: nC diff --git a/src/GF/UG0F2.f90 b/src/GF/UG0F2.f90 index 0389f93..be7d4c1 100644 --- a/src/GF/UG0F2.f90 +++ b/src/GF/UG0F2.f90 @@ -1,4 +1,4 @@ -subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & +subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF) ! Perform unrestricted G0W0 calculation @@ -18,6 +18,7 @@ subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta, logical,intent(in) :: spin_flip logical,intent(in) :: linearize double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) diff --git a/src/GF/evGF2.f90 b/src/GF/evGF2.f90 index 2203533..f55ae32 100644 --- a/src/GF/evGF2.f90 +++ b/src/GF/evGF2.f90 @@ -1,5 +1,5 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, & - linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform eigenvalue self-consistent second-order Green function calculation @@ -20,6 +20,8 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta + logical,intent(in) :: regularize + integer,intent(in) :: nBas integer,intent(in) :: nO integer,intent(in) :: nC diff --git a/src/GF/evUGF2.f90 b/src/GF/evUGF2.f90 index 5f4d0f3..07749af 100644 --- a/src/GF/evUGF2.f90 +++ b/src/GF/evUGF2.f90 @@ -1,5 +1,5 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & - eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & + eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,cHF,eHF) ! Perform self-consistent eigenvalue-only GW calculation @@ -22,6 +22,7 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index db8e50d..3019091 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -1,5 +1,5 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, & - eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & + eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GF2 calculation @@ -20,6 +20,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 3295f6f..71d3c97 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -1,4 +1,4 @@ -subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta, & +subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize, & nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO, & ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) @@ -20,6 +20,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) diff --git a/src/MBPT/Bethe_Salpeter_Tmatrix.f90 b/src/GT/Bethe_Salpeter_Tmatrix.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter_Tmatrix.f90 rename to src/GT/Bethe_Salpeter_Tmatrix.f90 diff --git a/src/MBPT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 b/src/GT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 rename to src/GT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 diff --git a/src/MBPT/G0T0.f90 b/src/GT/G0T0.f90 similarity index 98% rename from src/MBPT/G0T0.f90 rename to src/GT/G0T0.f90 index 3054daa..8d058e7 100644 --- a/src/MBPT/G0T0.f90 +++ b/src/GT/G0T0.f90 @@ -1,5 +1,5 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet, & - linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0) + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -21,6 +21,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nBas integer,intent(in) :: nC diff --git a/src/MBPT/dynamic_Tmatrix_A.f90 b/src/GT/dynamic_Tmatrix_A.f90 similarity index 100% rename from src/MBPT/dynamic_Tmatrix_A.f90 rename to src/GT/dynamic_Tmatrix_A.f90 diff --git a/src/MBPT/evGT.f90 b/src/GT/evGT.f90 similarity index 99% rename from src/MBPT/evGT.f90 rename to src/GT/evGT.f90 index 12bf2a5..65dcfa4 100644 --- a/src/MBPT/evGT.f90 +++ b/src/GT/evGT.f90 @@ -1,5 +1,5 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & - BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas, & + BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nBas, & nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0T0) ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) @@ -24,6 +24,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nBas integer,intent(in) :: nC diff --git a/src/MBPT/excitation_density_Tmatrix.f90 b/src/GT/excitation_density_Tmatrix.f90 similarity index 100% rename from src/MBPT/excitation_density_Tmatrix.f90 rename to src/GT/excitation_density_Tmatrix.f90 diff --git a/src/MBPT/print_G0T0.f90 b/src/GT/print_G0T0.f90 similarity index 100% rename from src/MBPT/print_G0T0.f90 rename to src/GT/print_G0T0.f90 diff --git a/src/MBPT/print_evGT.f90 b/src/GT/print_evGT.f90 similarity index 100% rename from src/MBPT/print_evGT.f90 rename to src/GT/print_evGT.f90 diff --git a/src/MBPT/print_qsGT.f90 b/src/GT/print_qsGT.f90 similarity index 100% rename from src/MBPT/print_qsGT.f90 rename to src/GT/print_qsGT.f90 diff --git a/src/MBPT/qsGT.f90 b/src/GT/qsGT.f90 similarity index 98% rename from src/MBPT/qsGT.f90 rename to src/GT/qsGT.f90 index 1d8e32e..6a5aa5a 100644 --- a/src/MBPT/qsGT.f90 +++ b/src/GT/qsGT.f90 @@ -1,5 +1,5 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA, & - dBSE,dTDA,evDyn,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & + dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GT calculation @@ -24,6 +24,7 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) diff --git a/src/MBPT/renormalization_factor_Tmatrix.f90 b/src/GT/renormalization_factor_Tmatrix.f90 similarity index 100% rename from src/MBPT/renormalization_factor_Tmatrix.f90 rename to src/GT/renormalization_factor_Tmatrix.f90 diff --git a/src/MBPT/self_energy_Tmatrix.f90 b/src/GT/self_energy_Tmatrix.f90 similarity index 100% rename from src/MBPT/self_energy_Tmatrix.f90 rename to src/GT/self_energy_Tmatrix.f90 diff --git a/src/MBPT/self_energy_Tmatrix_diag.f90 b/src/GT/self_energy_Tmatrix_diag.f90 similarity index 100% rename from src/MBPT/self_energy_Tmatrix_diag.f90 rename to src/GT/self_energy_Tmatrix_diag.f90 diff --git a/src/MBPT/static_Tmatrix_TA.f90 b/src/GT/static_Tmatrix_TA.f90 similarity index 96% rename from src/MBPT/static_Tmatrix_TA.f90 rename to src/GT/static_Tmatrix_TA.f90 index 1c78e0d..9935919 100644 --- a/src/MBPT/static_Tmatrix_TA.f90 +++ b/src/GT/static_Tmatrix_TA.f90 @@ -50,7 +50,7 @@ subroutine static_Tmatrix_TA(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r do kl=1,nOO ! chi = chi + lambda*rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) - chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) + chi = chi - rho2(i,j,kl)*rho2(a,b,kl)*Omega2(kl)/(Omega2(kl)**2 + eta**2) enddo TA(ia,jb) = TA(ia,jb) - 1d0*lambda*chi diff --git a/src/MBPT/static_Tmatrix_TB.f90 b/src/GT/static_Tmatrix_TB.f90 similarity index 96% rename from src/MBPT/static_Tmatrix_TB.f90 rename to src/GT/static_Tmatrix_TB.f90 index d1f993c..d4707c8 100644 --- a/src/MBPT/static_Tmatrix_TB.f90 +++ b/src/GT/static_Tmatrix_TB.f90 @@ -50,7 +50,7 @@ subroutine static_Tmatrix_TB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,r do kl=1,nOO ! chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 - chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 + chi = chi - rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 enddo TB(ia,jb) = TB(ia,jb) - 1d0*lambda*chi diff --git a/src/MBPT/Bethe_Salpeter.f90 b/src/GW/Bethe_Salpeter.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter.f90 rename to src/GW/Bethe_Salpeter.f90 diff --git a/src/MBPT/Bethe_Salpeter_AB_matrix_dynamic.f90 b/src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter_AB_matrix_dynamic.f90 rename to src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90 diff --git a/src/MBPT/Bethe_Salpeter_A_matrix.f90 b/src/GW/Bethe_Salpeter_A_matrix.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter_A_matrix.f90 rename to src/GW/Bethe_Salpeter_A_matrix.f90 diff --git a/src/MBPT/Bethe_Salpeter_A_matrix_dynamic.f90 b/src/GW/Bethe_Salpeter_A_matrix_dynamic.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter_A_matrix_dynamic.f90 rename to src/GW/Bethe_Salpeter_A_matrix_dynamic.f90 diff --git a/src/MBPT/Bethe_Salpeter_B_matrix.f90 b/src/GW/Bethe_Salpeter_B_matrix.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter_B_matrix.f90 rename to src/GW/Bethe_Salpeter_B_matrix.f90 diff --git a/src/MBPT/Bethe_Salpeter_ZAB_matrix_dynamic.f90 b/src/GW/Bethe_Salpeter_ZAB_matrix_dynamic.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter_ZAB_matrix_dynamic.f90 rename to src/GW/Bethe_Salpeter_ZAB_matrix_dynamic.f90 diff --git a/src/MBPT/Bethe_Salpeter_ZA_matrix_dynamic.f90 b/src/GW/Bethe_Salpeter_ZA_matrix_dynamic.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter_ZA_matrix_dynamic.f90 rename to src/GW/Bethe_Salpeter_ZA_matrix_dynamic.f90 diff --git a/src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90 b/src/GW/Bethe_Salpeter_dynamic_perturbation.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter_dynamic_perturbation.f90 rename to src/GW/Bethe_Salpeter_dynamic_perturbation.f90 diff --git a/src/MBPT/Bethe_Salpeter_dynamic_perturbation_iterative.f90 b/src/GW/Bethe_Salpeter_dynamic_perturbation_iterative.f90 similarity index 100% rename from src/MBPT/Bethe_Salpeter_dynamic_perturbation_iterative.f90 rename to src/GW/Bethe_Salpeter_dynamic_perturbation_iterative.f90 diff --git a/src/MBPT/G0W0.f90 b/src/GW/G0W0.f90 similarity index 92% rename from src/MBPT/G0W0.f90 rename to src/GW/G0W0.f90 index 9065661..aa35dd4 100644 --- a/src/MBPT/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -1,5 +1,5 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & - dBSE,dTDA,evDyn,singlet,triplet,linearize,eta, & + dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0) ! Perform G0W0 calculation @@ -24,6 +24,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nBas integer,intent(in) :: nC @@ -124,13 +125,18 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & !------------------------! call self_energy_exchange_diag(nBas,cHF,PHF,ERI_AO,SigX) - call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) -!--------------------------------! -! Compute renormalization factor ! -!--------------------------------! + if(regularize) then - call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + + else + + call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + + end if !-----------------------------------! ! Solve the quasi-particle equation ! diff --git a/src/MBPT/G0W0_SOSEX.f90 b/src/GW/G0W0_SOSEX.f90 similarity index 100% rename from src/MBPT/G0W0_SOSEX.f90 rename to src/GW/G0W0_SOSEX.f90 diff --git a/src/MBPT/QP_graph.f90 b/src/GW/QP_graph.f90 similarity index 100% rename from src/MBPT/QP_graph.f90 rename to src/GW/QP_graph.f90 diff --git a/src/MBPT/QP_roots.f90 b/src/GW/QP_roots.f90 similarity index 100% rename from src/MBPT/QP_roots.f90 rename to src/GW/QP_roots.f90 diff --git a/src/MBPT/Sangalli_dynamic_perturbation.f90 b/src/GW/Sangalli_dynamic_perturbation.f90 similarity index 100% rename from src/MBPT/Sangalli_dynamic_perturbation.f90 rename to src/GW/Sangalli_dynamic_perturbation.f90 diff --git a/src/MBPT/SigmaC.f90 b/src/GW/SigmaC.f90 similarity index 100% rename from src/MBPT/SigmaC.f90 rename to src/GW/SigmaC.f90 diff --git a/src/MBPT/UG0W0.f90 b/src/GW/UG0W0.f90 similarity index 98% rename from src/MBPT/UG0W0.f90 rename to src/GW/UG0W0.f90 index 90d4ab1..9c7364d 100644 --- a/src/MBPT/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -1,5 +1,5 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & - linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & + linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb, & dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eGW) ! Perform unrestricted G0W0 calculation @@ -24,6 +24,7 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev logical,intent(in) :: spin_flip logical,intent(in) :: linearize double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) diff --git a/src/MBPT/USigmaC.f90 b/src/GW/USigmaC.f90 similarity index 100% rename from src/MBPT/USigmaC.f90 rename to src/GW/USigmaC.f90 diff --git a/src/MBPT/dSigmaC.f90 b/src/GW/dSigmaC.f90 similarity index 100% rename from src/MBPT/dSigmaC.f90 rename to src/GW/dSigmaC.f90 diff --git a/src/MBPT/dUSigmaC.f90 b/src/GW/dUSigmaC.f90 similarity index 100% rename from src/MBPT/dUSigmaC.f90 rename to src/GW/dUSigmaC.f90 diff --git a/src/MBPT/evGW.f90 b/src/GW/evGW.f90 similarity index 87% rename from src/MBPT/evGW.f90 rename to src/GW/evGW.f90 index c25f9ab..4790cf4 100644 --- a/src/MBPT/evGW.f90 +++ b/src/GW/evGW.f90 @@ -1,5 +1,5 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & - G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & + G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -29,6 +29,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta + logical,intent(in) :: regularize + integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -151,7 +153,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,OmRPA, & rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - endif + end if ! Compute spectral weights @@ -161,17 +163,33 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(G0W) then -! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) - call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) - call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + if(regularize) then + + call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + + else + + call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + + end if else -! call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) - call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) - call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + if(regularize) then - endif + call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) + call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + + else + + call self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + + end if + + end if ! Solve the quasi-particle equation @@ -198,9 +216,9 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, call DIIS_extrapolation(rcond,nBas,nBas,n_diis,error_diis,e_diis,eGW-eOld,eGW) else n_diis = 0 - endif + end if - endif + end if ! Save quasiparticles energy for next cycle @@ -210,7 +228,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, nSCF = nSCF + 1 - enddo + end do !------------------------------------------------------------------------ ! End main loop !------------------------------------------------------------------------ @@ -231,7 +249,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, stop - endif + end if ! Deallocate memory @@ -288,6 +306,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, end if - endif + end if end subroutine evGW diff --git a/src/MBPT/evUGW.f90 b/src/GW/evUGW.f90 similarity index 99% rename from src/MBPT/evUGW.f90 rename to src/GW/evUGW.f90 index ece48d8..3e5faec 100644 --- a/src/MBPT/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -1,5 +1,5 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & - G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nBas,nC,nO,nV,nR,nS,ENuc, & + G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc, & EUHF,S,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -29,6 +29,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) diff --git a/src/MBPT/exchange_matrix_MO_basis.f90 b/src/GW/exchange_matrix_MO_basis.f90 similarity index 100% rename from src/MBPT/exchange_matrix_MO_basis.f90 rename to src/GW/exchange_matrix_MO_basis.f90 diff --git a/src/MBPT/excitation_density.f90 b/src/GW/excitation_density.f90 similarity index 100% rename from src/MBPT/excitation_density.f90 rename to src/GW/excitation_density.f90 diff --git a/src/MBPT/excitation_density_SOSEX.f90 b/src/GW/excitation_density_SOSEX.f90 similarity index 100% rename from src/MBPT/excitation_density_SOSEX.f90 rename to src/GW/excitation_density_SOSEX.f90 diff --git a/src/MBPT/plot_GW.f90 b/src/GW/plot_GW.f90 similarity index 100% rename from src/MBPT/plot_GW.f90 rename to src/GW/plot_GW.f90 diff --git a/src/MBPT/print_G0W0.f90 b/src/GW/print_G0W0.f90 similarity index 100% rename from src/MBPT/print_G0W0.f90 rename to src/GW/print_G0W0.f90 diff --git a/src/MBPT/print_SOSEX.f90 b/src/GW/print_SOSEX.f90 similarity index 100% rename from src/MBPT/print_SOSEX.f90 rename to src/GW/print_SOSEX.f90 diff --git a/src/MBPT/print_UG0W0.f90 b/src/GW/print_UG0W0.f90 similarity index 100% rename from src/MBPT/print_UG0W0.f90 rename to src/GW/print_UG0W0.f90 diff --git a/src/MBPT/print_evGW.f90 b/src/GW/print_evGW.f90 similarity index 100% rename from src/MBPT/print_evGW.f90 rename to src/GW/print_evGW.f90 diff --git a/src/MBPT/print_evUGW.f90 b/src/GW/print_evUGW.f90 similarity index 100% rename from src/MBPT/print_evUGW.f90 rename to src/GW/print_evUGW.f90 diff --git a/src/MBPT/print_qsGW.f90 b/src/GW/print_qsGW.f90 similarity index 100% rename from src/MBPT/print_qsGW.f90 rename to src/GW/print_qsGW.f90 diff --git a/src/MBPT/print_qsUGW.f90 b/src/GW/print_qsUGW.f90 similarity index 100% rename from src/MBPT/print_qsUGW.f90 rename to src/GW/print_qsUGW.f90 diff --git a/src/MBPT/qsGW.f90 b/src/GW/qsGW.f90 similarity index 91% rename from src/MBPT/qsGW.f90 rename to src/GW/qsGW.f90 index 2b7b289..9a418aa 100644 --- a/src/MBPT/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -1,5 +1,5 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & - G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & + G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -27,6 +27,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) @@ -192,15 +193,31 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(G0W) then -! call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) - call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) - call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + if(regularize) then + + call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + + else + + call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + + end if else -! call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) - call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) - call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + if(regularize) then + + call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) + call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + + else + + call self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) + call renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + + endif endif diff --git a/src/MBPT/qsGW_PT.f90 b/src/GW/qsGW_PT.f90 similarity index 100% rename from src/MBPT/qsGW_PT.f90 rename to src/GW/qsGW_PT.f90 diff --git a/src/MBPT/qsUGW.f90 b/src/GW/qsUGW.f90 similarity index 99% rename from src/MBPT/qsUGW.f90 rename to src/GW/qsUGW.f90 index b9d055d..22745fe 100644 --- a/src/MBPT/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -1,5 +1,5 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & - G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & + G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_AO,dipole_int_aa, & dipole_int_bb,PHF,cHF,eHF) @@ -28,6 +28,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE logical,intent(in) :: spin_conserved logical,intent(in) :: spin_flip double precision,intent(in) :: eta + logical,intent(in) :: regularize integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) diff --git a/src/GW/regularized_renormalization_factor.f90 b/src/GW/regularized_renormalization_factor.f90 new file mode 100644 index 0000000..98f3084 --- /dev/null +++ b/src/GW/regularized_renormalization_factor.f90 @@ -0,0 +1,87 @@ +subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) + +! Compute the regularized version of the GW renormalization factor + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: COHSEX + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: i,a,p,jb + double precision :: eps + + double precision :: kappa + double precision :: fk,dfk + +! Output variables + + double precision,intent(out) :: Z(nBas) + +! Initialize + + Z(:) = 0d0 + +!-----------------------------------------! +! Parameters for regularized calculations ! +!-----------------------------------------! + + kappa = 1.1d0 + +! static COHSEX approximation + + if(COHSEX) then + + Z(:) = 1d0 + return + + else + + ! Occupied part of the correlation self-energy + + do p=nC+1,nBas-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 + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*dfk + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC+1,nBas-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 + dfk = - 1d0/eps + 2d0*kappa*exp(-kappa*abs(eps))/(1d0 - exp(-kappa*abs(eps))) + dfk = dfk*fk + Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2 + end do + end do + end do + + end if + +! Compute renormalization factor from derivative of SigC + + Z(:) = 1d0/(1d0 - Z(:)) + +end subroutine regularized_renormalization_factor diff --git a/src/MBPT/regularized_self_energy_correlation.f90 b/src/GW/regularized_self_energy_correlation.f90 similarity index 92% rename from src/MBPT/regularized_self_energy_correlation.f90 rename to src/GW/regularized_self_energy_correlation.f90 index 1a56dcd..dcb6ffc 100644 --- a/src/MBPT/regularized_self_energy_correlation.f90 +++ b/src/GW/regularized_self_energy_correlation.f90 @@ -9,7 +9,12 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e, logical,intent(in) :: COHSEX double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) diff --git a/src/MBPT/regularized_self_energy_correlation_diag.f90 b/src/GW/regularized_self_energy_correlation_diag.f90 similarity index 93% rename from src/MBPT/regularized_self_energy_correlation_diag.f90 rename to src/GW/regularized_self_energy_correlation_diag.f90 index ec7bb5f..2856e84 100644 --- a/src/MBPT/regularized_self_energy_correlation_diag.f90 +++ b/src/GW/regularized_self_energy_correlation_diag.f90 @@ -23,7 +23,6 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, integer :: i,a,p,q,jb double precision :: eps - double precision,external :: SigC_dcgw double precision :: kappa double precision :: fk @@ -37,9 +36,9 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, SigC(:) = 0d0 -!---------------------------------------------! -! Parameters for regularized MP2 calculations ! -!---------------------------------------------! +!-----------------------------------------! +! Parameters for regularized calculations ! +!-----------------------------------------! kappa = 1.1d0 diff --git a/src/MBPT/renormalization_factor.f90 b/src/GW/renormalization_factor.f90 similarity index 63% rename from src/MBPT/renormalization_factor.f90 rename to src/GW/renormalization_factor.f90 index 10f2849..eb699bf 100644 --- a/src/MBPT/renormalization_factor.f90 +++ b/src/GW/renormalization_factor.f90 @@ -9,14 +9,19 @@ subroutine renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) logical,intent(in) :: COHSEX double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) ! Local variables - integer :: i,j,a,b,x,jb + integer :: p,i,a,jb double precision :: eps ! Output variables @@ -38,30 +43,22 @@ subroutine renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) ! Occupied part of the correlation self-energy - do x=nC+1,nBas-nR + do p=nC+1,nBas-nR do i=nC+1,nO - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = e(x) - e(i) + Omega(jb) - Z(x) = Z(x) - 2d0*rho(x,i,jb)**2*(eps/(eps**2 + eta**2))**2 - end do + do jb=1,nS + eps = e(p) - e(i) + Omega(jb) + Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2 end do end do end do ! Virtual part of the correlation self-energy - do x=nC+1,nBas-nR + do p=nC+1,nBas-nR do a=nO+1,nBas-nR - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = e(x) - e(a) - Omega(jb) - Z(x) = Z(x) - 2d0*rho(x,a,jb)**2*(eps/(eps**2 + eta**2))**2 - end do + do jb=1,nS + eps = e(p) - e(a) - Omega(jb) + Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2 end do end do end do diff --git a/src/MBPT/renormalization_factor_SOSEX.f90 b/src/GW/renormalization_factor_SOSEX.f90 similarity index 100% rename from src/MBPT/renormalization_factor_SOSEX.f90 rename to src/GW/renormalization_factor_SOSEX.f90 diff --git a/src/MBPT/self_energy_correlation.f90 b/src/GW/self_energy_correlation.f90 similarity index 91% rename from src/MBPT/self_energy_correlation.f90 rename to src/GW/self_energy_correlation.f90 index 369ee32..4d62b43 100644 --- a/src/MBPT/self_energy_correlation.f90 +++ b/src/GW/self_energy_correlation.f90 @@ -9,7 +9,12 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec logical,intent(in) :: COHSEX double precision,intent(in) :: eta - integer,intent(in) :: nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) diff --git a/src/MBPT/self_energy_correlation_SOSEX_diag.f90 b/src/GW/self_energy_correlation_SOSEX_diag.f90 similarity index 100% rename from src/MBPT/self_energy_correlation_SOSEX_diag.f90 rename to src/GW/self_energy_correlation_SOSEX_diag.f90 diff --git a/src/MBPT/self_energy_correlation_diag.f90 b/src/GW/self_energy_correlation_diag.f90 similarity index 98% rename from src/MBPT/self_energy_correlation_diag.f90 rename to src/GW/self_energy_correlation_diag.f90 index 59aab3d..9adc747 100644 --- a/src/MBPT/self_energy_correlation_diag.f90 +++ b/src/GW/self_energy_correlation_diag.f90 @@ -23,7 +23,6 @@ subroutine self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,r integer :: i,a,p,q,jb double precision :: eps - double precision,external :: SigC_dcgw ! Output variables diff --git a/src/MBPT/self_energy_exchange.f90 b/src/GW/self_energy_exchange.f90 similarity index 100% rename from src/MBPT/self_energy_exchange.f90 rename to src/GW/self_energy_exchange.f90 diff --git a/src/MBPT/self_energy_exchange_diag.f90 b/src/GW/self_energy_exchange_diag.f90 similarity index 100% rename from src/MBPT/self_energy_exchange_diag.f90 rename to src/GW/self_energy_exchange_diag.f90 diff --git a/src/MBPT/static_screening_WA.f90 b/src/GW/static_screening_WA.f90 similarity index 100% rename from src/MBPT/static_screening_WA.f90 rename to src/GW/static_screening_WA.f90 diff --git a/src/MBPT/static_screening_WB.f90 b/src/GW/static_screening_WB.f90 similarity index 100% rename from src/MBPT/static_screening_WB.f90 rename to src/GW/static_screening_WB.f90 diff --git a/src/MBPT/ufBSE.f90 b/src/GW/ufBSE.f90 similarity index 97% rename from src/MBPT/ufBSE.f90 rename to src/GW/ufBSE.f90 index 9ee38c6..426b4bc 100644 --- a/src/MBPT/ufBSE.f90 +++ b/src/GW/ufBSE.f90 @@ -1,4 +1,4 @@ -subroutine ufBSE(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) +subroutine ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Unfold BSE@GW equations @@ -7,7 +7,6 @@ subroutine ufBSE(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eGW) ! Input variables - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO diff --git a/src/MBPT/ufG0W0.f90 b/src/GW/ufG0W0.f90 similarity index 97% rename from src/MBPT/ufG0W0.f90 rename to src/GW/ufG0W0.f90 index 8211821..51ac629 100644 --- a/src/MBPT/ufG0W0.f90 +++ b/src/GW/ufG0W0.f90 @@ -1,4 +1,4 @@ -subroutine ufG0W0(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +subroutine ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Unfold G0W0 equations @@ -7,7 +7,6 @@ subroutine ufG0W0(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Input variables - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO diff --git a/src/MBPT/ufGW.f90 b/src/GW/ufGW.f90 similarity index 97% rename from src/MBPT/ufGW.f90 rename to src/GW/ufGW.f90 index f55d571..d64807e 100644 --- a/src/MBPT/ufGW.f90 +++ b/src/GW/ufGW.f90 @@ -1,4 +1,4 @@ -subroutine ufGW(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) +subroutine ufGW(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Unfold GW equations @@ -7,7 +7,6 @@ subroutine ufGW(eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Input variables - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO diff --git a/src/MBPT/unrestricted_Bethe_Salpeter.f90 b/src/GW/unrestricted_Bethe_Salpeter.f90 similarity index 100% rename from src/MBPT/unrestricted_Bethe_Salpeter.f90 rename to src/GW/unrestricted_Bethe_Salpeter.f90 diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix.f90 b/src/GW/unrestricted_Bethe_Salpeter_A_matrix.f90 similarity index 100% rename from src/MBPT/unrestricted_Bethe_Salpeter_A_matrix.f90 rename to src/GW/unrestricted_Bethe_Salpeter_A_matrix.f90 diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 b/src/GW/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 similarity index 100% rename from src/MBPT/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 rename to src/GW/unrestricted_Bethe_Salpeter_A_matrix_dynamic.f90 diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_B_matrix.f90 b/src/GW/unrestricted_Bethe_Salpeter_B_matrix.f90 similarity index 100% rename from src/MBPT/unrestricted_Bethe_Salpeter_B_matrix.f90 rename to src/GW/unrestricted_Bethe_Salpeter_B_matrix.f90 diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 b/src/GW/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 similarity index 100% rename from src/MBPT/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 rename to src/GW/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 b/src/GW/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 similarity index 100% rename from src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 rename to src/GW/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 diff --git a/src/MBPT/unrestricted_QP_graph.f90 b/src/GW/unrestricted_QP_graph.f90 similarity index 100% rename from src/MBPT/unrestricted_QP_graph.f90 rename to src/GW/unrestricted_QP_graph.f90 diff --git a/src/MBPT/unrestricted_excitation_density.f90 b/src/GW/unrestricted_excitation_density.f90 similarity index 100% rename from src/MBPT/unrestricted_excitation_density.f90 rename to src/GW/unrestricted_excitation_density.f90 diff --git a/src/GW/unrestricted_regularized_renormalization_factor.f90 b/src/GW/unrestricted_regularized_renormalization_factor.f90 new file mode 100644 index 0000000..34a8a90 --- /dev/null +++ b/src/GW/unrestricted_regularized_renormalization_factor.f90 @@ -0,0 +1,111 @@ +subroutine unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,Z) + +! Compute the renormalization factor in the unrestricted formalism + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nSt + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: Omega(nSt) + double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) + +! Local variables + + integer :: i,a,p,jb + double precision :: eps + + double precision :: kappa + double precision :: fk,dfk + +! Output variables + + double precision,intent(out) :: Z(nBas,nspin) + +! Initialize + + Z(:,:) = 0d0 + +!-----------------------------------------! +! Parameters for regularized calculations ! +!-----------------------------------------! + + kappa = 1.1d0 + +!--------------! +! Spin-up part ! +!--------------! + + ! Occupied part of the renormalization factor + + do p=nC(1)+1,nBas-nR(1) + do i=nC(1)+1,nO(1) + do jb=1,nSt + eps = e(p,1) - e(i,1) + 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 + Z(p,1) = Z(p,1) + rho(p,i,jb,1)**2*dfk + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do a=nO(1)+1,nBas-nR(1) + do jb=1,nSt + eps = e(p,1) - e(a,1) - 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 + Z(p,1) = Z(p,1) + rho(p,a,jb,1)**2*dfk + end do + end do + end do + +!----------------! +! Spin-down part ! +!----------------! + + ! Occupied part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do i=nC(2)+1,nO(2) + do jb=1,nSt + eps = e(p,2) - e(i,2) + 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 + Z(p,2) = Z(p,2) + rho(p,i,jb,2)**2*dfk + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do a=nO(2)+1,nBas-nR(2) + do jb=1,nSt + eps = e(p,2) - e(a,2) - 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 + Z(p,2) = Z(p,2) + rho(p,a,jb,2)**2*dfk + end do + end do + end do + +! Final rescaling + + Z(:,:) = 1d0/(1d0 + Z(:,:)) + +end subroutine unrestricted_regularized_renormalization_factor diff --git a/src/GW/unrestricted_regularized_self_energy_correlation.f90 b/src/GW/unrestricted_regularized_self_energy_correlation.f90 new file mode 100644 index 0000000..9553e84 --- /dev/null +++ b/src/GW/unrestricted_regularized_self_energy_correlation.f90 @@ -0,0 +1,133 @@ +subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nSt + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: Omega(nSt) + double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) + +! Local variables + + integer :: i,a,p,q,jb + double precision :: eps + + double precision :: kappa + double precision :: fk + +! Output variables + + double precision,intent(out) :: SigC(nBas,nBas,nspin) + double precision :: EcGM(nspin) + +! Initialize + + SigC(:,:,:) = 0d0 + EcGM(:) = 0d0 + +!-----------------------------------------! +! Parameters for regularized calculations ! +!-----------------------------------------! + + kappa = 1.1d0 + +!--------------! +! Spin-up part ! +!--------------! + + ! Occupied part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do q=nC(1)+1,nBas-nR(1) + do i=nC(1)+1,nO(1) + do jb=1,nSt + eps = e(p,1) - e(i,1) + Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + SigC(p,q,1) = SigC(p,q,1) + rho(p,i,jb,1)*rho(q,i,jb,1)*eps/(eps**2 + eta**2) + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do q=nC(1)+1,nBas-nR(1) + do a=nO(1)+1,nBas-nR(1) + do jb=1,nSt + eps = e(p,1) - e(a,1) - Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + SigC(p,q,1) = SigC(p,q,1) + rho(p,a,jb,1)*rho(q,a,jb,1)*eps/(eps**2 + eta**2) + end do + end do + end do + end do + + ! GM correlation energy + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do jb=1,nSt + eps = e(a,1) - e(i,1) + Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*eps/(eps**2 + eta**2) + end do + end do + end do + +!----------------! +! Spin-down part ! +!----------------! + + ! Occupied part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do q=nC(2)+1,nBas-nR(2) + do i=nC(2)+1,nO(2) + do jb=1,nSt + eps = e(p,2) - e(i,2) + Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + SigC(p,q,2) = SigC(p,q,2) + rho(p,i,jb,2)*rho(q,i,jb,2)*eps/(eps**2 + eta**2) + end do + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do q=nC(2)+1,nBas-nR(2) + do a=nO(2)+1,nBas-nR(2) + do jb=1,nSt + eps = e(p,2) - e(a,2) - Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + SigC(p,q,2) = SigC(p,q,2) + rho(p,a,jb,2)*rho(q,a,jb,2)*eps/(eps**2 + eta**2) + end do + end do + end do + end do + + ! GM correlation energy + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do jb=1,nSt + eps = e(a,2) - e(i,2) + Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*eps/(eps**2 + eta**2) + end do + end do + end do + +end subroutine unrestricted_self_energy_correlation diff --git a/src/GW/unrestricted_regularized_self_energy_correlation_diag.f90 b/src/GW/unrestricted_regularized_self_energy_correlation_diag.f90 new file mode 100644 index 0000000..02de065 --- /dev/null +++ b/src/GW/unrestricted_regularized_self_energy_correlation_diag.f90 @@ -0,0 +1,126 @@ +subroutine unrestricted_regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM) + +! Compute diagonal of the correlation part of the self-energy + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nSt + double precision,intent(in) :: e(nBas,nspin) + double precision,intent(in) :: Omega(nSt) + double precision,intent(in) :: rho(nBas,nBas,nSt,nspin) + +! Local variables + + integer :: i,a,p,q,jb + double precision :: eps + + double precision :: kappa + double precision :: fk + +! Output variables + + double precision,intent(out) :: SigC(nBas,nspin) + double precision :: EcGM(nspin) + +! Initialize + + SigC(:,:) = 0d0 + EcGM(:) = 0d0 + +!-----------------------------------------! +! Parameters for regularized calculations ! +!-----------------------------------------! + + kappa = 1.1d0 + +!--------------! +! Spin-up part ! +!--------------! + + ! Occupied part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do i=nC(1)+1,nO(1) + do jb=1,nSt + eps = e(p,1) - e(i,1) + Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + SigC(p,1) = SigC(p,1) + rho(p,i,jb,1)**2*fk + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(1)+1,nBas-nR(1) + do a=nO(1)+1,nBas-nR(1) + do jb=1,nSt + eps = e(p,1) - e(a,1) - Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + SigC(p,1) = SigC(p,1) + rho(p,a,jb,1)**2*fk + end do + end do + end do + + ! GM correlation energy + + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + do jb=1,nSt + eps = e(a,1) - e(i,1) + Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + EcGM(1) = EcGM(1) - rho(a,i,jb,1)**2*fk + end do + end do + end do + +!----------------! +! Spin-down part ! +!----------------! + + ! Occupied part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do i=nC(2)+1,nO(2) + do jb=1,nSt + eps = e(p,2) - e(i,2) + Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + SigC(p,2) = SigC(p,2) + rho(p,i,jb,2)**2*fk + end do + end do + end do + + ! Virtual part of the correlation self-energy + + do p=nC(2)+1,nBas-nR(2) + do a=nO(2)+1,nBas-nR(2) + do jb=1,nSt + eps = e(p,2) - e(a,2) - Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + SigC(p,2) = SigC(p,2) + rho(p,a,jb,2)**2*fk + end do + end do + end do + + ! GM correlation energy + + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + do jb=1,nSt + eps = e(a,2) - e(i,2) + Omega(jb) + fk = (1d0 - exp(-kappa*abs(eps)))**2/eps + EcGM(2) = EcGM(2) - rho(a,i,jb,2)**2*fk + end do + end do + end do + + +end subroutine unrestricted_regularized_self_energy_correlation_diag diff --git a/src/MBPT/unrestricted_renormalization_factor.f90 b/src/GW/unrestricted_renormalization_factor.f90 similarity index 97% rename from src/MBPT/unrestricted_renormalization_factor.f90 rename to src/GW/unrestricted_renormalization_factor.f90 index f175d93..147b483 100644 --- a/src/MBPT/unrestricted_renormalization_factor.f90 +++ b/src/GW/unrestricted_renormalization_factor.f90 @@ -20,7 +20,7 @@ subroutine unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nSt,e,Omega, ! Local variables - integer :: i,j,a,b,p,q,jb + integer :: i,a,p,jb double precision :: eps ! Output variables diff --git a/src/MBPT/unrestricted_self_energy_correlation.f90 b/src/GW/unrestricted_self_energy_correlation.f90 similarity index 98% rename from src/MBPT/unrestricted_self_energy_correlation.f90 rename to src/GW/unrestricted_self_energy_correlation.f90 index d6c5f98..de1ae28 100644 --- a/src/MBPT/unrestricted_self_energy_correlation.f90 +++ b/src/GW/unrestricted_self_energy_correlation.f90 @@ -20,7 +20,7 @@ subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega ! Local variables - integer :: i,j,a,b,p,q,jb + integer :: i,a,p,q,jb double precision :: eps ! Output variables diff --git a/src/MBPT/unrestricted_self_energy_correlation_diag.f90 b/src/GW/unrestricted_self_energy_correlation_diag.f90 similarity index 98% rename from src/MBPT/unrestricted_self_energy_correlation_diag.f90 rename to src/GW/unrestricted_self_energy_correlation_diag.f90 index d683380..8290d9d 100644 --- a/src/MBPT/unrestricted_self_energy_correlation_diag.f90 +++ b/src/GW/unrestricted_self_energy_correlation_diag.f90 @@ -20,7 +20,7 @@ subroutine unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nSt,e, ! Local variables - integer :: i,j,a,b,p,q,jb + integer :: i,a,p,q,jb double precision :: eps ! Output variables diff --git a/src/MBPT/Makefile b/src/MBPT/Makefile deleted file mode 100644 index 6ca514e..0000000 --- a/src/MBPT/Makefile +++ /dev/null @@ -1,10 +0,0 @@ -default: - ninja - make -C .. - -clean: - ninja -t clean - -debug: - ninja -t clean - make -C .. debug diff --git a/src/MBPT/Sangalli_A_matrix_dynamic.f90.x b/src/MBPT/Sangalli_A_matrix_dynamic.f90.x deleted file mode 100644 index 75f8ffd..0000000 --- a/src/MBPT/Sangalli_A_matrix_dynamic.f90.x +++ /dev/null @@ -1,78 +0,0 @@ -subroutine Sangalli_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,OmRPA,OmBSE,rho,A_dyn) - -! Compute the dynamic part of the Bethe-Salpeter equation matrices - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas,nC,nO,nV,nR,nS - double precision,intent(in) :: eta - double precision,intent(in) :: lambda - double precision,intent(in) :: eGW(nBas) - double precision,intent(in) :: OmRPA(nS) - double precision,intent(in) :: OmBSE - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: maxS - double precision :: chi - double precision :: eps - integer :: i,j,a,b,ia,jb,kc - -! Output variables - - double precision,intent(out) :: A_dyn(nS,nS) - -! Initialization - - A_dyn(:,:) = 0d0 - -! Number of poles taken into account - - maxS = nS - -! Build dynamic A matrix - - do - - - (ERI(i,k,j,c)*KroneckerDelta(a,b) + ERI(k,a,c,b)*KroneckerDelta(i,j))*(X) - - ia = 0 - do i=nC+1,nO - do a=nO+1,nBas-nR - ia = ia + 1 - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - - chi = 0d0 - do kc=1,maxS - - chi = chi + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/(OmRPA(kc)**2 + eta**2) - - enddo - - A_dyn(ia,jb) = A_dyn(ia,jb) - 4d0*lambda*chi - - chi = 0d0 - do kc=1,maxS - do ld=1,maxS - - eps = OmBSE - (OmRPA(kc) + OmRPA(ld)) - chi = chi + cRPA(ia,kc,ld)*cRPA(jb,kc,ld)*eps/(eps**2 + (2d0*eta)**2) - - enddo - - A_dyn(ia,jb) = A_dyn(ia,jb) - 2d0*lambda*chi - - enddo - enddo - enddo - enddo - -end subroutine Sangalli_A_matrix_dynamic diff --git a/src/MBPT/excitation_density_RI.f90 b/src/MBPT/excitation_density_RI.f90 deleted file mode 100644 index 967f0bb..0000000 --- a/src/MBPT/excitation_density_RI.f90 +++ /dev/null @@ -1,65 +0,0 @@ -subroutine excitation_density_RI(nBas,nC,nO,nR,nS,c,G,XpY,rho) - -! Compute excitation densities in the RI approximation - - implicit none - -! Input variables - - integer,intent(in) :: nBas,nC,nO,nR,nS - double precision,intent(in) :: c(nBas,nBas),G(nBas,nBas,nBas,nBas),XpY(nS,nS) - -! Local variables - - double precision,allocatable :: scr(:,:,:) - integer :: mu,nu,la,si,ia,jb,x,y,j,b - -! Output variables - - double precision,intent(out) :: rho(nBas,nBas,nS) - -! Memory allocation - allocate(scr(nBas,nBas,nS)) - - rho(:,:,:) = 0d0 - do nu=1,nBas - do si=1,nBas - do ia=1,nS - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - rho(nu,si,ia) = rho(nu,si,ia) + c(nu,j)*XpY(ia,jb)*c(si,b) - enddo - enddo - enddo - enddo - enddo - - scr(:,:,:) = 0d0 - do mu=1,nBas - do la=1,nBas - do ia=1,nS - do nu=1,nBas - do si=1,nBas - scr(mu,la,ia) = scr(mu,la,ia) + G(mu,nu,la,si)*rho(nu,si,ia) - enddo - enddo - enddo - enddo - enddo - - rho(:,:,:) = 0d0 - do ia=1,nS - do x=nC+1,nBas-nR - do y=nC+1,nBas-nR - do mu=1,nBas - do la=1,nBas - rho(x,y,ia) = rho(x,y,ia) + c(mu,x)*scr(mu,la,ia)*c(la,y) - enddo - enddo - enddo - enddo - enddo - -end subroutine excitation_density_RI diff --git a/src/MBPT/excitation_density_SOSEX_RI.f90 b/src/MBPT/excitation_density_SOSEX_RI.f90 deleted file mode 100644 index 8f3f6b5..0000000 --- a/src/MBPT/excitation_density_SOSEX_RI.f90 +++ /dev/null @@ -1,65 +0,0 @@ -subroutine excitation_density_SOSEX_RI(nBas,nC,nO,nR,nS,c,G,XpY,rho) - -! Compute excitation densities - - implicit none - -! Input variables - - integer,intent(in) :: nBas,nC,nO,nR,nS - double precision,intent(in) :: c(nBas,nBas),G(nBas,nBas,nBas,nBas),XpY(nS,nS) - -! Local variables - - double precision,allocatable :: scr(:,:,:) - integer :: mu,nu,la,si,ia,jb,x,y,j,b - -! Output variables - - double precision,intent(out) :: rho(nBas,nBas,nS) - -! Memory allocation - allocate(scr(nBas,nBas,nS)) - - rho(:,:,:) = 0d0 - do nu=1,nBas - do si=1,nBas - do ia=1,nS - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - rho(nu,si,ia) = rho(nu,si,ia) + c(nu,j)*XpY(ia,jb)*c(si,b) - enddo - enddo - enddo - enddo - enddo - - scr(:,:,:) = 0d0 - do mu=1,nBas - do la=1,nBas - do ia=1,nS - do nu=1,nBas - do si=1,nBas - scr(mu,la,ia) = scr(mu,la,ia) + G(mu,nu,la,si)*rho(nu,si,ia) - enddo - enddo - enddo - enddo - enddo - - rho(:,:,:) = 0d0 - do ia=1,nS - do x=nC+1,nBas-nR - do y=nC+1,nBas-nR - do mu=1,nBas - do la=1,nBas - rho(x,y,ia) = rho(x,y,ia) + c(mu,x)*scr(mu,la,ia)*c(la,y) - enddo - enddo - enddo - enddo - enddo - -end subroutine excitation_density_SOSEX_RI diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index c6b0444..4adadbd 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -124,14 +124,19 @@ program QuAcK integer :: maxSCF_GF,n_diis_GF,renormGF double precision :: thresh_GF - logical :: DIIS_GF,linGF + logical :: DIIS_GF,linGF,regGF double precision :: eta_GF integer :: maxSCF_GW,n_diis_GW double precision :: thresh_GW - logical :: DIIS_GW,COHSEX,SOSEX,TDA_W,G0W,GW0,linGW + logical :: DIIS_GW,COHSEX,SOSEX,TDA_W,G0W,GW0,linGW,regGW double precision :: eta_GW + integer :: maxSCF_GT,n_diis_GT + double precision :: thresh_GT + logical :: DIIS_GT,TDA_T,linGT,regGT + double precision :: eta_GT + logical :: BSE,dBSE,dTDA,evDyn integer :: nMC,nEq,nWalk,nPrint,iSeed @@ -176,9 +181,10 @@ program QuAcK call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, & 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, & - maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, & + maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, & COHSEX,SOSEX,TDA_W,G0W,GW0, & + maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & doACFDT,exchange_kernel,doXBS, & BSE,dBSE,dTDA,evDyn, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) @@ -868,13 +874,14 @@ program QuAcK if(unrestricted) then - call UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, & - S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb,eHF) + call UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linGF,eta_GF,regGF, & + nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & + dipole_int_aa,dipole_int_bb,eHF) else - call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, & - eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF,eta_GF,regGF, & + nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) end if @@ -897,13 +904,13 @@ program QuAcK if(unrestricted) then call evUGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & - eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & + eta_GF,regGF,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & dipole_int_aa,dipole_int_bb,cHF,eHF) else call evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF_GF,thresh_GF,n_diis_GF, & - singlet,triplet,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & + singlet,triplet,linGF,eta_GF,regGF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & ERI_MO,dipole_int_MO,eHF) end if @@ -926,13 +933,13 @@ program QuAcK if(unrestricted) then - call qsUGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GF, & + call qsUGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GF,regGF, & nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO, & ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_AO,dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) else - call qsGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GF,nNuc,ZNuc,rNuc,ENuc, & + call qsGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GF,regGF,nNuc,ZNuc,rNuc,ENuc, & nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) end if @@ -989,7 +996,7 @@ program QuAcK if(unrestricted) then call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & - linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & + linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb, & dipole_int_aa,dipole_int_bb,PHF,cHF,eHF,Vxc,eG0W0) else @@ -997,10 +1004,10 @@ program QuAcK if(SOSEX) then call G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & - eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) + eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) else call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & - linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) + linGW,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) end if end if @@ -1023,14 +1030,14 @@ program QuAcK if(unrestricted) then call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & - G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc, & + G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,regGW,nBas,nC,nO,nV,nR,nS,ENuc, & EUHF,S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb, & PHF,cHF,eHF,Vxc,eG0W0) else call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX, & - BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & + BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,regGW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0W0) end if call cpu_time(end_evGW) @@ -1052,14 +1059,14 @@ program QuAcK if(unrestricted) then call qsUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & - G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & + G0W,GW0,dBSE,dTDA,evDyn,spin_conserved,spin_flip,eta_GW,regGW,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO, & nV,nR,nS,EUHF,S,X,T,V,Hc,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_AO, & dipole_int_aa,dipole_int_bb,PHF,cHF,eHF) else call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX, & - BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,nNuc,ZNuc,rNuc,ENuc, & + BSE,TDA_W,TDA,G0W,GW0,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,regGW,nNuc,ZNuc,rNuc,ENuc, & nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) end if @@ -1079,7 +1086,7 @@ program QuAcK if(doufG0W0) then call cpu_time(start_ufGW) - call ufG0W0(eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_ufGW) t_ufGW = end_ufGW - start_ufGW @@ -1095,14 +1102,14 @@ program QuAcK if(doufGW) then call cpu_time(start_ufGW) - call ufGW(eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufGW(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_ufGW) t_ufGW = end_ufGW - start_ufGW write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_ufGW,' seconds' write(*,*) -! call ufBSE(eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) +! call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) end if @@ -1115,8 +1122,8 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) - call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet, & - linGW,eta_GW,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO, & + call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet, & + linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO, & PHF,cHF,eHF,Vxc,eG0T0) call cpu_time(end_G0T0) @@ -1133,8 +1140,8 @@ program QuAcK if(doevGT) then call cpu_time(start_evGT) - call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & - BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW, & + call evGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & + BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO, & PHF,cHF,eHF,Vxc,eG0T0) call cpu_time(end_evGT) @@ -1153,9 +1160,10 @@ program QuAcK call cpu_time(start_qsGT) - call qsGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & - BSE,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GW,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + call qsGT(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & + BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT, & + nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & + ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) call cpu_time(end_qsGT) diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 2d7aa4f..382ec95 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -1,9 +1,10 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, & 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, & - maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW, & + maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, & COHSEX,SOSEX,TDA_W,G0W,GW0, & + maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & doACFDT,exchange_kernel,doXBS, & BSE,dBSE,dTDA,evDyn, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) @@ -41,6 +42,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: linGF integer,intent(out) :: renormGF double precision,intent(out) :: eta_GF + logical,intent(out) :: regGF integer,intent(out) :: maxSCF_GW double precision,intent(out) :: thresh_GW @@ -53,6 +55,16 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: GW0 logical,intent(out) :: linGW double precision,intent(out) :: eta_GW + logical,intent(out) :: regGW + + integer,intent(out) :: maxSCF_GT + double precision,intent(out) :: thresh_GT + logical,intent(out) :: DIIS_GT + integer,intent(out) :: n_diis_GT + logical,intent(out) :: TDA_T + logical,intent(out) :: linGT + double precision,intent(out) :: eta_GT + logical,intent(out) :: regGT logical,intent(out) :: doACFDT logical,intent(out) :: exchange_kernel @@ -73,7 +85,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t ! Local variables - character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6,answer7 + character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6,answer7,answer8 ! Open file with method specification @@ -135,7 +147,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t if(answer4 == 'T') spin_conserved = .true. if(answer5 == 'T') spin_flip = .true. -! Read Green function options +! Read GF options maxSCF_GF = 64 thresh_GF = 1d-5 @@ -144,12 +156,14 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t linGF = .false. eta_GF = 0d0 renormGF = 0 + regGF = .false. read(1,*) - read(1,*) maxSCF_GF,thresh_GF,answer1,n_diis_GF,answer2,eta_GF,renormGF + read(1,*) maxSCF_GF,thresh_GF,answer1,n_diis_GF,answer2,eta_GF,renormGF,answer3 if(answer1 == 'T') DIIS_GF = .true. if(answer2 == 'T') linGF = .true. + if(answer3 == 'T') regGF = .true. if(.not.DIIS_GF) n_diis_GF = 1 ! Read GW options @@ -160,6 +174,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t n_diis_GW = 5 linGW = .false. eta_GW = 0d0 + regGW = .false. COHSEX = .false. SOSEX = .false. TDA_W = .false. @@ -168,7 +183,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t read(1,*) read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2,eta_GW, & - answer3,answer4,answer5,answer6,answer7 + answer3,answer4,answer5,answer6,answer7,answer8 if(answer1 == 'T') DIIS_GW = .true. if(answer2 == 'T') linGW = .true. @@ -177,8 +192,30 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t if(answer5 == 'T') TDA_W = .true. if(answer6 == 'T') G0W = .true. if(answer7 == 'T') GW0 = .true. + if(answer8 == 'T') regGW = .true. if(.not.DIIS_GW) n_diis_GW = 1 +! Read GF options + + maxSCF_GF = 64 + thresh_GF = 1d-5 + DIIS_GF = .false. + n_diis_GF = 5 + linGF = .false. + eta_GF = 0d0 + regGF = .false. + TDA_T = .false. + + read(1,*) + read(1,*) maxSCF_GT,thresh_GT,answer1,n_diis_GT,answer2,eta_GT, & + answer3,answer4 + + if(answer1 == 'T') DIIS_GT = .true. + if(answer2 == 'T') linGT = .true. + if(answer3 == 'T') TDA_T = .true. + if(answer4 == 'T') regGT = .true. + if(.not.DIIS_GT) n_diis_GT = 1 + ! Options for adiabatic connection doACFDT = .false.