From 545b5fb5e2013bf45dfc641d0dcf00f78a1319ac Mon Sep 17 00:00:00 2001 From: Antoine Marie Date: Fri, 30 Jun 2023 16:47:26 +0200 Subject: [PATCH] trying to fix diverging version --- basis/aug-cc-pvtz | 177 ++++++++++++++++++ input/methods | 2 +- input/options | 2 +- int/ERI.Hu.dat | 2 - int/Kin.Hu.dat | 2 - int/Nuc.Hu.dat | 2 - int/Ov.Hu.dat | 2 - int/x.Hu.dat | 0 int/y.Hu.dat | 0 int/z.Hu.dat | 0 mol/h2.xyz | 4 - src/GF/G0F2.f90 | 15 +- src/GT/G0T0.f90 | 18 +- src/GW/G0W0.f90 | 6 +- src/GW/QP_graph.f90 | 7 +- src/GW/QP_roots.f90 | 2 +- src/GW/SRG_qsGW.f90 | 2 +- src/GW/SigmaC.f90 | 66 +++++-- src/GW/dSigmaC.f90 | 66 ++++--- src/GW/evGW.f90 | 33 +++- ...gularized_self_energy_correlation_diag.f90 | 25 +-- src/GW/renormalization_factor_SRG.f90 | 28 ++- src/HF/MOM.f90 | 2 +- src/HF/MOM_overlap.f90 | 5 +- src/LR/print_excitation.f90 | 2 +- src/QuAcK/QuAcK.f90 | 6 +- src/RPA/sort_ppRPA.f90 | 14 +- src/make_ninja.py | 14 +- src/utils/read_basis.f90 | 172 ++++++++--------- 29 files changed, 485 insertions(+), 191 deletions(-) delete mode 100644 int/ERI.Hu.dat delete mode 100644 int/Kin.Hu.dat delete mode 100644 int/Nuc.Hu.dat delete mode 100644 int/Ov.Hu.dat delete mode 100644 int/x.Hu.dat delete mode 100644 int/y.Hu.dat delete mode 100644 int/z.Hu.dat delete mode 100644 mol/h2.xyz diff --git a/basis/aug-cc-pvtz b/basis/aug-cc-pvtz index 4e085c0..76ba2b2 100644 --- a/basis/aug-cc-pvtz +++ b/basis/aug-cc-pvtz @@ -1224,6 +1224,183 @@ F 1 F 1 1 0.4060000 1.0000000 +POTASSIUM +S 23 +1 2.709796E+07 1.300000E-05 +2 3.592856E+06 4.100000E-05 +3 7.222638E+05 1.210000E-04 +4 1.833976E+05 3.350000E-04 +5 5.464812E+04 9.160000E-04 +6 1.830913E+04 2.488000E-03 +7 6.712501E+03 6.698000E-03 +8 2.643719E+03 1.765000E-02 +9 1.103835E+03 4.442700E-02 +10 4.836777E+02 1.025880E-01 +11 2.206245E+02 2.038040E-01 +12 1.041117E+02 3.138510E-01 +13 5.052749E+01 3.074340E-01 +14 2.485247E+01 1.398430E-01 +15 1.128102E+01 1.624100E-02 +16 5.510166E+00 -2.630000E-04 +17 2.706914E+00 5.030000E-04 +18 1.122807E+00 -2.220000E-04 +19 5.236385E-01 6.500000E-05 +20 2.355180E-01 -3.600000E-05 +21 4.157977E-02 1.900000E-05 +22 2.555456E-02 -2.100000E-05 +23 1.471495E-02 7.000000E-06 +S 23 +1 2.709796E+07 -4.000000E-06 +2 3.592856E+06 -1.200000E-05 +3 7.222638E+05 -3.500000E-05 +4 1.833976E+05 -9.700000E-05 +5 5.464812E+04 -2.660000E-04 +6 1.830913E+04 -7.240000E-04 +7 6.712501E+03 -1.961000E-03 +8 2.643719E+03 -5.206000E-03 +9 1.103835E+03 -1.339600E-02 +10 4.836777E+02 -3.218800E-02 +11 2.206245E+02 -7.000500E-02 +12 1.041117E+02 -1.277050E-01 +13 5.052749E+01 -1.733640E-01 +14 2.485247E+01 -8.573700E-02 +15 1.128102E+01 2.571140E-01 +16 5.510166E+00 5.577460E-01 +17 2.706914E+00 3.224970E-01 +18 1.122807E+00 3.251000E-02 +19 5.236385E-01 -1.916000E-03 +20 2.355180E-01 1.300000E-03 +21 4.157977E-02 -4.890000E-04 +22 2.555456E-02 5.300000E-04 +23 1.471495E-02 -1.780000E-04 +S 23 +1 2.709796E+07 1.000000E-06 +2 3.592856E+06 4.000000E-06 +3 7.222638E+05 1.200000E-05 +4 1.833976E+05 3.200000E-05 +5 5.464812E+04 8.800000E-05 +6 1.830913E+04 2.390000E-04 +7 6.712501E+03 6.480000E-04 +8 2.643719E+03 1.716000E-03 +9 1.103835E+03 4.444000E-03 +10 4.836777E+02 1.066500E-02 +11 2.206245E+02 2.357300E-02 +12 1.041117E+02 4.340000E-02 +13 5.052749E+01 6.202600E-02 +14 2.485247E+01 3.059000E-02 +15 1.128102E+01 -1.053100E-01 +16 5.510166E+00 -3.334020E-01 +17 2.706914E+00 -2.967080E-01 +18 1.122807E+00 3.175360E-01 +19 5.236385E-01 6.707400E-01 +20 2.355180E-01 2.761620E-01 +21 4.157977E-02 8.694000E-03 +22 2.555456E-02 -6.306000E-03 +23 1.471495E-02 1.866000E-03 +S 1 +1 2.555456E-02 1.000000E+00 +S 22 +1 3.592856E+06 -1.000000E-06 +2 7.222638E+05 -2.000000E-06 +3 1.833976E+05 -6.000000E-06 +4 5.464812E+04 -1.700000E-05 +5 1.830913E+04 -4.600000E-05 +6 6.712501E+03 -1.250000E-04 +7 2.643719E+03 -3.330000E-04 +8 1.103835E+03 -8.590000E-04 +9 4.836777E+02 -2.070000E-03 +10 2.206245E+02 -4.559000E-03 +11 1.041117E+02 -8.459000E-03 +12 5.052749E+01 -1.203100E-02 +13 2.485247E+01 -6.123000E-03 +14 1.128102E+01 2.118800E-02 +15 5.510166E+00 6.841600E-02 +16 2.706914E+00 6.476500E-02 +17 1.122807E+00 -8.370900E-02 +18 5.236385E-01 -1.749990E-01 +19 2.355180E-01 -2.269470E-01 +20 4.157977E-02 5.221420E-01 +21 2.555456E-02 3.354410E-01 +22 1.471495E-02 2.796490E-01 +S 1 +1 1.471495E-02 1.000000E+00 +S 1 +1 8.470000E-03 1.000000E+00 +P 16 +1 9.533309E+03 8.200000E-05 +2 1.852081E+03 6.320000E-04 +3 5.464251E+02 3.611000E-03 +4 1.966947E+02 1.596000E-02 +5 7.964112E+01 5.553600E-02 +6 3.488324E+01 1.484500E-01 +7 1.604479E+01 2.883660E-01 +8 7.575632E+00 3.744160E-01 +9 3.672548E+00 2.533910E-01 +10 1.771130E+00 6.028300E-02 +11 8.515199E-01 2.503000E-03 +12 4.016017E-01 6.680000E-04 +13 1.851933E-01 -1.490000E-04 +14 5.523390E-02 4.900000E-05 +15 2.431085E-02 -3.500000E-05 +16 1.080339E-02 1.200000E-05 +P 16 +1 9.533309E+03 -2.500000E-05 +2 1.852081E+03 -1.930000E-04 +3 5.464251E+02 -1.103000E-03 +4 1.966947E+02 -4.924000E-03 +5 7.964112E+01 -1.734600E-02 +6 3.488324E+01 -4.789600E-02 +7 1.604479E+01 -9.572200E-02 +8 7.575632E+00 -1.333680E-01 +9 3.672548E+00 -7.324200E-02 +10 1.771130E+00 1.654520E-01 +11 8.515199E-01 4.055900E-01 +12 4.016017E-01 4.107060E-01 +13 1.851933E-01 1.733080E-01 +14 5.523390E-02 1.262400E-02 +15 2.431085E-02 -3.299000E-03 +16 1.080339E-02 1.050000E-03 +P 1 +1 5.523390E-02 1.000000E+00 +P 16 +1 9.533309E+03 3.000000E-06 +2 1.852081E+03 2.700000E-05 +3 5.464251E+02 1.520000E-04 +4 1.966947E+02 6.760000E-04 +5 7.964112E+01 2.393000E-03 +6 3.488324E+01 6.595000E-03 +7 1.604479E+01 1.327600E-02 +8 7.575632E+00 1.843400E-02 +9 3.672548E+00 1.010700E-02 +10 1.771130E+00 -2.662400E-02 +11 8.515199E-01 -6.080600E-02 +12 4.016017E-01 -7.470000E-02 +13 1.851933E-01 -4.857200E-02 +14 5.523390E-02 2.248180E-01 +15 2.431085E-02 5.838670E-01 +16 1.080339E-02 2.978970E-01 +P 1 +1 1.080339E-02 1.000000E+00 +P 1 +1 4.800000E-03 1.000000E+00 +D 1 +1 5.516013E-02 1.000000E+00 +D 6 +1 5.670903E+00 8.130000E-03 +2 1.459042E+00 2.288600E-02 +3 5.431827E-01 4.807800E-02 +4 1.346875E-01 7.257100E-02 +5 5.516013E-02 2.382010E-01 +6 1.501876E-02 8.262780E-01 +D 1 +1 1.501876E-02 1.000000E+00 +D 1 +1 6.010000E-03 1.000000E+00 +F 1 +1 8.493750E-02 1.000000E+00 +F 1 +1 3.398000E-02 1.000000E+00 + SCANDIUM S 20 1 2.715278E+06 8.147221E-06 diff --git a/input/methods b/input/methods index 951d243..f706c1e 100644 --- a/input/methods +++ b/input/methods @@ -15,5 +15,5 @@ # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW T F F F F F # G0T0 evGT qsGT ehG0T0 - F F F T + F F F F # * unrestricted version available diff --git a/input/options b/input/options index 71de6ce..5973b7a 100644 --- a/input/options +++ b/input/options @@ -9,7 +9,7 @@ # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 0 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W reg - 256 0.00001 T 5 T 0.0 F F F F + 256 0.00001 T 5 F 1000 F F F T # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg 10 0.00001 T 5 T 0.0 F F # ACFDT: AC Kx XBS diff --git a/int/ERI.Hu.dat b/int/ERI.Hu.dat deleted file mode 100644 index e7b1d13..0000000 --- a/int/ERI.Hu.dat +++ /dev/null @@ -1,2 +0,0 @@ -1 1 1 1 1. -2 2 2 2 1. diff --git a/int/Kin.Hu.dat b/int/Kin.Hu.dat deleted file mode 100644 index f31176f..0000000 --- a/int/Kin.Hu.dat +++ /dev/null @@ -1,2 +0,0 @@ -1 2 -1.0 -2 1 -1.0 diff --git a/int/Nuc.Hu.dat b/int/Nuc.Hu.dat deleted file mode 100644 index 48f4be0..0000000 --- a/int/Nuc.Hu.dat +++ /dev/null @@ -1,2 +0,0 @@ -1 1 0.0 -2 2 0.0 diff --git a/int/Ov.Hu.dat b/int/Ov.Hu.dat deleted file mode 100644 index f52fb30..0000000 --- a/int/Ov.Hu.dat +++ /dev/null @@ -1,2 +0,0 @@ -1 1 1.0 -2 2 1.0 diff --git a/int/x.Hu.dat b/int/x.Hu.dat deleted file mode 100644 index e69de29..0000000 diff --git a/int/y.Hu.dat b/int/y.Hu.dat deleted file mode 100644 index e69de29..0000000 diff --git a/int/z.Hu.dat b/int/z.Hu.dat deleted file mode 100644 index e69de29..0000000 diff --git a/mol/h2.xyz b/mol/h2.xyz deleted file mode 100644 index 8128fb2..0000000 --- a/mol/h2.xyz +++ /dev/null @@ -1,4 +0,0 @@ -2 - -H 0. 0. 0. -H 0. 0. 1.4 diff --git a/src/GF/G0F2.f90 b/src/GF/G0F2.f90 index 0744763..0429863 100644 --- a/src/GF/G0F2.f90 +++ b/src/GF/G0F2.f90 @@ -35,6 +35,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize double precision :: Ec double precision :: EcBSE(nspin) double precision,allocatable :: eGF2(:) + double precision,allocatable :: eGF2lin(:) double precision,allocatable :: SigC(:) double precision,allocatable :: Z(:) @@ -48,7 +49,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize ! Memory allocation - allocate(SigC(nBas),Z(nBas),eGF2(nBas)) + allocate(SigC(nBas),Z(nBas),eGF2(nBas),eGF2lin(nBas)) if(linearize) then @@ -69,13 +70,19 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,regularize end if + eGF2lin(:) = eHF(:) + Z(:)*SigC(:) + if(linearize) then - eGF2(:) = eHF(:) + Z(:)*SigC(:) + eGF2(:) = eGF2lin(:) - else + else + + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) + + call QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2) - eGF2(:) = eHF(:) + SigC(:) end if diff --git a/src/GT/G0T0.f90 b/src/GT/G0T0.f90 index 3668b25..daf8b1f 100644 --- a/src/GT/G0T0.f90 +++ b/src/GT/G0T0.f90 @@ -62,6 +62,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,ppBS double precision,allocatable :: SigX(:) double precision,allocatable :: SigT(:) double precision,allocatable :: Z(:) + double precision,allocatable :: eGTlin(:) ! Output variables @@ -91,7 +92,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,ppBS Om1aa(nVVaa),X1aa(nVVaa,nVVaa),Y1aa(nOOaa,nVVaa), & Om2aa(nOOaa),X2aa(nVVaa,nOOaa),Y2aa(nOOaa,nOOaa), & rho1aa(nBas,nBas,nVVaa),rho2aa(nBas,nBas,nOOaa), & - SigX(nBas),SigT(nBas),Z(nBas)) + SigX(nBas),SigT(nBas),Z(nBas),eGTlin(nBas)) !---------------------------------------------- ! alpha-beta block @@ -175,13 +176,22 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,ppBS ! Solve the quasi-particle equation !---------------------------------------------- +eGTlin(:) = 0d0 +eGTlin(:) = eHF(:) + Z(:)*(SigX(:) + SigT(:) - Vxc(:)) + if(linearize) then - eGT(:) = eHF(:) + Z(:)*(SigX(:) + SigT(:) - Vxc(:)) + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + + eGT(:) = eGTlin(:) else - - eGT(:) = eHF(:) + SigX(:) + SigT(:) - Vxc(:) + + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) + + call QP_graph_GT(eta,nBas,nC,nO,nV,nR,nOOaa,nVVaa,eHF,Om1aa,rho1aa,Om2aa,rho2aa,eGTlin,eGT) end if diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 4627a3d..19fbe0e 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -170,11 +170,11 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' write(*,*) - call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eGWlin,eGW) + call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eGWlin,eGW,regularize) ! Find all the roots of the QP equation if necessary - ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGWlin) + ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,OmRPA,rho_RPA,eGWlin) end if @@ -195,7 +195,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD ! Plot stuff -! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA) +! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA) ! Perform BSE calculation diff --git a/src/GW/QP_graph.f90 b/src/GW/QP_graph.f90 index 69628e8..93faf9d 100644 --- a/src/GW/QP_graph.f90 +++ b/src/GW/QP_graph.f90 @@ -1,4 +1,4 @@ -subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho,eGWlin,eGW) +subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho,eGWlin,eGW,regularize) ! Compute the graphical solution of the QP equation @@ -21,6 +21,7 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho,eGWlin,eGW) double precision,intent(in) :: rho(nBas,nBas,nS) double precision,intent(in) :: eGWlin(nBas) + logical,intent(in) :: regularize ! Local variables @@ -54,8 +55,8 @@ subroutine QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Omega,rho,eGWlin,eGW) nIt = nIt + 1 - sigC = SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) - dsigC = dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho) + sigC = SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize) + dsigC = dSigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,eHF,Omega,rho,regularize) f = w - eHF(p) - SigX(p) - sigC + Vxc(p) df = 1d0 - dsigC diff --git a/src/GW/QP_roots.f90 b/src/GW/QP_roots.f90 index fc1cde5..99a848a 100644 --- a/src/GW/QP_roots.f90 +++ b/src/GW/QP_roots.f90 @@ -21,7 +21,7 @@ subroutine QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Omega,rho,eGW) ! Local variables integer :: i,j,a,b,x,jb,g - integer,parameter :: nGrid = 100000 + integer,parameter :: nGrid = 1000000 double precision,parameter :: wmin = -50d0 double precision,parameter :: wmax = +50d0 double precision :: dw diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index 293217f..7745817 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -72,7 +72,7 @@ subroutine SRG_qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE double precision,external :: trace_matrix double precision :: dipole(ncart) - logical :: print_W = .false. + logical :: print_W = .true. double precision,allocatable :: error_diis(:,:) double precision,allocatable :: F_diis(:,:) double precision,allocatable :: OmRPA(:) diff --git a/src/GW/SigmaC.f90 b/src/GW/SigmaC.f90 index 10eddc2..70b32fb 100644 --- a/src/GW/SigmaC.f90 +++ b/src/GW/SigmaC.f90 @@ -1,4 +1,4 @@ -double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) +double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regularize) ! Compute diagonal of the correlation part of the self-energy @@ -19,32 +19,66 @@ double precision function SigmaC(p,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) + logical,intent(in) :: regularize ! Local variables integer :: i,a,jb double precision :: eps + double precision :: Dpijb,Dpajb ! Initialize SigmaC = 0d0 -! Occupied part of the correlation self-energy - do i=nC+1,nO - do jb=1,nS - eps = w - e(i) + Omega(jb) - SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2) - enddo - enddo + if (regularize) then + ! Occupied part of the correlation self-energy + do i=nC+1,nO + do jb=1,nS + eps = w - e(i) + Omega(jb) + Dpijb = e(p) - e(i) + Omega(jb) + SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*(1d0-exp(-2d0*eta*Dpijb*Dpijb))/eps + enddo + enddo + ! Virtual part of the correlation self-energy + do a=nO+1,nBas-nR + do jb=1,nS + eps = w - e(a) - Omega(jb) + Dpajb = e(p) - e(a) - Omega(jb) + SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*(1d0-exp(-2d0*eta*Dpajb*Dpajb))/eps + enddo + enddo -! Virtual part of the correlation self-energy - - do a=nO+1,nBas-nR - do jb=1,nS - eps = w - e(a) - Omega(jb) - SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2) - enddo - enddo + ! We add the static SRG term in the self-energy directly + ! do i=nC+1,nO + ! do jb=1,nS + ! Dpijb = e(p) - e(i) + Omega(jb) + ! SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*(exp(-2d0*eta*Dpijb*Dpijb)/Dpijb) + ! enddo + ! enddo + ! do a=nO+1,nBas-nR + ! do jb=1,nS + ! Dpajb = e(p) - e(a) - Omega(jb) + ! SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*(exp(-2d0*eta*Dpajb*Dpajb)/Dpajb) + ! enddo + ! enddo + + else + ! Occupied part of the correlation self-energy + do i=nC+1,nO + do jb=1,nS + eps = w - e(i) + Omega(jb) + SigmaC = SigmaC + 2d0*rho(p,i,jb)**2*eps/(eps**2 + eta**2) + enddo + enddo + ! Virtual part of the correlation self-energy + do a=nO+1,nBas-nR + do jb=1,nS + eps = w - e(a) - Omega(jb) + SigmaC = SigmaC + 2d0*rho(p,a,jb)**2*eps/(eps**2 + eta**2) + enddo + enddo + end if end function SigmaC diff --git a/src/GW/dSigmaC.f90 b/src/GW/dSigmaC.f90 index 6323de9..f0b902d 100644 --- a/src/GW/dSigmaC.f90 +++ b/src/GW/dSigmaC.f90 @@ -1,4 +1,4 @@ -double precision function dSigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) +double precision function dSigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,regularize) ! Compute the derivative of the correlation part of the self-energy @@ -19,40 +19,60 @@ double precision function dSigmaC(x,w,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho) double precision,intent(in) :: e(nBas) double precision,intent(in) :: Omega(nS) double precision,intent(in) :: rho(nBas,nBas,nS) + logical,intent(in) :: regularize ! Local variables integer :: i,j,a,b,p,jb double precision :: eps + double precision :: Dpijb,Dpajb ! Initialize dSigmaC = 0d0 -! Occupied part of the correlation self-energy + if (regularize) then + ! Occupied part of the correlation self-energy + do i=nC+1,nO + do jb=1,nS + eps = w - e(i) + Omega(jb) + Dpijb = e(p) - e(i) + Omega(jb) + dSigmaC = dSigmaC - 2d0*rho(p,i,jb)**2*(1d0-exp(-2*eta*Dpijb*Dpijb))/(eps**2) + enddo + enddo + ! Virtual part of the correlation self-energy + do a=nO+1,nBas-nR + do jb=1,nS + eps = w - e(a) - Omega(jb) + Dpajb = e(p) - e(a) - Omega(jb) + dSigmaC = dSigmaC - 2d0*rho(p,a,jb)**2*(1d0-exp(-2*eta*Dpajb*Dpajb))/(eps**2) + enddo + enddo - do i=nC+1,nO - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = w - e(i) + Omega(jb) - dSigmaC = dSigmaC - 2d0*rho(x,i,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + else + ! Occupied part of the correlation self-energy + do i=nC+1,nO + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + eps = w - e(i) + Omega(jb) + dSigmaC = dSigmaC - 2d0*rho(x,i,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo + enddo + enddo ! Virtual part of the correlation self-energy - - do a=nO+1,nBas-nR - jb = 0 - do j=nC+1,nO - do b=nO+1,nBas-nR - jb = jb + 1 - eps = w - e(a) - Omega(jb) - dSigmaC = dSigmaC - 2d0*rho(x,a,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - enddo - enddo - enddo + do a=nO+1,nBas-nR + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + eps = w - e(a) - Omega(jb) + dSigmaC = dSigmaC - 2d0*rho(x,a,jb)**2*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + enddo + enddo + enddo + end if end function dSigmaC diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 74dd127..a0fc67c 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -1,5 +1,5 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, & - singlet,triplet,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eG0W0) + 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 self-consistent eigenvalue-only GW calculation @@ -27,6 +27,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, logical,intent(in) :: ppBSE logical,intent(in) :: singlet logical,intent(in) :: triplet + logical,intent(in) :: linearize double precision,intent(in) :: eta logical,intent(in) :: regularize @@ -51,6 +52,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, integer :: ispin integer :: nSCF integer :: n_diis + integer :: i,a,jb,p double precision :: rcond double precision :: Conv double precision :: EcRPA @@ -59,6 +61,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, double precision :: EcppBSE(nspin) double precision :: EcGM double precision :: alpha + double precision :: Dpijb,Dpajb double precision,allocatable :: error_diis(:,:) double precision,allocatable :: e_diis(:,:) double precision,allocatable :: eGW(:) @@ -70,6 +73,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:) double precision,allocatable :: rho_RPA(:,:,:) + + double precision,allocatable :: eGWlin(:) integer :: nBas2 integer :: nC2 @@ -117,7 +122,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, ! Memory allocation allocate(eGW(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & - rho_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) + rho_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eGWlin(nBas)) ! Compute the exchange part of the self-energy @@ -136,6 +141,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, Z(:) = 1d0 rcond = 0d0 + + !------------------------------------------------------------------------ ! Main loop !------------------------------------------------------------------------ @@ -156,7 +163,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(regularize) then 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) + call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) else @@ -167,7 +174,25 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, ! Solve the quasi-particle equation - eGW(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:) + eGWlin(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:) + + ! Linearized or graphical solution? + + if(linearize) then + + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + + eGW(:) = eGWlin(:) + + else + + write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' + write(*,*) + + call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eGWlin,eGW,regularize) + + end if ! Convergence criteria diff --git a/src/GW/regularized_self_energy_correlation_diag.f90 b/src/GW/regularized_self_energy_correlation_diag.f90 index 7f4ea9e..545c73a 100644 --- a/src/GW/regularized_self_energy_correlation_diag.f90 +++ b/src/GW/regularized_self_energy_correlation_diag.f90 @@ -22,10 +22,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, ! Local variables integer :: i,a,p,q,jb - double precision :: eps - - double precision :: kappa - double precision :: fk + double precision :: Dpijb,Dpajb ! Output variables @@ -36,12 +33,6 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, SigC(:) = 0d0 -!-----------------------------------------! -! Parameters for regularized calculations ! -!-----------------------------------------! - - kappa = 1d0 - !----------------------------- ! COHSEX static self-energy !----------------------------- @@ -86,9 +77,8 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, 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(-2d0*eps**2/kappa**2))/eps - SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*fk + Dpijb = e(p) - e(i) + Omega(jb) + SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*(1d0 - exp(-2d0*eta*Dpijb*Dpijb))/Dpijb end do end do end do @@ -98,9 +88,8 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, 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(-2d0*eps**2/kappa**2))/eps - SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*fk + Dpajb = e(p) - e(a) - Omega(jb) + SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*(1d0 - exp(-2d0*eta*Dpajb*Dpajb))/Dpajb end do end do end do @@ -111,9 +100,7 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, do i=nC+1,nO do a=nO+1,nBas-nR do jb=1,nS - eps = e(a) - e(i) + Omega(jb) - fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps - EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk + EcGM = EcGM - 4d0*rho(a,i,jb)**2 end do end do end do diff --git a/src/GW/renormalization_factor_SRG.f90 b/src/GW/renormalization_factor_SRG.f90 index f292b18..cff5735 100644 --- a/src/GW/renormalization_factor_SRG.f90 +++ b/src/GW/renormalization_factor_SRG.f90 @@ -20,8 +20,8 @@ subroutine renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) ! Local variables - integer :: p,i,a,jb - double precision :: eps + integer :: p,i,a,m + double precision :: Dpim,Dpam ! Output variables @@ -30,9 +30,31 @@ subroutine renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) ! Initialize Z(:) = 0d0 + + do p=nC+1,nBas-nR + do i=nC+1,nO + do m=1,nS + Dpim = e(p) - e(i) + Omega(m) + Z(p) = Z(p) - 2d0*rho(p,i,m)**2*(1d0-dexp(-2d0*eta*Dpim*Dpim))/Dpim**2 + 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 m=1,nS + Dpam = e(p) - e(a) - Omega(m) + Z(p) = Z(p) - 2d0*rho(p,a,m)**2*(1d0-dexp(-2d0*eta*Dpam*Dpam))/Dpam**2 + end do + end do + end do + + ! Compute renormalization factor from derivative of SigC Z(:) = 1d0/(1d0 - Z(:)) - + end subroutine renormalization_factor_SRG diff --git a/src/HF/MOM.f90 b/src/HF/MOM.f90 index 3b720e4..6de6630 100644 --- a/src/HF/MOM.f90 +++ b/src/HF/MOM.f90 @@ -102,7 +102,7 @@ subroutine MOM(maxSCF,thresh,max_diis,nBas,nO,S,T,V,Hc,ERI,X,ENuc,ERHF,c,e,P) call Coulomb_matrix_AO_basis(nBas,P,ERI,J) call exchange_matrix_AO_basis(nBas,P,ERI,K) - F(:,:) = Hc(:,:) + J(:,:) + K(:,:) + F(:,:) = Hc(:,:) + J(:,:) + 0.5*K(:,:) ! Check convergence diff --git a/src/HF/MOM_overlap.f90 b/src/HF/MOM_overlap.f90 index f4f4c61..3ca91a0 100644 --- a/src/HF/MOM_overlap.f90 +++ b/src/HF/MOM_overlap.f90 @@ -26,7 +26,7 @@ subroutine MOM_overlap(nBas,nO,S,cG,c,ON) do i=1,nBas do j=1,nBas - pOv(j) = pOv(j) + ON(i)*Ov(i,j)**2 + pOv(j) = pOv(j) + Ov(i,j)**2 enddo enddo @@ -41,7 +41,8 @@ subroutine MOM_overlap(nBas,nO,S,cG,c,ON) ploc = maxloc(pOv,nBas) ON(ploc) = 1d0 pOv(ploc) = 0d0 - enddo + enddo + ! print*,'--- Occupation numbers ---' ! call matout(nBas,1,ON) diff --git a/src/LR/print_excitation.f90 b/src/LR/print_excitation.f90 index 898f9f3..acf7045 100644 --- a/src/LR/print_excitation.f90 +++ b/src/LR/print_excitation.f90 @@ -14,7 +14,7 @@ subroutine print_excitation(method,ispin,nS,Omega) ! Local variables character*14 :: spin_manifold - integer,parameter :: maxS = 20 + integer,parameter :: maxS = 50 integer :: ia if(ispin == 1) spin_manifold = 'singlet' diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index b89450c..4327923 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -335,7 +335,7 @@ program QuAcK else - call RMOM(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, & + call MOM(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, & nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc) end if @@ -1017,7 +1017,7 @@ program QuAcK else call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX, & - BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet,eta_GW,regGW, & + BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet,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 call cpu_time(end_evGW) @@ -1094,7 +1094,7 @@ program QuAcK if(doufG0W0) then call cpu_time(start_ufGW) - call ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,TDA_W) call cpu_time(end_ufGW) t_ufGW = end_ufGW - start_ufGW diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index f24743f..355f71f 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -205,8 +205,11 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) call dgemm ('N', 'N', nOO+nVV, nVV, nOO+nVV, 1d0, M, nOO+nVV, Z1, nOO+nVV, 0d0, tmp1, nOO+nVV) call dgemm ('T', 'N', nVV , nVV, nOO+nVV, 1d0, Z1, nOO+nVV, tmp1, nOO+nVV, 0d0, S1, nVV) - - S2 = - matmul(transpose(Z2),matmul(M,Z2)) + !S1 = + matmul(transpose(Z1),matmul(M,Z1)) + + call dgemm ('N', 'N', nOO+nVV, nOO, nOO+nVV, 1d0, M, nOO+nVV, -1d0*Z2, nOO+nVV, 0d0, tmp2, nOO+nVV) + call dgemm ('T', 'N', nOO , nOO, nOO+nVV, 1d0, Z2, nOO+nVV, tmp2, nOO+nVV, 0d0, S2, nOO) + ! S2 = - matmul(transpose(Z2),matmul(M,Z2)) if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1) if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2) @@ -214,8 +217,11 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) write (*,*) 'OK SO FAR' - call dgemm ('N', 'N', nVV, nVV, nVV, 1d0, Z1, nVV, O1, nVV, 0d0, Z1, nVV) - Z2 = matmul(Z2,O2) + !Z1 = matmul(Z1,O1) + call dgemm ('N', 'N', nOO+nVV,nVV,nVV, 1d0, Z1, nOO+nVV, O1, nVV,0d0, tmp1, nOO+nVV) + Z1 = tmp1 + call dgemm ('N', 'N', nOO+nVV,nOO,nOO, 1d0, Z2, nOO+nVV, O2, nOO,0d0, tmp2, nOO+nVV) + Z2 = tmp2 ! Define submatrices X1, Y1, X2, & Y2 diff --git a/src/make_ninja.py b/src/make_ninja.py index fd31a89..fa39fdf 100755 --- a/src/make_ninja.py +++ b/src/make_ninja.py @@ -41,6 +41,17 @@ CXX = g++ LAPACK=-lblas -llapack STDCXX=-lstdc++ FIX_ORDER_OF_LIBS=-Wl,--start-group +""" + + compile_ifort_linux = """ +FC = ifort -mkl=parallel -qopenmp +AR = ar crs +FFLAGS = -I$IDIR -g -Ofast -traceback +CC = icc +CXX = icpc +LAPACK= +STDCXX=-lstdc++ +FIX_ORDER_OF_LIBS=-Wl,--start-group """ else: compile_gfortran_mac = """ @@ -67,7 +78,8 @@ FIX_ORDER_OF_LIBS=-Wl,--start-group if sys.platform in ["linux", "linux2"]: - compiler = compile_gfortran_linux +# compiler = compile_gfortran_linux + compiler = compile_ifort_linux elif sys.platform == "darwin": compiler = compile_gfortran_mac else: diff --git a/src/utils/read_basis.f90 b/src/utils/read_basis.f90 index 75a9f49..88f64b5 100644 --- a/src/utils/read_basis.f90 +++ b/src/utils/read_basis.f90 @@ -52,131 +52,135 @@ subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KSh !------------------------------------------------------------------------ ! Loop over atoms !------------------------------------------------------------------------ - do i=1,nNuc + ! do i=1,nNuc - read(2,*) iNuc,nShAt - write(*,'(A28,1X,I16)') 'Atom n. ',iNuc - write(*,'(A28,1X,I16)') 'number of shells ',nShAt - write(*,'(A28)') '------------------' +! read(2,*) iNuc,nShAt +! write(*,'(A28,1X,I16)') 'Atom n. ',iNuc +! write(*,'(A28,1X,I16)') 'number of shells ',nShAt +! write(*,'(A28)') '------------------' -!------------------------------------------------------------------------ -! Loop over shells -!------------------------------------------------------------------------ - do j=1,nShAt +! !------------------------------------------------------------------------ +! ! Loop over shells +! !------------------------------------------------------------------------ +! do j=1,nShAt - nShell = nShell + 1 +! nShell = nShell + 1 - ! Basis function centers +! ! Basis function centers - do k=1,ncart - CenterShell(nShell,k) = rNuc(iNuc,k) - enddo +! do k=1,ncart +! CenterShell(nShell,k) = rNuc(iNuc,k) +! enddo - ! Shell type and contraction degree +! ! Shell type and contraction degree - read(2,*) shelltype,KShell(nShell) +! read(2,*) shelltype,KShell(nShell) - select case (shelltype) - case ("S") +! select case (shelltype) +! case ("S") - TotAngMomShell(nShell) = 0 - write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell) +! TotAngMomShell(nShell) = 0 +! write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell) - case ("P") +! case ("P") - TotAngMomShell(nShell) = 1 - write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) +! TotAngMomShell(nShell) = 1 +! write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) - case ("D") +! case ("D") - TotAngMomShell(nShell) = 2 - write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell) +! TotAngMomShell(nShell) = 2 +! write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell) - case ("F") +! case ("F") - TotAngMomShell(nShell) = 3 - write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell) +! TotAngMomShell(nShell) = 3 +! write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell) - case ("G") +! case ("G") - TotAngMomShell(nShell) = 4 - write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell) +! TotAngMomShell(nShell) = 4 +! write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell) - case ("H") +! case ("H") - TotAngMomShell(nShell) = 5 - write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell) +! TotAngMomShell(nShell) = 5 +! write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell) - case ("I") +! case ("I") - TotAngMomShell(nShell) = 6 - write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) +! TotAngMomShell(nShell) = 6 +! write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) - case ("J") +! case ("J") - TotAngMomShell(nShell) = 7 - write(*,'(A28,1X,I16)') 'j-type shell with K = ',KShell(nShell) +! TotAngMomShell(nShell) = 7 +! write(*,'(A28,1X,I16)') 'j-type shell with K = ',KShell(nShell) - case default +! case default - call print_warning('!!! Angular momentum too high !!!') - stop +! call print_warning('!!! Angular momentum too high !!!') +! stop - end select +! end select -! Read exponents and contraction coefficients +! ! Read exponents and contraction coefficients - write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' - do k=1,Kshell(nShell) - read(2,*) kk,ExpShell(nShell,k),DShell(nShell,k) - write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) - enddo +! write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' +! do k=1,Kshell(nShell) +! read(2,*) kk,ExpShell(nShell,k),DShell(nShell,k) +! write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) +! enddo - min_exponent(iNuc,TotAngMomShell(nShell)+1) & - = min(min_exponent(iNuc,TotAngMomShell(nShell)+1),minval(ExpShell(nShell,1:KShell(nShell)))) - max_exponent(iNuc) = max(max_exponent(iNuc),maxval(ExpShell(nShell,:))) - max_ang_mom(iNuc) = max(max_ang_mom(iNuc),TotAngMomShell(nShell)) +! min_exponent(iNuc,TotAngMomShell(nShell)+1) & +! = min(min_exponent(iNuc,TotAngMomShell(nShell)+1),minval(ExpShell(nShell,1:KShell(nShell)))) +! max_exponent(iNuc) = max(max_exponent(iNuc),maxval(ExpShell(nShell,:))) +! max_ang_mom(iNuc) = max(max_ang_mom(iNuc),TotAngMomShell(nShell)) - enddo -!------------------------------------------------------------------------ -! End loop over shells -!------------------------------------------------------------------------ +! enddo +! !------------------------------------------------------------------------ +! ! End loop over shells +! !------------------------------------------------------------------------ - write(*,'(A28)') '------------------' +! write(*,'(A28)') '------------------' -! print*,'maximum angular momemtum for atom n. ',iNuc,' = ' -! print*,max_ang_mom(iNuc) -! print*,'maximum exponent for atom n. ',iNuc,' = ' -! print*,max_exponent(iNuc) -! print*,'minimum exponent for atom n. ',iNuc,' = ' -! print*,min_exponent(iNuc,1:max_ang_mom(iNuc)+1) +! ! print*,'maximum angular momemtum for atom n. ',iNuc,' = ' +! ! print*,max_ang_mom(iNuc) +! ! print*,'maximum exponent for atom n. ',iNuc,' = ' +! ! print*,max_exponent(iNuc) +! ! print*,'minimum exponent for atom n. ',iNuc,' = ' +! ! print*,min_exponent(iNuc,1:max_ang_mom(iNuc)+1) - enddo -!------------------------------------------------------------------------ -! End loop over atoms -!------------------------------------------------------------------------ +! enddo +! !------------------------------------------------------------------------ +! ! End loop over atoms +! !------------------------------------------------------------------------ -! Total number of shells +! ! Total number of shells - write(*,'(A28,1X,I16)') 'Number of shells',nShell - write(*,'(A28)') '------------------' - write(*,*) +! write(*,'(A28,1X,I16)') 'Number of shells',nShell +! write(*,'(A28)') '------------------' +! write(*,*) -! Close file with basis set specification +! ! Close file with basis set specification - close(unit=2) +! close(unit=2) ! Calculate number of basis functions - nBas = 0 - do iShell=1,nShell - nBas = nBas + (TotAngMomShell(iShell)*TotAngMomShell(iShell) + 3*TotAngMomShell(iShell) + 2)/2 - enddo + ! nBas = 0 + ! do iShell=1,nShell + ! nBas = nBas + (TotAngMomShell(iShell)*TotAngMomShell(iShell) + 3*TotAngMomShell(iShell) + 2)/2 + ! enddo - write(*,'(A28)') '------------------' - write(*,'(A28,1X,I16)') 'Number of basis functions',NBas - write(*,'(A28)') '------------------' - write(*,*) +! write(*,'(A28)') '------------------' +! write(*,'(A28,1X,I16)') 'Number of basis functions',NBas +! write(*,'(A28)') '------------------' +! write(*,*) + + open(unit=3,file='int/nBas.dat') + read(3,*) nBas + close(unit=3) ! Number of virtual orbitals