diff --git a/input/methods b/input/methods index 26ff8ab..89fb8f0 100644 --- a/input/methods +++ b/input/methods @@ -9,11 +9,11 @@ # CIS* CIS(D) CID CISD F F F F # RPA* RPAx* ppRPA - T F F -# G0F2 evGF2 G0F3 evGF3 - F F F F + F F F +# G0F2 evGF2 qsGF2 G0F3 evGF3 + F F T F F # G0W0* evGW* qsGW* - F F T + F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/mol/benzene.xyz b/mol/benzene.xyz index 4b49a6a..08d5f47 100644 --- a/mol/benzene.xyz +++ b/mol/benzene.xyz @@ -1,5 +1,5 @@ 12 -Benzene,^1A_{1g},CC3,aug-cc-pVTZ + C 0.00000000 1.39250319 0.00000000 C -1.20594314 0.69625160 0.00000000 C -1.20594314 -0.69625160 0.00000000 @@ -11,4 +11,4 @@ H -2.14171677 -1.23652075 0.00000000 H 0.00000000 -2.47304151 0.00000000 H 2.14171677 -1.23652075 0.00000000 H 2.14171677 1.23652075 0.00000000 -H 0.00000000 2.47304151 0.00000000 \ No newline at end of file +H 0.00000000 2.47304151 0.00000000 diff --git a/src/MBPT/G0F2.f90 b/src/MBPT/G0F2.f90 index 8204ad1..db22b4b 100644 --- a/src/MBPT/G0F2.f90 +++ b/src/MBPT/G0F2.f90 @@ -1,5 +1,4 @@ -subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & - linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) +subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform a one-shot second-order Green function calculation @@ -13,8 +12,8 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta integer,intent(in) :: nBas @@ -116,7 +115,7 @@ subroutine G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold, & if(BSE) then - call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) + call BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) end if diff --git a/src/MBPT/evGF2.f90 b/src/MBPT/evGF2.f90 index 168a762..14bd86f 100644 --- a/src/MBPT/evGF2.f90 +++ b/src/MBPT/evGF2.f90 @@ -1,4 +1,4 @@ -subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold,triplet_manifold, & +subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, & linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform eigenvalue self-consistent second-order Green function calculation @@ -16,8 +16,8 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold integer,intent(in) :: maxSCF double precision,intent(in) :: thresh integer,intent(in) :: max_diis - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet logical,intent(in) :: linearize double precision,intent(in) :: eta integer,intent(in) :: nBas @@ -169,7 +169,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet_manifold if(BSE) then - call BSE2(TDA,dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) + call BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF2,EcBSE) end if diff --git a/src/MBPT/print_qsGF2.f90 b/src/MBPT/print_qsGF2.f90 new file mode 100644 index 0000000..c78100b --- /dev/null +++ b/src/MBPT/print_qsGF2.f90 @@ -0,0 +1,122 @@ +subroutine print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigC,Z,EqsGF2,dipole) + +! Print one-electron energies and other stuff for qsGF2 + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nO + integer,intent(in) :: nSCF + double precision,intent(in) :: ENuc + double precision,intent(in) :: Conv + double precision,intent(in) :: thresh + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: eGF2(nBas) + double precision,intent(in) :: c(nBas) + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas),V(nBas,nBas) + double precision,intent(in) :: J(nBas,nBas),K(nBas,nBas),F(nBas,nBas) + double precision,intent(in) :: Z(nBas),SigC(nBas,nBas) + double precision,intent(in) :: dipole(ncart) + +! Local variables + + integer :: q,ixyz,HOMO,LUMO + double precision :: Gap,ET,EV,EJ,Ex,Ec + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: EqsGF2 + +! HOMO and LUMO + + HOMO = nO + LUMO = HOMO + 1 + Gap = eGF2(LUMO)-eGF2(HOMO) + +! Compute energies + + ET = trace_matrix(nBas,matmul(P,T)) + EV = trace_matrix(nBas,matmul(P,V)) + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) + Ec = 0.50d0*trace_matrix(nBas,matmul(P,SigC)) + EqsGF2 = ET + EV + EJ + Ex + Ec + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + if(nSCF < 10) then + write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + else + write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + endif + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & + '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + do q=1,nBas + write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & + '|',q,'|',eHF(q)*HaToeV,'|',SigC(q,q)*HaToeV,'|',Z(q),'|',eGF2(q)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A10,I3)') 'Iteration ',nSCF + write(*,'(2X,A19,F15.5)')'max(|FPS - SPF|) = ',Conv + write(*,*)'-------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') 'qsGF2 HOMO energy:',eGF2(HOMO)*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'qsGF2 LUMO energy:',eGF2(LUMO)*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'qsGF2 HOMO-LUMO gap :',Gap*HaToeV,' eV' + write(*,*)'-------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') ' qsGF2 total energy:',EqsGF2 + ENuc,' au' + write(*,'(2X,A30,F15.6,A3)') ' qsGF2 exchange energy:',Ex,' au' + write(*,*)'-------------------------------------------' + write(*,*) + +! Dump results for final iteration + + if(Conv < thresh) then + + write(*,*) + write(*,'(A50)') '---------------------------------------' + write(*,'(A32)') ' Summary ' + write(*,'(A50)') '---------------------------------------' + write(*,'(A32,1X,F16.10,A3)') ' One-electron energy: ',ET + EV,' au' + write(*,'(A32,1X,F16.10,A3)') ' Kinetic energy: ',ET,' au' + write(*,'(A32,1X,F16.10,A3)') ' Potential energy: ',EV,' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A32,1X,F16.10,A3)') ' Two-electron energy: ',EJ + Ex,' au' + write(*,'(A32,1X,F16.10,A3)') ' Hartree energy: ',EJ,' au' + write(*,'(A32,1X,F16.10,A3)') ' Exchange energy: ',Ex,' au' + write(*,'(A32,1X,F16.10,A3)') ' Correlation energy: ',Ec,' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EqsGF2,' au' + write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' + write(*,'(A32,1X,F16.10,A3)') ' qsGF2 energy: ',ENuc + EqsGF2,' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A35)') ' Dipole moment (Debye) ' + write(*,'(10X,4A10)') 'X','Y','Z','Tot.' + write(*,'(10X,4F10.6)') (dipole(ixyz)*auToD,ixyz=1,ncart),norm2(dipole)*auToD + write(*,'(A50)') '-----------------------------------------' + write(*,*) + + write(*,'(A50)') '---------------------------------------' + write(*,'(A32)') ' qsGF2 MO coefficients' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,nBas,c) + write(*,*) + write(*,'(A50)') '---------------------------------------' + write(*,'(A32)') ' qsGF2 MO energies' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,1,eGF2) + write(*,*) + + endif + + +end subroutine print_qsGF2 diff --git a/src/MBPT/qsGF2.f90 b/src/MBPT/qsGF2.f90 new file mode 100644 index 0000000..0248685 --- /dev/null +++ b/src/MBPT/qsGF2.f90 @@ -0,0 +1,222 @@ +subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, & + eta,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF, & + S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + +! Perform a quasiparticle self-consistent GF2 calculation + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + logical,intent(in) :: BSE + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: evDyn + logical,intent(in) :: singlet + logical,intent(in) :: triplet + double precision,intent(in) :: eta + + integer,intent(in) :: nNuc + double precision,intent(in) :: ZNuc(nNuc) + double precision,intent(in) :: rNuc(nNuc,ncart) + double precision,intent(in) :: ENuc + + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: cHF(nBas,nBas) + double precision,intent(in) :: PHF(nBas,nBas) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nBas,nBas,ncart) + +! Local variables + + integer :: nSCF + integer :: nBasSq + integer :: ispin + integer :: n_diis + double precision :: EqsGF2 + double precision :: Conv + double precision :: rcond + double precision,external :: trace_matrix + double precision :: dipole(ncart) + double precision :: EcBSE(nspin) + + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: F_diis(:,:) + double precision,allocatable :: c(:,:) + double precision,allocatable :: cp(:,:) + double precision,allocatable :: eGF2(:) + double precision,allocatable :: P(:,:) + double precision,allocatable :: F(:,:) + double precision,allocatable :: Fp(:,:) + double precision,allocatable :: J(:,:) + double precision,allocatable :: K(:,:) + double precision,allocatable :: SigC(:,:) + double precision,allocatable :: SigCp(:,:) + double precision,allocatable :: SigCm(:,:) + double precision,allocatable :: Z(:) + double precision,allocatable :: error(:,:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| Self-consistent qsGF2 calculation |' + write(*,*)'************************************************' + write(*,*) + +! Warning + + write(*,*) '!! ERIs in MO basis will be overwritten in qsGF2 !!' + write(*,*) + +! Stuff + + nBasSq = nBas*nBas + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Memory allocation + + allocate(eGF2(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & + J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & + error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + +! Initialization + + nSCF = -1 + n_diis = 0 + ispin = 1 + Conv = 1d0 + P(:,:) = PHF(:,:) + eGF2(:) = eHF(:) + c(:,:) = cHF(:,:) + F_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF <= maxSCF) + + ! Increment + + nSCF = nSCF + 1 + + ! Buid Coulomb matrix + + call Coulomb_matrix_AO_basis(nBas,P,ERI_AO,J) + + ! Compute exchange part of the self-energy + + call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) + + ! AO to MO transformation of two-electron integrals + + call AOtoMO_integral_transform(1,1,1,1,nBas,c,ERI_AO,ERI_MO) + + ! Compute self-energy and renormalization factor + + call self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI_MO,SigC,Z) + + ! Make correlation self-energy Hermitian and transform it back to AO basis + + SigCp = 0.5d0*(SigC + transpose(SigC)) + SigCm = 0.5d0*(SigC - transpose(SigC)) + + call MOtoAO_transform(nBas,S,c,SigCp) + + ! Solve the quasi-particle equation + + F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:) + + ! Compute commutator and convergence criteria + + error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) + Conv = maxval(abs(error)) + + ! DIIS extrapolation + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + +! Reset DIIS if required + + if(abs(rcond) < 1d-15) n_diis = 0 + + ! Diagonalize Hamiltonian in AO basis + + Fp = matmul(transpose(X),matmul(F,X)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nBas,cp,eGF2) + c = matmul(X,cp) + + ! Compute new density matrix in the AO basis + + P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + + ! Print results + + call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsGF2(nBas,nO,nSCF,Conv,thresh,eHF,eGF2,c,ENuc,P,T,V,J,K,F,SigCp,Z,EqsGF2,dipole) + + enddo +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF+1) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +! Deallocate memory + + deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,error,error_diis,F_diis) + +! Perform BSE calculation + + if(BSE) then + + call BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGF2,EcBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy (triplet) =',EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 correlation energy =',EcBSE(1) + EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGF2 total energy =',ENuc + EqsGF2 + EcBSE(1) + EcBSE(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + +end subroutine qsGF2 diff --git a/src/MBPT/self_energy_GF2.f90 b/src/MBPT/self_energy_GF2.f90 new file mode 100644 index 0000000..e4dcdd1 --- /dev/null +++ b/src/MBPT/self_energy_GF2.f90 @@ -0,0 +1,73 @@ +subroutine self_energy_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI,SigC,Z) + +! Compute GF2 self-energy and its renormalization factor + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: eGF2(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: i,j,a,b + integer :: p,q + double precision :: eps + double precision :: num + +! Output variables + + double precision,intent(out) :: SigC(nBas,nBas) + double precision,intent(out) :: Z(nBas) + +! Initialize + + SigC(:,:) = 0d0 + Z(:) = 0d0 + +! Compute GF2 self-energy + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do i=nC+1,nO + do j=nC+1,nO + do a=nO+1,nBas-nR + + eps = eGF2(p) + eHF(a) - eHF(i) - eHF(j) + num = (2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j) + + SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + end do + end do + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do i=nC+1,nO + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + + eps = eGF2(p) + eHF(i) - eHF(a) - eHF(b) + num = (2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b) + + SigC(p,q) = SigC(p,q) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + end do + end do + + Z(:) = 1d0/(1d0 - Z(:)) + +end subroutine self_energy_GF2 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 177f7e4..9647ef1 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -14,7 +14,7 @@ program QuAcK logical :: doCIS,doCIS_D,doCID,doCISD logical :: doRPA,doRPAx,doppRPA logical :: doADC - logical :: doG0F2,doevGF2,doG0F3,doevGF3 + logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical :: doG0W0,doevGW,doqsGW logical :: doG0T0,doevGT,doqsGT logical :: doMCMP2,doMinMCMP2 @@ -164,7 +164,8 @@ program QuAcK do_drCCD,do_rCCD,do_lCCD,do_pCCD, & doCIS,doCIS_D,doCID,doCISD, & doRPA,doRPAx,doppRPA, & - doG0F2,doevGF2,doG0F3,doevGF3, & + doG0F2,doevGF2,doqsGF2, & + doG0F3,doevGF3, & doG0W0,doevGW,doqsGW, & doG0T0,doevGT,doqsGT, & doMCMP2) @@ -839,6 +840,25 @@ program QuAcK end if +!------------------------------------------------------------------------ +! Perform qsGF2 calculation +!------------------------------------------------------------------------ + + if(doqsGF2) then + + call cpu_time(start_GF2) + + call qsGF2(maxSCF_GF,thresh_GF,n_diis_GF,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GF,nNuc,ZNuc,rNuc,ENuc, & + nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) + + call cpu_time(end_GF2) + + t_GF2 = end_GF2 - start_GF2 + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF2,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Compute G0F3 electronic binding energies !------------------------------------------------------------------------ diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index ab0d2d6..1401cc9 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -4,7 +4,8 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & do_drCCD,do_rCCD,do_lCCD,do_pCCD, & doCIS,doCIS_D,doCID,doCISD, & doRPA,doRPAx,doppRPA, & - doG0F2,doevGF2,doG0F3,doevGF3, & + doG0F2,doevGF2,doqsGF2, & + doG0F3,doevGF3, & doG0W0,doevGW,doqsGW, & doG0T0,doevGT,doqsGT, & doMCMP2) @@ -21,7 +22,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD logical,intent(out) :: doRPA,doRPAx,doppRPA - logical,intent(out) :: doG0F2,doevGF2,doG0F3,doevGF3 + logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW logical,intent(out) :: doG0T0,doevGT,doqsGT logical,intent(out) :: doMCMP2 @@ -66,6 +67,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & doG0F2 = .false. doevGF2 = .false. + doqsGF2 = .false. doG0F3 = .false. doevGF3 = .false. @@ -132,11 +134,12 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Read Green function methods read(1,*) - read(1,*) answer1,answer2,answer3,answer4 + read(1,*) answer1,answer2,answer3,answer4,answer5 if(answer1 == 'T') doG0F2 = .true. if(answer2 == 'T') doevGF2 = .true. - if(answer3 == 'T') doG0F3 = .true. - if(answer4 == 'T') doevGF3 = .true. + if(answer3 == 'T') doqsGF2 = .true. + if(answer4 == 'T') doG0F3 = .true. + if(answer5 == 'T') doevGF3 = .true. ! Read GW methods