diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 02d5058..02087a0 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -2,4 +2,4 @@ 2 1 1 0 0 # Znuc x y z H 0. 0. 0. - H 0. 0. 1.435 + H 0. 0. 10.1 diff --git a/examples/molecule.H2O b/examples/molecule.H2O index 12a9ac0..d22e22e 100644 --- a/examples/molecule.H2O +++ b/examples/molecule.H2O @@ -1,10 +1,6 @@ # nAt nEla nElb nCore nRyd 3 5 5 0 0 # Znuc x y z -O 0.0 0.0 0.226200217 -H 0.0 1.43921542 -0.904422923 -H 0.0 -1.43921542 -0.904422923 - O 0. 0. -0.13209669 H 0. 1.43152878 0.97970006 H 0. -1.43152878 0.97970006 diff --git a/input/basis b/input/basis index 93f1bae..80356bb 100644 --- a/input/basis +++ b/input/basis @@ -1,73 +1,64 @@ -1 15 +1 6 S 9 - 1 33980.0000000 0.0000910 - 2 5089.0000000 0.0007040 - 3 1157.0000000 0.0036930 - 4 326.6000000 0.0153600 - 5 106.1000000 0.0529290 - 6 38.1100000 0.1470430 - 7 14.7500000 0.3056310 - 8 6.0350000 0.3993450 - 9 2.5300000 0.2170510 +1 9.046000E+03 7.000000E-04 +2 1.357000E+03 5.389000E-03 +3 3.093000E+02 2.740600E-02 +4 8.773000E+01 1.032070E-01 +5 2.856000E+01 2.787230E-01 +6 1.021000E+01 4.485400E-01 +7 3.838000E+00 2.782380E-01 +8 7.466000E-01 1.544000E-02 +9 2.248000E-01 -2.864000E-03 S 9 - 1 33980.0000000 -0.0000190 - 2 5089.0000000 -0.0001510 - 3 1157.0000000 -0.0007850 - 4 326.6000000 -0.0033240 - 5 106.1000000 -0.0115120 - 6 38.1100000 -0.0341600 - 7 14.7500000 -0.0771730 - 8 6.0350000 -0.1414930 - 9 2.5300000 -0.1180190 +1 9.046000E+03 -1.530000E-04 +2 1.357000E+03 -1.208000E-03 +3 3.093000E+02 -5.992000E-03 +4 8.773000E+01 -2.454400E-02 +5 2.856000E+01 -6.745900E-02 +6 1.021000E+01 -1.580780E-01 +7 3.838000E+00 -1.218310E-01 +8 7.466000E-01 5.490030E-01 +9 2.248000E-01 5.788150E-01 S 1 - 1 0.7355000 1.0000000 -S 1 - 1 0.2905000 1.0000000 -S 1 - 1 0.1111000 1.0000000 -P 3 - 1 34.5100000 0.0053780 - 2 7.9150000 0.0361320 - 3 2.3680000 0.1424930 +1 2.248000E-01 1.000000E+00 +P 4 +1 1.355000E+01 3.991900E-02 +2 2.917000E+00 2.171690E-01 +3 7.973000E-01 5.103190E-01 +4 2.185000E-01 4.622140E-01 P 1 - 1 0.8132000 1.0000000 -P 1 - 1 0.2890000 1.0000000 -P 1 - 1 0.1007000 1.0000000 +1 2.185000E-01 1.000000E+00 D 1 - 1 1.8480000 1.0000000 -D 1 - 1 0.6490000 1.0000000 -D 1 - 1 0.2280000 1.0000000 -F 1 - 1 1.4190000 1.0000000 -F 1 - 1 0.4850000 1.0000000 -G 1 - 1 1.0110000 1.0000000 -2 10 -S 3 - 1 82.6400000 0.0020060 - 2 12.4100000 0.0153430 - 3 2.8240000 0.0755790 +1 8.170000E-01 1.0000000 +2 6 +S 9 +1 9.046000E+03 7.000000E-04 +2 1.357000E+03 5.389000E-03 +3 3.093000E+02 2.740600E-02 +4 8.773000E+01 1.032070E-01 +5 2.856000E+01 2.787230E-01 +6 1.021000E+01 4.485400E-01 +7 3.838000E+00 2.782380E-01 +8 7.466000E-01 1.544000E-02 +9 2.248000E-01 -2.864000E-03 +S 9 +1 9.046000E+03 -1.530000E-04 +2 1.357000E+03 -1.208000E-03 +3 3.093000E+02 -5.992000E-03 +4 8.773000E+01 -2.454400E-02 +5 2.856000E+01 -6.745900E-02 +6 1.021000E+01 -1.580780E-01 +7 3.838000E+00 -1.218310E-01 +8 7.466000E-01 5.490030E-01 +9 2.248000E-01 5.788150E-01 S 1 - 1 0.7977000 1.0000000 -S 1 - 1 0.2581000 1.0000000 -S 1 - 1 0.0898900 1.0000000 +1 2.248000E-01 1.000000E+00 +P 4 +1 1.355000E+01 3.991900E-02 +2 2.917000E+00 2.171690E-01 +3 7.973000E-01 5.103190E-01 +4 2.185000E-01 4.622140E-01 P 1 - 1 2.2920000 1.0000000 -P 1 - 1 0.8380000 1.0000000 -P 1 - 1 0.2920000 1.0000000 +1 2.185000E-01 1.000000E+00 D 1 - 1 2.0620000 1.0000000 -D 1 - 1 0.6620000 1.0000000 -F 1 - 1 1.3970000 1.0000000 - +1 8.170000E-01 1.0000000 diff --git a/input/methods b/input/methods index 296048d..c48703f 100644 --- a/input/methods +++ b/input/methods @@ -7,9 +7,9 @@ # drCCD rCCD lCCD pCCD F F F F # CIS CID CISD - T F F + F F F # RPA RPAx ppRPA - F T F + F F F # G0F2 evGF2 G0F3 evGF3 T F F F # G0W0 evGW qsGW diff --git a/input/molecule b/input/molecule index 8715e4f..76ebcdf 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd - 2 3 3 0 0 + 2 7 7 0 0 # Znuc x y z - C 0. 0. 0. - H 0. 0. 2.13709127 + N 0. 0. -1.04008632 + N 0. 0. +1.04008632 diff --git a/input/molecule.xyz b/input/molecule.xyz index a901cdf..e1773f0 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,4 @@ 2 - C 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 1.1309000791 + N 0.0000000000 0.0000000000 -0.5503900175 + N 0.0000000000 0.0000000000 0.5503900175 diff --git a/input/options b/input/options index cd0edbf..1ed8348 100644 --- a/input/options +++ b/input/options @@ -6,10 +6,10 @@ 64 0.0000001 T 5 # spin: singlet triplet T T -# GF: maxSCF thresh DIIS n_diis lin renorm - 256 0.00001 T 5 T 3 +# GF: maxSCF thresh DIIS n_diis lin renorm BSE TDA eta + 256 0.00001 T 5 T 3 T T 0.00367493 # GW/GT: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta - 256 0.00001 T 10 F F T T F F T 0.0 + 256 0.00001 T 5 F F T F F F T 0.00367493 # ACFDT: AC Kx XBS F F T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/src/QuAcK/BSE2_dynamic_perturbation.f90 b/src/QuAcK/BSE2_dynamic_perturbation.f90 index f24b13e..b06fbbe 100644 --- a/src/QuAcK/BSE2_dynamic_perturbation.f90 +++ b/src/QuAcK/BSE2_dynamic_perturbation.f90 @@ -47,11 +47,17 @@ subroutine BSE2_dynamic_perturbation(ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,O if(.not.dTDA) allocate(B_dyn(nS,nS),ZB_dyn(nS,nS)) + ! Print main components of transition vectors + + call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) + gapGF = eGF(nO+1) - eGF(nO) write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) ' First-order dynamical correction to static 2nd-order Bethe-Salpeter excitation energies ' write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(A58,F10.6,A3)') ' BSE neutral excitation must be lower than the GF2 gap = ',gapGF*HaToeV,' eV' + write(*,*) '---------------------------------------------------------------------------------------------------' write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' @@ -95,8 +101,6 @@ subroutine BSE2_dynamic_perturbation(ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF,eGF,O write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,ZDyn(ia) - if(OmBSE(ia) > gapGF) write(*,*) ' !!! BSE2 neutral excitation larger than the GF2 gap !!! ' - end do write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 index ee5cd3f..07d1fc2 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 @@ -25,8 +25,9 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O ! Local variables - logical :: dTDA = .true. integer :: ia + + logical :: dTDA = .true. integer,parameter :: maxS = 10 double precision :: gapGW @@ -47,11 +48,17 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O if(.not.dTDA) allocate(B_dyn(nS,nS),ZB_dyn(nS,nS)) + ! Print main components of transition vectors + + call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) + gapGW = eGW(nO+1) - eGW(nO) write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies ' write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GW gap = ',gapGW*HaToeV,' eV' + write(*,*) '---------------------------------------------------------------------------------------------------' write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' @@ -107,8 +114,6 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,ZDyn(ia) - if(OmBSE(ia) > gapGW) write(*,*) ' !!! BSE neutral excitation larger than the GW gap !!! ' - end do write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) diff --git a/src/QuAcK/G0F2.f90 b/src/QuAcK/G0F2.f90 index ce8d8ff..3d0905a 100644 --- a/src/QuAcK/G0F2.f90 +++ b/src/QuAcK/G0F2.f90 @@ -66,8 +66,8 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC, eps = eHF(p) + eHF(a) - eHF(i) - eHF(j) V = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j) - Sig(p) = Sig(p) + V/eps - Z(p) = Z(p) + V/eps**2 + Sig(p) = Sig(p) + V*eps/(eps**2 + eta**2) + Z(p) = Z(p) - V*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do @@ -81,15 +81,15 @@ subroutine G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linearize,eta,nBas,nC, eps = eHF(p) + eHF(i) - eHF(a) - eHF(b) V = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b) - Sig(p) = Sig(p) + V/eps - Z(p) = Z(p) + V/eps**2 + Sig(p) = Sig(p) + V*eps/(eps**2 + eta**2) + Z(p) = Z(p) - V*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do end do end do - Z(:) = 1d0/(1d0 + Z(:)) + Z(:) = 1d0/(1d0 - Z(:)) if(linearize) then diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 62f7723..718f0b6 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -96,12 +96,13 @@ program QuAcK integer :: maxSCF_GF,n_diis_GF,renormGF double precision :: thresh_GF - logical :: DIIS_GF,linGF + logical :: DIIS_GF,linGF,BSE_GF,TDA_GF + double precision :: eta_GF integer :: maxSCF_GW,n_diis_GW double precision :: thresh_GW - logical :: DIIS_GW,COHSEX,SOSEX,BSE,TDA,G0W,GW0,linGW - double precision :: eta + logical :: DIIS_GW,COHSEX,SOSEX,BSE_GW,TDA_GW,G0W,GW0,linGW + double precision :: eta_GW integer :: nMC,nEq,nWalk,nPrint,iSeed double precision :: dt @@ -144,8 +145,9 @@ program QuAcK maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & singlet_manifold,triplet_manifold, & maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,renormGF, & + BSE_GF,TDA_GF,eta_GF, & maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW, & - COHSEX,SOSEX,BSE,TDA,G0W,GW0,linGW,eta, & + COHSEX,SOSEX,BSE_GW,TDA_GW,G0W,GW0,linGW,eta_GW, & doACFDT,exchange_kernel,doXBS, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) @@ -533,7 +535,7 @@ program QuAcK if(doRPA) then call cpu_time(start_RPA) - call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & + call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_RPA) @@ -550,7 +552,7 @@ program QuAcK if(doRPAx) then call cpu_time(start_RPAx) - call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, & + call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,0d0, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_RPAx) @@ -601,8 +603,8 @@ program QuAcK if(doG0F2) then call cpu_time(start_GF2) - call G0F2(BSE,TDA,singlet_manifold,triplet_manifold,linGF, & - eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF) + call G0F2(BSE_GF,TDA_GF,singlet_manifold,triplet_manifold,linGF, & + eta_GF,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 @@ -618,8 +620,8 @@ program QuAcK if(doevGF2) then call cpu_time(start_GF2) - call evGF2(BSE,TDA,maxSCF_GF,thresh_GF,n_diis_GF,singlet_manifold,triplet_manifold,linGF, & - eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF) + call evGF2(BSE_GF,TDA_GF,maxSCF_GF,thresh_GF,n_diis_GF,singlet_manifold,triplet_manifold,linGF, & + eta_GF,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF) call cpu_time(end_GF2) t_GF2 = end_GF2 - start_GF2 @@ -669,8 +671,8 @@ program QuAcK if(doG0W0) then call cpu_time(start_G0W0) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, & - singlet_manifold,triplet_manifold,linGW,eta, & + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_GW, & + singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) @@ -687,8 +689,8 @@ program QuAcK if(doevGW) then call cpu_time(start_evGW) - call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, & - singlet_manifold,triplet_manifold,eta, & + call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & + BSE_GW,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_evGW) @@ -705,8 +707,8 @@ program QuAcK if(doqsGW) then call cpu_time(start_qsGW) - call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, & - singlet_manifold,triplet_manifold,eta, & + call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & + BSE_GW,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF) call cpu_time(end_qsGW) @@ -725,8 +727,8 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) - call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, & - singlet_manifold,triplet_manifold,linGW,eta, & + call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_GW, & + singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_G0T0) @@ -743,8 +745,9 @@ 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,singlet_manifold,triplet_manifold, & - eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0) + call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & + BSE_GW,TDA_GW,singlet_manifold,triplet_manifold,eta_GW, & + nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_evGT) t_evGT = end_evGT - start_evGT @@ -859,8 +862,8 @@ program QuAcK ! Long-range G0W0 calculation call cpu_time(start_G0W0) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, & - singlet_manifold,triplet_manifold,linGW,eta, & + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_GW, & + singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) @@ -873,8 +876,8 @@ program QuAcK ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:) call cpu_time(start_G0T0) - call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, & - singlet_manifold,triplet_manifold,linGW,eta, & + call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_GW, & + singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0) call cpu_time(end_G0T0) diff --git a/src/QuAcK/evGF2.f90 b/src/QuAcK/evGF2.f90 index 98fb616..fa99903 100644 --- a/src/QuAcK/evGF2.f90 +++ b/src/QuAcK/evGF2.f90 @@ -1,5 +1,5 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold,linearize, & - eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,V,e0) + eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,V,eHF) ! Perform eigenvalue self-consistent second-order Green function calculation @@ -25,7 +25,7 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol integer,intent(in) :: nS double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: e0(nBas) + double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: V(nBas,nBas,nBas,nBas) ! Local variables @@ -63,8 +63,8 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol n_diis = 0 e_diis(:,:) = 0d0 error_diis(:,:) = 0d0 - eGF2(:) = e0(:) - eOld(:) = e0(:) + eGF2(:) = eHF(:) + eOld(:) = eHF(:) !------------------------------------------------------------------------ ! Main SCF loop @@ -81,10 +81,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol do j=nC+1,nO do a=nO+1,nBas-nR - eps = eGF2(p) + e0(a) - e0(i) - e0(j) + eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j) Sig(p) = Sig(p) & - + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps + + (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)*eps/(eps**2 + eta**2) end do end do @@ -96,10 +96,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol do a=nO+1,nBas-nR do b=nO+1,nBas-nR - eps = eGF2(p) + e0(i) - e0(a) - e0(b) + eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b) Sig(p) = Sig(p) & - + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps + + (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)*eps/(eps**2 + eta**2) end do end do @@ -115,9 +115,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol do j=nC+1,nO do a=nO+1,nBas-nR - eps = eGF2(p) + e0(a) - e0(i) - e0(j) + eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j) - Z(p) = Z(p) - (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)/eps**2 + Z(p) = Z(p) & + - (2d0*V(p,a,i,j) - V(p,a,j,i))*V(p,a,i,j)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do @@ -129,9 +130,10 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol do a=nO+1,nBas-nR do b=nO+1,nBas-nR - eps = eGF2(p) + e0(i) - e0(a) - e0(b) + eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b) - Z(p) = Z(p) - (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)/eps**2 + Z(p) = Z(p) & + - (2d0*V(p,i,a,b) - V(p,i,b,a))*V(p,i,a,b)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do end do @@ -142,11 +144,11 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol if(linearize) then - eGF2(:) = e0(:) + Z(:)*Sig(:) + eGF2(:) = eHF(:) + Z(:)*Sig(:) else - eGF2(:) = e0(:) + Sig(:) + eGF2(:) = eHF(:) + Sig(:) end if @@ -154,7 +156,7 @@ subroutine evGF2(BSE,TDA,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifol ! Print results - call print_evGF2(nBas,nO,nSCF,Conv,e0,eGF2) + call print_evGF2(nBas,nO,nSCF,Conv,eHF,eGF2) ! DIIS extrapolation diff --git a/src/QuAcK/print_transition_vectors.f90 b/src/QuAcK/print_transition_vectors.f90 new file mode 100644 index 0000000..0dfaf8c --- /dev/null +++ b/src/QuAcK/print_transition_vectors.f90 @@ -0,0 +1,66 @@ +subroutine print_transition_vectors(nBas,nC,nO,nV,nR,nS,Omega,XpY,XmY) + +! Print transition vectors for linear response calculation + + implicit none + include 'parameters.h' + +! Input variables + + 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) :: Omega(nS) + double precision,intent(in) :: XpY(nS,nS) + double precision,intent(in) :: XmY(nS,nS) + +! Local variables + + integer :: ia,jb,i,j,a,b + integer,parameter :: maxS = 10 + double precision,parameter :: thres_vec = 0.1d0 + double precision,allocatable :: X(:) + double precision,allocatable :: Y(:) + +! Memory allocation + + allocate(X(nS),Y(nS)) + + write(*,*) + do ia=1,min(nS,maxS) + + X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) + Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) + + print*,'--------------------------------' + write(*,'(A15,I3,A2,F10.6,A3)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV' + print*,'--------------------------------' + print*,'Resonant vectors' + print*,'--------------------------------' + + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A3,F10.6)') j,' --> ',b,' = ',X(jb) + end do + end do + + print*,'--------------------------------' + print*,'Anti-resonant vectors' + print*,'--------------------------------' + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A3,F10.6)') j,' --> ',b,' = ',Y(jb) + end do + end do + write(*,*) + + end do + +end subroutine print_transition_vectors diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 8ac01f2..3beec7c 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -2,8 +2,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & singlet_manifold,triplet_manifold, & maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,renormGF, & + BSE_GF,TDA_GF,eta_GF, & maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW, & - COHSEX,SOSEX,BSE,TDA,G0W,GW0,linGW,eta, & + COHSEX,SOSEX,BSE_GW,TDA_GW,G0W,GW0,linGW,eta_GW, & doACFDT,exchange_kernel,doXBS, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) @@ -34,6 +35,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t integer,intent(out) :: n_diis_GF logical,intent(out) :: linGF integer,intent(out) :: renormGF + logical,intent(out) :: BSE_GF + logical,intent(out) :: TDA_GF + double precision,intent(out) :: eta_GF integer,intent(out) :: maxSCF_GW double precision,intent(out) :: thresh_GW @@ -41,12 +45,12 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t integer,intent(out) :: n_diis_GW logical,intent(out) :: COHSEX logical,intent(out) :: SOSEX - logical,intent(out) :: BSE - logical,intent(out) :: TDA + logical,intent(out) :: BSE_GW + logical,intent(out) :: TDA_GW logical,intent(out) :: G0W logical,intent(out) :: GW0 logical,intent(out) :: linGW - double precision,intent(out) :: eta + double precision,intent(out) :: eta_GW logical,intent(out) :: doACFDT logical,intent(out) :: exchange_kernel @@ -120,14 +124,19 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t thresh_GF = 1d-5 DIIS_GF = .false. n_diis_GF = 5 - linGF = .false. - renormGF = 0 + linGF = .false. + renormGF = 0 + BSE_GF = .false. + TDA_GF = .false. + eta_GF = 0d0 read(1,*) - read(1,*) maxSCF_GF,thresh_GW,answer1,n_diis_GF,answer2,renormGF + read(1,*) maxSCF_GF,thresh_GW,answer1,n_diis_GF,answer2,renormGF,answer3,answer4,eta_GF if(answer1 == 'T') DIIS_GF = .true. if(answer2 == 'T') linGF = .true. + if(answer3 == 'T') BSE_GF = .true. + if(answer4 == 'T') TDA_GF = .true. if(.not.DIIS_GF) n_diis_GF = 1 ! Read GW options @@ -138,25 +147,25 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t n_diis_GW = 5 COHSEX = .false. SOSEX = .false. - BSE = .false. - TDA = .false. + BSE_GW = .false. + TDA_GW = .false. G0W = .false. GW0 = .false. linGW = .false. - eta = 0d0 + eta_GW = 0d0 read(1,*) read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2, & - answer3,answer4,answer5,answer6,answer7,answer8,eta + answer3,answer4,answer5,answer6,answer7,answer8,eta_GW if(answer1 == 'T') DIIS_GW = .true. - if(answer2 == 'T') COHSEX = .true. - if(answer3 == 'T') SOSEX = .true. - if(answer4 == 'T') BSE = .true. - if(answer5 == 'T') TDA = .true. - if(answer6 == 'T') G0W = .true. - if(answer7 == 'T') GW0 = .true. - if(answer8 == 'T') linGW = .true. + if(answer2 == 'T') COHSEX = .true. + if(answer3 == 'T') SOSEX = .true. + if(answer4 == 'T') BSE_GW = .true. + if(answer5 == 'T') TDA_GW = .true. + if(answer6 == 'T') G0W = .true. + if(answer7 == 'T') GW0 = .true. + if(answer8 == 'T') linGW = .true. if(.not.DIIS_GW) n_diis_GW = 1 ! Options for adiabatic connection