From a61fd351f99a1de8c6e9e960a66f6cc5ef2bdb36 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 17 Jun 2020 22:10:08 +0200 Subject: [PATCH] fix option mismatch --- input/methods | 4 +-- input/options | 8 ++--- src/QuAcK/BSE2_B_matrix_dynamic.f90 | 12 +------ src/QuAcK/BSE2_dynamic_perturbation.f90 | 10 ++---- src/QuAcK/Bethe_Salpeter.f90 | 29 +++------------ src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 | 4 +-- .../Bethe_Salpeter_ZA_matrix_dynamic.f90 | 4 +-- .../Bethe_Salpeter_dynamic_perturbation.f90 | 36 +++++++++---------- src/QuAcK/G0W0.f90 | 5 ++- src/QuAcK/read_options.f90 | 8 ++--- 10 files changed, 44 insertions(+), 76 deletions(-) diff --git a/input/methods b/input/methods index c48703f..63c03d4 100644 --- a/input/methods +++ b/input/methods @@ -11,9 +11,9 @@ # RPA RPAx ppRPA F F F # G0F2 evGF2 G0F3 evGF3 - T F F F + F F F F # G0W0 evGW qsGW - F F F + T F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index add4505..ea311e6 100644 --- a/input/options +++ b/input/options @@ -6,10 +6,10 @@ 64 0.0000001 T 5 # spin: singlet triplet TDA T T F -# GF: maxSCF thresh DIIS n_diis lin eta renorm - 256 0.00001 T 5 T 0.0 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 T F F +# 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.00367493 F F T F F # ACFDT: AC Kx XBS F F T # BSE: BSE dBSE dTDA evDyn diff --git a/src/QuAcK/BSE2_B_matrix_dynamic.f90 b/src/QuAcK/BSE2_B_matrix_dynamic.f90 index 8df510a..aefbea7 100644 --- a/src/QuAcK/BSE2_B_matrix_dynamic.f90 +++ b/src/QuAcK/BSE2_B_matrix_dynamic.f90 @@ -1,4 +1,4 @@ -subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_dyn,ZB_dyn) +subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_dyn) ! Compute the anti-resonant part of the dynamic BSE2 matrix @@ -24,12 +24,10 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_ ! Output variables double precision,intent(out) :: B_dyn(nS,nS) - double precision,intent(out) :: ZB_dyn(nS,nS) ! Initialization B_dyn(:,:) = 0d0 - ZB_dyn(:,:) = 0d0 ! Second-order correlation kernel for the block A of the singlet manifold @@ -53,14 +51,12 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_ - ERI(b,k,c,i)*ERI(a,c,j,k) + 2d0*ERI(b,k,c,i)*ERI(a,c,k,j) B_dyn(ia,jb) = B_dyn(ia,jb) - num*dem/(dem**2 + eta**2) - ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 dem = + eGF(i) - eGF(c) + eGF(k) - eGF(b) num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) & - ERI(b,c,k,i)*ERI(a,k,j,c) + 2d0*ERI(b,c,k,i)*ERI(a,k,c,j) B_dyn(ia,jb) = B_dyn(ia,jb) - num*dem/(dem**2 + eta**2) - ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 end do end do @@ -73,7 +69,6 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_ - ERI(a,b,d,c)*ERI(c,d,i,j) + 2d0*ERI(a,b,d,c)*ERI(c,d,j,i) B_dyn(ia,jb) = B_dyn(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) - ZB_dyn(ia,jb) = ZB_dyn(ia,jb) - 0.5d0*num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 end do end do @@ -86,7 +81,6 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_ - ERI(a,b,l,k)*ERI(k,l,i,j) + 2d0*ERI(a,b,l,k)*ERI(k,l,j,i) B_dyn(ia,jb) = B_dyn(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) - ZB_dyn(ia,jb) = ZB_dyn(ia,jb) - 0.5d0*num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 end do end do @@ -120,13 +114,11 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_ num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) - ERI(b,k,c,i)*ERI(a,c,j,k) B_dyn(ia,jb) = B_dyn(ia,jb) - num*dem/(dem**2 + eta**2) - ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 dem = + eGF(i) - eGF(c) + eGF(k) - eGF(b) num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) - ERI(b,c,k,i)*ERI(a,k,j,c) B_dyn(ia,jb) = B_dyn(ia,jb) - num*dem/(dem**2 + eta**2) - ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 end do end do @@ -138,7 +130,6 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_ num = ERI(a,b,c,d)*ERI(c,d,j,i) + ERI(a,b,d,c)*ERI(c,d,i,j) B_dyn(ia,jb) = B_dyn(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) - ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + 0.5d0*num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 end do end do @@ -150,7 +141,6 @@ subroutine BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_ num = ERI(a,b,k,l)*ERI(k,l,j,i) + ERI(a,b,l,k)*ERI(k,l,i,j) B_dyn(ia,jb) = B_dyn(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) - ZB_dyn(ia,jb) = ZB_dyn(ia,jb) + 0.5d0*num*(dem**2 - eta**2)/(dem**2 + eta**2)**2 end do end do diff --git a/src/QuAcK/BSE2_dynamic_perturbation.f90 b/src/QuAcK/BSE2_dynamic_perturbation.f90 index 703684f..9a69147 100644 --- a/src/QuAcK/BSE2_dynamic_perturbation.f90 +++ b/src/QuAcK/BSE2_dynamic_perturbation.f90 @@ -41,13 +41,12 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF, double precision,allocatable :: ZAm_dyn(:,:) double precision,allocatable :: B_dyn(:,:) - double precision,allocatable :: ZB_dyn(:,:) ! Memory allocation allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) - if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),B_dyn(nS,nS),ZB_dyn(nS,nS)) + if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),B_dyn(nS,nS)) ! Print main components of transition vectors @@ -82,14 +81,11 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,eHF, ! Second part of the resonant and anti-resonant part of the BSE correction (frequency independent) call BSE2_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,-OmBSE(ia),Am_dyn,ZAm_dyn) - ZAm_dyn(:,:) = - ZAm_dyn(:,:) - call BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn,ZB_dyn) + call BSE2_B_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_dyn) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & - - dot_product(Y,matmul(ZAm_dyn,Y)) & - + dot_product(X,matmul(ZB_dyn,Y)) & - - dot_product(Y,matmul(ZB_dyn,X)) + + dot_product(Y,matmul(ZAm_dyn,Y)) OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & - dot_product(Y,matmul(Am_dyn,Y)) & diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 index 88ad8fe..d01f978 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -29,13 +29,11 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man ! Local variables - logical :: W_BSE = .false. integer :: ispin integer :: isp_W double precision,allocatable :: OmBSE(:,:) double precision,allocatable :: XpY_BSE(:,:,:) double precision,allocatable :: XmY_BSE(:,:,:) - double precision,allocatable :: rho_BSE(:,:,:,:) ! Output variables @@ -45,7 +43,6 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man ! Memory allocation allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) - if(W_BSE) allocate(rho_BSE(nBas,nBas,nS,nspin)) !------------------- ! Singlet manifold @@ -77,8 +74,6 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man if(dBSE) then - if(W_BSE) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) - ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) if(evDyn) then @@ -87,14 +82,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) else - if(W_BSE) then - call Bethe_Salpeter_dynamic_perturbation(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(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) - end if - + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) end if end if @@ -131,24 +120,16 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_man if(dBSE) then - if(W_BSE) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_BSE(:,:,ispin),rho_BSE(:,:,:,ispin)) - ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) if(evDyn) then - call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & + call Bethe_Salpeter_dynamic_perturbation_iterative(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(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(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW(:),OmRPA(:,ispin),OmBSE(:,ispin), & - XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) - end if - + call Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA(:,ispin),OmBSE(:,ispin), & + XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin),rho_RPA(:,:,:,ispin)) end if end if diff --git a/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 index 1d8f891..61dfa01 100644 --- a/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 @@ -57,10 +57,10 @@ subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,Om chi = 0d0 do kc=1,maxS - eps = OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) + eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2) - eps = OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) + eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) chi = chi + rho(i,j,kc)*rho(a,b,kc)*eps/(eps**2 + eta**2) enddo diff --git a/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 index 3e2f245..eb71cb4 100644 --- a/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 +++ b/src/QuAcK/Bethe_Salpeter_ZA_matrix_dynamic.f90 @@ -48,10 +48,10 @@ subroutine Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O chi = 0d0 do kc=1,maxS - eps = OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) + eps = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(j)) chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - eps = OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) + eps = + OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) chi = chi + rho(i,j,kc)*rho(a,b,kc)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 enddo diff --git a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 index 71028b9..d01e611 100644 --- a/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/QuAcK/Bethe_Salpeter_dynamic_perturbation.f90 @@ -84,42 +84,40 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW, ! Resonant part of the BSE correction for dynamical TDA - call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & - Ap_dyn(:,:)) + call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,Ap_dyn) ! Renormalization factor of the resonant parts for dynamical TDA - call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & - ZAp_dyn(:,:)) + call Bethe_Salpeter_ZA_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho,ZAp_dyn) - ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) - OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) + ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) + OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X)) else ! Resonant and anti-resonant part of the BSE correction - call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & - Ap_dyn(:,:),Am_dyn(:,:),Bp_dyn(:,:),Bm_dyn(:,:)) + call Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, & + Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) ! Renormalization factor of the resonant and anti-resonant parts - call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW(:),OmRPA(:),OmBSE(ia),rho(:,:,:), & - ZAp_dyn(:,:),ZAm_dyn(:,:),ZBp_dyn(:,:),ZBm_dyn(:,:)) + call Bethe_Salpeter_ZAB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,OmRPA,OmBSE(ia),rho, & + ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn) - ZDyn(ia) = dot_product(X(:),matmul(ZAp_dyn(:,:),X(:))) & - - dot_product(Y(:),matmul(ZAm_dyn(:,:),Y(:))) & - + dot_product(X(:),matmul(ZBp_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(ZBm_dyn(:,:),X(:))) + ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & + - dot_product(Y,matmul(ZAm_dyn,Y)) & + + dot_product(X,matmul(ZBp_dyn,Y)) & + - dot_product(Y,matmul(ZBm_dyn,X)) - OmDyn(ia) = dot_product(X(:),matmul(Ap_dyn(:,:),X(:))) & - - dot_product(Y(:),matmul(Am_dyn(:,:),Y(:))) & - + dot_product(X(:),matmul(Bp_dyn(:,:),Y(:))) & - - dot_product(Y(:),matmul(Bm_dyn(:,:),X(:))) + OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & + - dot_product(Y,matmul(Am_dyn,Y)) & + + dot_product(X,matmul(Bp_dyn,Y)) & + - dot_product(Y,matmul(Bm_dyn,X)) end if - ZDyn(ia) = 1d0/(1d0 - ZDyn(ia)) + ZDyn(ia) = 1d0/(1d0 - ZDyn(ia)) OmDyn(ia) = ZDyn(ia)*OmDyn(ia) write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & diff --git a/src/QuAcK/G0W0.f90 b/src/QuAcK/G0W0.f90 index a330cdf..96b12dc 100644 --- a/src/QuAcK/G0W0.f90 +++ b/src/QuAcK/G0W0.f90 @@ -66,6 +66,10 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, write(*,*)'************************************************' write(*,*) +! Initialization + + EcRPA(:) = 0d0 + ! SOSEX correction if(SOSEX) write(*,*) 'SOSEX correction activated!' @@ -148,7 +152,6 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1) diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 005cd22..d4d80e8 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -147,13 +147,13 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t thresh_GW = 1d-5 DIIS_GW = .false. n_diis_GW = 5 + linGW = .false. + eta_GW = 0d0 COHSEX = .false. SOSEX = .false. TDA_W = .false. G0W = .false. GW0 = .false. - linGW = .false. - eta_GW = 0d0 read(1,*) read(1,*) maxSCF_GW,thresh_GW,answer1,n_diis_GW,answer2,eta_GW, & @@ -163,9 +163,9 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t if(answer2 == 'T') linGW = .true. if(answer3 == 'T') COHSEX = .true. if(answer4 == 'T') SOSEX = .true. - if(answer5 == 'T') TDA_W = .true. - if(answer6 == 'T') G0W = .true. + if(answer5 == 'T') G0W = .true. if(answer7 == 'T') GW0 = .true. + if(answer7 == 'T') TDA_W = .true. if(.not.DIIS_GW) n_diis_GW = 1 ! Options for adiabatic connection