From f7e90a9f99819770b3b5dd3ec014c11ecf9dcf61 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 14 Jun 2020 13:04:16 +0200 Subject: [PATCH] spin in W --- input/basis | 69 +++++++++++++++++-- input/molecule | 5 +- input/molecule.xyz | 5 +- input/options | 4 +- src/QuAcK/ACFDT.f90 | 7 +- src/QuAcK/AOtoMO_integral_transform.f90 | 2 +- src/QuAcK/Bethe_Salpeter.f90 | 22 +++--- .../Bethe_Salpeter_ZAB_matrix_dynamic.f90 | 12 ++-- .../Bethe_Salpeter_dynamic_perturbation.f90 | 15 +++- src/QuAcK/G0T0.f90 | 5 +- src/QuAcK/G0W0.f90 | 7 +- src/QuAcK/QuAcK.f90 | 20 +++--- src/QuAcK/evGT.f90 | 5 +- src/QuAcK/evGW.f90 | 5 +- src/QuAcK/qsGW.f90 | 5 +- src/QuAcK/read_options.f90 | 31 +++++---- 16 files changed, 152 insertions(+), 67 deletions(-) diff --git a/input/basis b/input/basis index 0a75001..80356bb 100644 --- a/input/basis +++ b/input/basis @@ -1,7 +1,64 @@ -1 2 -S 3 - 1 38.4216340 0.0237660 - 2 5.7780300 0.1546790 - 3 1.2417740 0.4696300 +1 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.2979640 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.185000E-01 1.000000E+00 +D 1 +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 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.185000E-01 1.000000E+00 +D 1 +1 8.170000E-01 1.0000000 diff --git a/input/molecule b/input/molecule index c78e87e..76ebcdf 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,5 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 2 7 7 0 0 # Znuc x y z - He 0.0 0.0 0.0 + N 0. 0. -1.04008632 + N 0. 0. +1.04008632 diff --git a/input/molecule.xyz b/input/molecule.xyz index 797b5fc..e1773f0 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,4 @@ - 1 + 2 - He 0.0000000000 0.0000000000 0.0000000000 + N 0.0000000000 0.0000000000 -0.5503900175 + N 0.0000000000 0.0000000000 0.5503900175 diff --git a/input/options b/input/options index ae799d8..b4c3817 100644 --- a/input/options +++ b/input/options @@ -8,8 +8,8 @@ T T # 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_W TDA G0W GW0 lin eta - 256 0.00001 T 5 F F T T T F F T 0.00367493 +# GW/GT: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA_W TDA dTDA G0W GW0 lin eta + 256 0.00001 T 5 F F T F F 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/ACFDT.f90 b/src/QuAcK/ACFDT.f90 index 89b77b2..1b0c42e 100644 --- a/src/QuAcK/ACFDT.f90 +++ b/src/QuAcK/ACFDT.f90 @@ -32,6 +32,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl ! Local variables integer :: ispin + integer :: isp_W integer :: iAC double precision :: lambda double precision,allocatable :: Ec(:,:) @@ -62,6 +63,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl if(singlet_manifold) then ispin = 1 + isp_W = 1 write(*,*) '--------------' write(*,*) 'Singlet states' @@ -78,7 +80,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl if(doXBS) then - call linear_response(ispin,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & + call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) @@ -108,6 +110,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl if(triplet_manifold) then ispin = 2 + isp_W = 1 write(*,*) '--------------' write(*,*) 'Triplet states' @@ -124,7 +127,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet_manifold,tripl if(doXBS) then - call linear_response(1,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & + call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) diff --git a/src/QuAcK/AOtoMO_integral_transform.f90 b/src/QuAcK/AOtoMO_integral_transform.f90 index 8246fe7..17d4119 100644 --- a/src/QuAcK/AOtoMO_integral_transform.f90 +++ b/src/QuAcK/AOtoMO_integral_transform.f90 @@ -73,7 +73,7 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis) do nu=1,nBas ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j)*scr(i,nu,k,l) enddo - print*,i,k,j,l,ERI_MO_basis(i,j,k,l) +! print*,i,k,j,l,ERI_MO_basis(i,j,k,l) enddo enddo enddo diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 index 52e28b4..5b827ac 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & +subroutine Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,EcRPA,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -10,6 +10,7 @@ subroutine Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dTDA logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold @@ -29,6 +30,7 @@ subroutine Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & logical :: evDyn = .false. logical :: W_BSE = .false. integer :: ispin + integer :: isp_W double precision,allocatable :: OmBSE(:,:) double precision,allocatable :: XpY_BSE(:,:,:) double precision,allocatable :: XmY_BSE(:,:,:) @@ -51,11 +53,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & if(singlet_manifold) then ispin = 1 + isp_W = 1 EcBSE(ispin) = 0d0 ! Compute (singlet) RPA screening - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin)) @@ -75,15 +78,15 @@ subroutine Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation_iterative(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) else if(W_BSE) then - call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) end if @@ -98,11 +101,12 @@ subroutine Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & if(triplet_manifold) then ispin = 2 + isp_W = 1 EcBSE(ispin) = 0d0 ! Compute (singlet) RPA screening - call linear_response(1,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA(:,:,ispin),rho_RPA(:,:,:,ispin)) @@ -122,15 +126,15 @@ subroutine Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation_iterative(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) else if(W_BSE) then - call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) else - call Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) end if diff --git a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 index 9ed4416..14cd120 100644 --- a/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_ZAB_matrix_dynamic.f90 @@ -68,10 +68,10 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, eps_Ap = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) chi_Ap = chi_Ap + rho(i,j,kc)*rho(a,b,kc)*(eps_Ap**2 - eta**2)/(eps_Ap**2 + eta**2)**2 - eps_Am = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) + eps_Am = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 - eps_Am = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) + eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) chi_Am = chi_Am + rho(i,j,kc)*rho(a,b,kc)*(eps_Am**2 - eta**2)/(eps_Am**2 + eta**2)**2 eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) @@ -80,19 +80,19 @@ subroutine Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW, eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) chi_Bp = chi_Bp + rho(i,b,kc)*rho(a,j,kc)*(eps_Bp**2 - eta**2)/(eps_Bp**2 + eta**2)**2 - eps_Bm = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) + eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 - eps_Bm = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) + eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) chi_Bm = chi_Bm + rho(i,b,kc)*rho(a,j,kc)*(eps_Bm**2 - eta**2)/(eps_Bm**2 + eta**2)**2 enddo ZAp(ia,jb) = ZAp(ia,jb) + 2d0*lambda*chi_Ap - ZAm(ia,jb) = ZAm(ia,jb) + 2d0*lambda*chi_Am + ZAm(ia,jb) = ZAm(ia,jb) - 2d0*lambda*chi_Am ZBp(ia,jb) = ZBp(ia,jb) + 2d0*lambda*chi_Bp - ZBm(ia,jb) = ZBm(ia,jb) + 2d0*lambda*chi_Bm + ZBm(ia,jb) = ZBm(ia,jb) - 2d0*lambda*chi_Bm enddo enddo diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 index 9113519..214cbf2 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) +subroutine Bethe_Salpeter_dynamic_perturbation(TDA,dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,OmBSE,XpY,XmY,rho) ! Compute dynamical effects via perturbation theory for BSE @@ -8,6 +8,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O ! Input variables logical,intent(in) :: TDA + logical,intent(in) :: dTDA double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -27,7 +28,6 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O integer :: ia - logical :: dTDA = .true. integer,parameter :: maxS = 10 double precision :: gapGW @@ -54,10 +54,16 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS)) - ! Print main components of transition vectors +! Print main components of transition vectors call print_transition_vectors(nBas,nC,nO,nV,nR,nS,OmBSE,XpY,XmY) + if(dTDA) then + write(*,*) + write(*,*) '*** dynamical TDA activated ***' + write(*,*) + end if + gapGW = eGW(nO+1) - eGW(nO) write(*,*) '---------------------------------------------------------------------------------------------------' @@ -73,6 +79,9 @@ subroutine Bethe_Salpeter_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,eGW,O X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) +! call matout(nS,1,X) +! call matout(nS,1,Y) + ! First-order correction if(dTDA) then diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index 48a2e98..c36d140 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -1,4 +1,4 @@ -subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,singlet_manifold,triplet_manifold, & +subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold, & linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -14,6 +14,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,singlet_manifold,tri logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dTDA logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold logical,intent(in) :: linearize @@ -211,7 +212,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,singlet_manifold,tri allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin)) - call Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & + call Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI,eHF,eG0T0,Omega,XpY,XmY,rho,EcRPA,EcBSE) if(exchange_kernel) then diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index 53f9105..4c3e87d 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -1,5 +1,5 @@ -subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & - singlet_manifold,triplet_manifold,linearize,eta, & +subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dTDA, & + singlet_manifold,triplet_manifold,linearize,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eGW) ! Perform G0W0 calculation @@ -18,6 +18,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dTDA logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold logical,intent(in) :: linearize @@ -163,7 +164,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, & if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & + call Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) if(exchange_kernel) then diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index bfc5be4..d1c3472 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -101,7 +101,7 @@ program QuAcK integer :: maxSCF_GW,n_diis_GW double precision :: thresh_GW - logical :: DIIS_GW,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,G0W,GW0,linGW + logical :: DIIS_GW,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,dTDA,G0W,GW0,linGW double precision :: eta_GW integer :: nMC,nEq,nWalk,nPrint,iSeed @@ -147,7 +147,7 @@ program QuAcK 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_GW,TDA_W,TDA_GW,G0W,GW0,linGW,eta_GW, & + COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,dTDA,G0W,GW0,linGW,eta_GW, & doACFDT,exchange_kernel,doXBS, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) @@ -671,8 +671,8 @@ program QuAcK if(doG0W0) then call cpu_time(start_G0W0) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW, & - singlet_manifold,triplet_manifold,linGW,eta_GW, & + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,dTDA, & + singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) @@ -690,7 +690,7 @@ program QuAcK call cpu_time(start_evGW) call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & - BSE_GW,TDA_W,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & + BSE_GW,TDA_W,TDA_GW,dTDA,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_evGW) @@ -708,7 +708,7 @@ program QuAcK call cpu_time(start_qsGW) call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX, & - BSE_GW,TDA_W,TDA_GW,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & + BSE_GW,TDA_W,TDA_GW,dTDA,G0W,GW0,singlet_manifold,triplet_manifold,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF) call cpu_time(end_qsGW) @@ -727,7 +727,7 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) - call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_W,TDA_GW, & + call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_W,TDA_GW,dTDA, & singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_G0T0) @@ -746,7 +746,7 @@ program QuAcK call cpu_time(start_evGT) call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS, & - BSE_GW,TDA_W,TDA_GW,singlet_manifold,triplet_manifold,eta_GW, & + BSE_GW,TDA_W,TDA_GW,dTDA,singlet_manifold,triplet_manifold,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0T0) call cpu_time(end_evGT) @@ -862,7 +862,7 @@ program QuAcK ! Long-range G0W0 calculation call cpu_time(start_G0W0) - call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW, & + call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,dTDA, & singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0) call cpu_time(end_G0W0) @@ -876,7 +876,7 @@ program QuAcK ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:) call cpu_time(start_G0T0) - call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_W,TDA_GW, & + call G0T0(doACFDT,exchange_kernel,doXBS,BSE_GW,TDA_W,TDA_GW,dTDA, & singlet_manifold,triplet_manifold,linGW,eta_GW, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0) call cpu_time(end_G0T0) diff --git a/src/QuAcK/evGT.f90 b/src/QuAcK/evGT.f90 index 4166c85..436379a 100644 --- a/src/QuAcK/evGT.f90 +++ b/src/QuAcK/evGT.f90 @@ -1,5 +1,5 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & - BSE,TDA_W,TDA,singlet_manifold,triplet_manifold, & + BSE,TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold, & eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0) ! Perform eigenvalue self-consistent calculation with a T-matrix self-energy (evGT) @@ -18,6 +18,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dTDA logical,intent(in) :: singlet_manifold logical,intent(in) :: triplet_manifold double precision,intent(in) :: eta @@ -259,7 +260,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin),rho(nBas,nBas,nS,nspin)) - call Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & + call Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI,eGT,eGT,Omega,XpY,XmY,rho,EcRPA,EcBSE) if(exchange_kernel) then diff --git a/src/QuAcK/evGW.f90 b/src/QuAcK/evGW.f90 index 240c3ed..2d25f4e 100644 --- a/src/QuAcK/evGW.f90 +++ b/src/QuAcK/evGW.f90 @@ -1,4 +1,4 @@ -subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,G0W,GW0, & +subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA,dTDA,G0W,GW0, & singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI,PHF,cHF,eHF,eG0W0) ! Perform self-consistent eigenvalue-only GW calculation @@ -21,6 +21,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dTDA logical,intent(in) :: G0W logical,intent(in) :: GW0 logical,intent(in) :: singlet_manifold @@ -227,7 +228,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & + call Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) write(*,*) diff --git a/src/QuAcK/qsGW.f90 b/src/QuAcK/qsGW.f90 index 0339249..0e61ea7 100644 --- a/src/QuAcK/qsGW.f90 +++ b/src/QuAcK/qsGW.f90 @@ -1,5 +1,5 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & - COHSEX,SOSEX,BSE,TDA_W,TDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, & + COHSEX,SOSEX,BSE,TDA_W,TDA,dTDA,G0W,GW0,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -20,6 +20,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, logical,intent(in) :: BSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA + logical,intent(in) :: dTDA logical,intent(in) :: G0W logical,intent(in) :: GW0 logical,intent(in) :: singlet_manifold @@ -270,7 +271,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, if(BSE) then - call Bethe_Salpeter(TDA_W,TDA,singlet_manifold,triplet_manifold,eta, & + call Bethe_Salpeter(TDA_W,TDA,dTDA,singlet_manifold,triplet_manifold,eta, & nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) write(*,*) diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index cfd2f7c..9f4bc05 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -4,7 +4,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t 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_GW,TDA_W,TDA_GW,G0W,GW0,linGW,eta_GW, & + COHSEX,SOSEX,BSE_GW,TDA_W,TDA_GW,dTDA,G0W,GW0,linGW,eta_GW, & doACFDT,exchange_kernel,doXBS, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) @@ -48,6 +48,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t logical,intent(out) :: BSE_GW logical,intent(out) :: TDA_W logical,intent(out) :: TDA_GW + logical,intent(out) :: dTDA logical,intent(out) :: G0W logical,intent(out) :: GW0 logical,intent(out) :: linGW @@ -67,7 +68,8 @@ 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,answer8,answer9 + character(len=1) :: answer1,answer2,answer3,answer4,answer5, & + answer6,answer7,answer8,answer9,answer10 ! Open file with method specification @@ -151,24 +153,27 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t BSE_GW = .false. TDA_W = .false. TDA_GW = .false. + dTDA = .false. G0W = .false. GW0 = .false. linGW = .false. eta_GW = 0d0 read(1,*) - read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2, & - answer3,answer4,answer5,answer6,answer7,answer8,answer9,eta_GW + read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2, & + answer3,answer4,answer5,answer6,answer7,answer8, & + answer9,answer10,eta_GW - if(answer1 == 'T') DIIS_GW = .true. - if(answer2 == 'T') COHSEX = .true. - if(answer3 == 'T') SOSEX = .true. - if(answer4 == 'T') BSE_GW = .true. - if(answer5 == 'T') TDA_W = .true. - if(answer6 == 'T') TDA_GW = .true. - if(answer7 == 'T') G0W = .true. - if(answer8 == 'T') GW0 = .true. - if(answer9 == 'T') linGW = .true. + if(answer1 == 'T') DIIS_GW = .true. + if(answer2 == 'T') COHSEX = .true. + if(answer3 == 'T') SOSEX = .true. + if(answer4 == 'T') BSE_GW = .true. + if(answer5 == 'T') TDA_W = .true. + if(answer6 == 'T') TDA_GW = .true. + if(answer7 == 'T') dTDA = .true. + if(answer8 == 'T') G0W = .true. + if(answer9 == 'T') GW0 = .true. + if(answer10 == 'T') linGW = .true. if(.not.DIIS_GW) n_diis_GW = 1 ! Options for adiabatic connection