From 2fc348ebe4d4335633fefc4792e5837ef3a2ebe9 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 4 Jul 2023 09:27:35 +0200 Subject: [PATCH] qsGTeh --- input/int | 14 -- input/methods | 2 +- input/options | 2 +- src/GT/GTeh_self_energy.f90 | 94 ++++++++++ src/GT/print_qsGTeh.f90 | 120 +++++++++++++ src/GT/qsGTeh.f90 | 341 ++++++++++++++++++++++++++++++++++++ src/QuAcK/QuAcK.f90 | 320 +++++++++++++++++---------------- 7 files changed, 713 insertions(+), 180 deletions(-) delete mode 100644 input/int create mode 100644 src/GT/GTeh_self_energy.f90 create mode 100644 src/GT/print_qsGTeh.f90 create mode 100644 src/GT/qsGTeh.f90 diff --git a/input/int b/input/int deleted file mode 100644 index 86aa017..0000000 --- a/input/int +++ /dev/null @@ -1,14 +0,0 @@ -# Debuggin mode? - F -# Chemist notation for two-electron integral? - T -# Exposant of the Slater geminal - 1.0 -# One-electron integrals: Ov Kin Nuc - T T T -# Two-electron integrals: ERI F12 Yuk Erf - T F F F -# Three-electron integrals: Type1 Type2 Type3 - F F F -# Four-electron integrals: Type1 Type2 Type3 - F F F diff --git a/input/methods b/input/methods index fc58a00..bb503e1 100644 --- a/input/methods +++ b/input/methods @@ -15,5 +15,5 @@ # G0W0* evGW* qsGW* SRG-qsGW ufG0W0 ufGW T F F F F F # G0T0pp evGTpp qsGTpp G0T0eh evGTeh qsGTeh - F F F T T F + F F F F F T # * unrestricted version available diff --git a/input/options b/input/options index 7d6351a..f5c5171 100644 --- a/input/options +++ b/input/options @@ -11,7 +11,7 @@ # 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 # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg - 10 0.00001 T 5 T 0.1 T F + 256 0.00001 T 5 T 0.0 T F # ACFDT: AC Kx XBS F T T # BSE: BSE dBSE dTDA evDyn ppBSE BSE2 diff --git a/src/GT/GTeh_self_energy.f90 b/src/GT/GTeh_self_energy.f90 new file mode 100644 index 0000000..af06e80 --- /dev/null +++ b/src/GT/GTeh_self_energy.f90 @@ -0,0 +1,94 @@ +subroutine GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rhoL,rhoR,EcGM,SigC) + +! Compute correlation part of the self-energy for GTeh + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + 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) :: e(nBas) + double precision,intent(in) :: Om(nS) + double precision,intent(in) :: rhoL(nBas,nBas,nS) + double precision,intent(in) :: rhoR(nBas,nBas,nS) + +! Local variables + + integer :: i,j,a,b + integer :: p,q,r + integer :: m + double precision :: eps + +! Output variables + + double precision,intent(out) :: EcGM + double precision,intent(out) :: SigC(nBas,nBas) + +! Initialize + + SigC(:,:) = 0d0 + +!----------------! +! GW self-energy ! +!----------------! + + ! Occupied part of the correlation self-energy + +!$OMP PARALLEL & +!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) & +!$OMP PRIVATE(m,i,q,p,eps) & +!$OMP DEFAULT(NONE) +!$OMP DO + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + do m=1,nS + do i=nC+1,nO + eps = e(p) - e(i) + Om(m) + SigC(p,q) = SigC(p,q) + rhoL(i,p,m)*rhoR(i,q,m)*eps/(eps**2 + eta**2) + end do + end do + end do + end do +!$OMP END DO +!$OMP END PARALLEL + + ! Virtual part of the correlation self-energy + +!$OMP PARALLEL & +!$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) & +!$OMP PRIVATE(m,a,q,p,eps) & +!$OMP DEFAULT(NONE) +!$OMP DO + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + do m=1,nS + do a=nO+1,nBas-nR + eps = e(p) - e(a) - Om(m) + SigC(p,q) = SigC(p,q) + rhoL(p,a,m)*rhoR(q,a,m)*eps/(eps**2 + eta**2) + end do + end do + end do + end do +!$OMP END DO +!$OMP END PARALLEL + + ! Galitskii-Migdal correlation energy + + EcGM = 0d0 + do m=1,nS + do a=nO+1,nBas-nR + do i=nC+1,nO + eps = e(a) - e(i) + Om(m) + EcGM = EcGM - rhoL(i,a,m)*rhoR(i,a,m)*eps/(eps**2 + eta**2) + end do + end do + end do + +end subroutine diff --git a/src/GT/print_qsGTeh.f90 b/src/GT/print_qsGTeh.f90 new file mode 100644 index 0000000..1bdf0b7 --- /dev/null +++ b/src/GT/print_qsGTeh.f90 @@ -0,0 +1,120 @@ +subroutine print_qsGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,dipole) + +! Print one-electron energies and other stuff for qsGTeh + + 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) :: ET + double precision,intent(in) :: EV + double precision,intent(in) :: EJ + double precision,intent(in) :: Ex + double precision,intent(in) :: EcGM + double precision,intent(in) :: EcRPA(nspin) + double precision,intent(in) :: Conv + double precision,intent(in) :: thresh + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: eGT(nBas) + double precision,intent(in) :: c(nBas) + double precision,intent(in) :: SigC(nBas,nBas) + double precision,intent(in) :: Z(nBas) + double precision,intent(in) :: dipole(ncart) + +! Local variables + + integer :: p,ixyz,HOMO,LUMO + double precision :: Gap + double precision,external :: trace_matrix + +! Output variables + + double precision,intent(out) :: EqsGT + +! HOMO and LUMO + + HOMO = nO + LUMO = HOMO + 1 + Gap = eGT(LUMO)-eGT(HOMO) + +! Compute energies + +! Dump results + + write(*,*)'-------------------------------------------------------------------------------' + if(nSCF < 10) then + write(*,'(1X,A21,I1,A3,I1,A12)')' Self-consistent qsG',nSCF,'Teh',nSCF,' calculation' + else + write(*,'(1X,A21,I2,A3,I2,A12)')' Self-consistent qsG',nSCF,'Teh',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_T (eV)','|','Z','|','e_QP (eV)','|' + write(*,*)'-------------------------------------------------------------------------------' + + do p=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)') & + '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' + enddo + + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A10,I3)') 'Iteration ',nSCF + write(*,'(2X,A14,F15.5)')'Convergence = ',Conv + write(*,*)'-------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') 'qsGTeh HOMO energy:',eGT(HOMO)*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'qsGTeh LUMO energy:',eGT(LUMO)*HaToeV,' eV' + write(*,'(2X,A30,F15.6,A3)') 'qsGTeh HOMO-LUMO gap :',Gap*HaToeV,' eV' + write(*,*)'-------------------------------------------' + write(*,'(2X,A30,F15.6,A3)') ' qsGTeh total energy:',ENuc + EqsGT,' au' + write(*,'(2X,A30,F15.6,A3)') ' qsGTeh exchange energy:',Ex,' au' + write(*,'(2X,A30,F15.6,A3)') ' GM@qsGTeh correlation energy:',EcGM,' au' + write(*,'(2X,A30,F15.6,A3)') 'ppRPA@qsGTeh correlation energy:',sum(EcRPA(:)),' 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: ',EcGM,' au' + write(*,'(A50)') '---------------------------------------' + write(*,'(A32,1X,F16.10,A3)') ' Electronic energy: ',EqsGT,' au' + write(*,'(A32,1X,F16.10,A3)') ' Nuclear repulsion: ',ENuc,' au' + write(*,'(A32,1X,F16.10,A3)') ' qsGTeh energy: ',ENuc + EqsGT,' 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)') ' qsGTeh MO coefficients' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,nBas,c) + write(*,*) + write(*,'(A50)') '---------------------------------------' + write(*,'(A32)') ' qsGTeh MO energies' + write(*,'(A50)') '---------------------------------------' + call matout(nBas,1,eGT) + write(*,*) + + endif + +end subroutine diff --git a/src/GT/qsGTeh.f90 b/src/GT/qsGTeh.f90 new file mode 100644 index 0000000..fbdf249 --- /dev/null +++ b/src/GT/qsGTeh.f90 @@ -0,0 +1,341 @@ +subroutine qsGTeh(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,BSE2,TDA_T,TDA, & + dBSE,dTDA,evDyn,singlet,triplet,eta,regularize,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 GTeh 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) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: BSE + logical,intent(in) :: BSE2 + logical,intent(in) :: TDA_T + 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 + logical,intent(in) :: regularize + + 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 + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: 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 :: ET + double precision :: EV + double precision :: EJ + double precision :: Ex + double precision :: EqsGT + double precision :: EcRPA + double precision :: EcBSE(nspin) + double precision :: EcAC(nspin) + double precision :: EcGM + double precision :: Conv + double precision :: rcond + double precision,external :: trace_matrix + double precision :: dipole(ncart) + + logical :: print_T = .true. + double precision,allocatable :: error_diis(:,:) + double precision,allocatable :: F_diis(:,:) + double precision,allocatable :: OmRPA(:) + double precision,allocatable :: XpY_RPA(:,:) + double precision,allocatable :: XmY_RPA(:,:) + double precision,allocatable :: rhoL_RPA(:,:,:) + double precision,allocatable :: rhoR_RPA(:,:,:) + double precision,allocatable :: c(:,:) + double precision,allocatable :: cp(:,:) + double precision,allocatable :: eGT(:) + double precision,allocatable :: eOld(:) + 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 qsGTeh calculation |' + write(*,*)'************************************************' + write(*,*) + +! Warning + + write(*,*) '!! ERIs in MO basis will be overwritten in qsGTeh !!' + write(*,*) + +! Stuff + + nBasSq = nBas*nBas + +! TDA for T + + if(TDA_T) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if + +! TDA + + if(TDA) then + write(*,*) 'Tamm-Dancoff approximation activated!' + write(*,*) + end if + +! Memory allocation + + allocate(eGT(nBas),eOld(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), & + OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rhoL_RPA(nBas,nBas,nS),rhoR_RPA(nBas,nBas,nS), & + error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + +! Initialization + + nSCF = -1 + n_diis = 0 + ispin = 2 + Conv = 1d0 + P(:,:) = PHF(:,:) + eGT(:) = eHF(:) + eOld(:) = eHF(:) + c(:,:) = cHF(:,:) + F_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + rcond = 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 linear response + + call linear_response(ispin,.false.,TDA_T,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) + if(print_T) call print_excitation('RPA@qsGTeh ',ispin,nS,OmRPA) + + ! Compute correlation part of the self-energy + + call GTeh_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA) + + if(regularize) then + +! call regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) +! call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z) + + else + + call GTeh_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,EcGM,SigC) + call GTeh_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGT,OmRPA,rhoL_RPA,rhoR_RPA,Z) + + endif + + ! 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) + + ! DIIS extrapolation + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + + end if + + ! Diagonalize Hamiltonian in AO basis + + Fp = matmul(transpose(X),matmul(F,X)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nBas,cp,eGT) + c = matmul(X,cp) + SigCp = matmul(transpose(c),matmul(SigCp,c)) + + ! Compute new density matrix in the AO basis + + P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + + ! Save quasiparticles energy for next cycle + + Conv = maxval(abs(error)) + eOld(:) = eGT(:) + + !------------------------------------------------------------------------ + ! Compute total energy + !------------------------------------------------------------------------ + + ! Kinetic energy + + ET = trace_matrix(nBas,matmul(P,T)) + + ! Potential energy + + EV = trace_matrix(nBas,matmul(P,V)) + + ! Coulomb energy + + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) + + ! Exchange energy + + Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) + + ! Total energy + + EqsGT = ET + EV + EJ + Ex + + ! Print results + + call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsGTeh(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigCp,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGT,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,OmRPA,XpY_RPA,XmY_RPA,rhoL_RPA,rhoR_RPA,error,error_diis,F_diis) + +! Perform BSE calculation + +! if(BSE) then + +! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & +! eGW,eGW,EcBSE) + +! if(exchange_kernel) then + +! EcBSE(1) = 0.5d0*EcBSE(1) +! EcBSE(2) = 1.5d0*EcBSE(2) + +! end if + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2) +! write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! Compute the BSE correlation energy via the adiabatic connection + +! if(doACFDT) then + +! write(*,*) '------------------------------------------------------' +! write(*,*) 'Adiabatic connection version of BSE correlation energy' +! write(*,*) '------------------------------------------------------' +! write(*,*) + +! if(doXBS) then + +! write(*,*) '*** scaled screening version (XBS) ***' +! write(*,*) + +! end if + +! call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) + +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy =',EcAC(1) + EcAC(2) +! write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcAC(1) + EcAC(2) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) + +! end if + +! end if + +end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index c74029e..07cdb93 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -40,19 +40,6 @@ program QuAcK logical :: exchange_kernel logical :: doXBS -! integer :: nShell -! integer,allocatable :: TotAngMomShell(:) -! integer,allocatable :: KShell(:) -! double precision,allocatable :: CenterShell(:,:) -! double precision,allocatable :: DShell(:,:) -! double precision,allocatable :: ExpShell(:,:) -! integer,allocatable :: max_ang_mom(:) -! double precision,allocatable :: min_exponent(:,:) -! double precision,allocatable :: max_exponent(:) - - integer :: TrialType - double precision,allocatable :: cTrial(:),gradient(:),hessian(:,:) - double precision,allocatable :: S(:,:) double precision,allocatable :: T(:,:) double precision,allocatable :: V(:,:) @@ -80,28 +67,14 @@ program QuAcK double precision :: start_HF ,end_HF ,t_HF double precision :: start_stab ,end_stab ,t_stab double precision :: start_KS ,end_KS ,t_KS - double precision :: start_MOM ,end_MOM ,t_MOM double precision :: start_AOtoMO ,end_AOtoMO ,t_AOtoMO - double precision :: start_CCD ,end_CCD ,t_CCD - double precision :: start_DCD ,end_DCD ,t_DCD - double precision :: start_CCSD ,end_CCSD ,t_CCSD - double precision :: start_CIS ,end_CIS ,t_CIS - double precision :: start_CID ,end_CID ,t_CID - double precision :: start_CISD ,end_CISD ,t_CISD - double precision :: start_FCI ,end_FCI ,t_FCI + double precision :: start_CC ,end_CC ,t_CC + double precision :: start_CI ,end_CI ,t_CI double precision :: start_RPA ,end_RPA ,t_RPA - double precision :: start_ADC ,end_ADC ,t_ADC - double precision :: start_GF2 ,end_GF2 ,t_GF2 - double precision :: start_GF3 ,end_GF3 ,t_GF3 - double precision :: start_G0W0 ,end_G0W0 ,t_G0W0 - double precision :: start_evGW ,end_evGW ,t_evGW - double precision :: start_qsGW ,end_qsGW ,t_qsGW - double precision :: start_ufGW ,end_ufGW ,t_ufGW - double precision :: start_G0T0 ,end_G0T0 ,t_G0T0 - double precision :: start_evGT ,end_evGT ,t_evGT - double precision :: start_qsGT ,end_qsGT ,t_qsGT - double precision :: start_MP2 ,end_MP2 ,t_MP2 - double precision :: start_MP3 ,end_MP3 ,t_MP3 + double precision :: start_GF ,end_GF ,t_GF + double precision :: start_GW ,end_GW ,t_GW + double precision :: start_GT ,end_GT ,t_GT + double precision :: start_MP ,end_MP ,t_MP integer :: maxSCF_HF,n_diis_HF double precision :: thresh_HF,level_shift @@ -136,10 +109,6 @@ program QuAcK logical :: BSE,dBSE,dTDA,evDyn,ppBSE,BSE2 - integer :: nMC,nEq,nWalk,nPrint,iSeed - double precision :: dt - logical :: doDrift - ! Hello World write(*,*) @@ -328,7 +297,7 @@ program QuAcK if(doMOM) then - call cpu_time(start_MOM) + call cpu_time(start_HF) if(unrestricted) then @@ -341,10 +310,10 @@ program QuAcK end if - call cpu_time(end_MOM) + call cpu_time(end_HF) - t_MOM = end_MOM - start_MOM - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM = ',t_MOM,' seconds' + t_HF = end_HF - start_HF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM = ',t_HF,' seconds' write(*,*) end if @@ -475,7 +444,7 @@ program QuAcK if(doMP2) then - call cpu_time(start_MP2) + call cpu_time(start_MP) if(unrestricted) then @@ -487,10 +456,10 @@ program QuAcK end if - call cpu_time(end_MP2) + call cpu_time(end_MP) - t_MP2 = end_MP2 - start_MP2 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP2,' seconds' + t_MP = end_MP - start_MP + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP2 = ',t_MP,' seconds' write(*,*) end if @@ -501,7 +470,7 @@ program QuAcK if(doMP3) then - call cpu_time(start_MP3) + call cpu_time(start_MP) if(unrestricted) then @@ -514,10 +483,10 @@ program QuAcK end if - call cpu_time(end_MP3) + call cpu_time(end_MP) - t_MP3 = end_MP3 - start_MP3 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP3 = ',t_MP3,' seconds' + t_MP = end_MP - start_MP + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MP3 = ',t_MP,' seconds' write(*,*) end if @@ -528,12 +497,12 @@ program QuAcK if(doCCD) then - call cpu_time(start_CCD) + call cpu_time(start_CC) call CCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) - call cpu_time(end_CCD) + call cpu_time(end_CC) - t_CCD = end_CCD - start_CCD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds' + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CC,' seconds' write(*,*) end if @@ -544,13 +513,13 @@ program QuAcK if(doDCD) then - call cpu_time(start_DCD) + call cpu_time(start_CC) call DCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & ERI_MO,ENuc,ERHF,eHF) - call cpu_time(end_DCD) + call cpu_time(end_CC) - t_DCD = end_DCD - start_DCD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_DCD,' seconds' + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for DCD = ',t_CC,' seconds' write(*,*) end if @@ -563,12 +532,12 @@ program QuAcK if(doCCSD) then - call cpu_time(start_CCSD) + call cpu_time(start_CC) call CCSD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) - call cpu_time(end_CCSD) + call cpu_time(end_CC) - t_CCSD = end_CCSD - start_CCSD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CCSD,' seconds' + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCSD or CCSD(T)= ',t_CC,' seconds' write(*,*) end if @@ -579,12 +548,12 @@ program QuAcK if(do_drCCD) then - call cpu_time(start_CCD) + call cpu_time(start_CC) call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) - call cpu_time(end_CCD) + call cpu_time(end_CC) - t_CCD = end_CCD - start_CCD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CCD,' seconds' + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for direct ring CCD = ',t_CC,' seconds' write(*,*) end if @@ -595,12 +564,12 @@ program QuAcK if(do_rCCD) then - call cpu_time(start_CCD) + call cpu_time(start_CC) call rCCD(.false.,maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,eHF) - call cpu_time(end_CCD) + call cpu_time(end_CC) - t_CCD = end_CCD - start_CCD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CCD,' seconds' + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for rCCD = ',t_CC,' seconds' write(*,*) end if @@ -611,13 +580,12 @@ program QuAcK if(do_crCCD) then - call cpu_time(start_CCD) + call cpu_time(start_CC) call crCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) + call cpu_time(end_CC) - call cpu_time(end_CCD) - - t_CCD = end_CCD - start_CCD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CCD,' seconds' + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for crossed-ring CCD = ',t_CC,' seconds' write(*,*) end if @@ -628,13 +596,13 @@ program QuAcK if(do_lCCD) then - call cpu_time(start_CCD) + call cpu_time(start_CC) call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR, & ERI_MO,ENuc,ERHF,eHF) - call cpu_time(end_CCD) + call cpu_time(end_CC) - t_CCD = end_CCD - start_CCD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CCD,' seconds' + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ladder CCD = ',t_CC,' seconds' write(*,*) end if @@ -645,12 +613,12 @@ program QuAcK if(dopCCD) then - call cpu_time(start_CCD) + call cpu_time(start_CC) call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) - call cpu_time(end_CCD) + call cpu_time(end_CC) - t_CCD = end_CCD - start_CCD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CCD,' seconds' + t_CC = end_CC - start_CC + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for pair CCD = ',t_CC,' seconds' write(*,*) end if @@ -661,7 +629,7 @@ program QuAcK if(doCIS) then - call cpu_time(start_CIS) + call cpu_time(start_CI) if(unrestricted) then @@ -674,10 +642,10 @@ program QuAcK end if - call cpu_time(end_CIS) + call cpu_time(end_CI) - t_CIS = end_CIS - start_CIS - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CIS = ',t_CIS,' seconds' + t_CI = end_CI - start_CI + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CIS = ',t_CI,' seconds' write(*,*) end if @@ -688,12 +656,12 @@ program QuAcK if(doCID) then - call cpu_time(start_CID) + call cpu_time(start_CI) call CID(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,F_MO,ERHF) - call cpu_time(end_CID) + call cpu_time(end_CI) - t_CID = end_CID - start_CID - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CID,' seconds' + t_CI = end_CI - start_CI + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CID = ',t_CI,' seconds' write(*,*) end if @@ -704,12 +672,12 @@ program QuAcK if(doCISD) then - call cpu_time(start_CISD) + call cpu_time(start_CI) call CISD(singlet,triplet,nBas,nC,nO,nV,nR,ERI_MO,F_MO,ERHF) - call cpu_time(end_CISD) + call cpu_time(end_CI) - t_CISD = end_CISD - start_CISD - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CISD,' seconds' + t_CI = end_CI - start_CI + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CISD = ',t_CI,' seconds' write(*,*) end if @@ -829,7 +797,7 @@ program QuAcK if(doG0F2) then - call cpu_time(start_GF2) + call cpu_time(start_GF) if(unrestricted) then @@ -844,10 +812,10 @@ program QuAcK end if - call cpu_time(end_GF2) + call cpu_time(end_GF) - t_GF2 = end_GF2 - start_GF2 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF2,' seconds' + t_GF = end_GF - start_GF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' write(*,*) end if @@ -858,7 +826,7 @@ program QuAcK if(doevGF2) then - call cpu_time(start_GF2) + call cpu_time(start_GF) if(unrestricted) then @@ -874,10 +842,10 @@ program QuAcK end if - call cpu_time(end_GF2) + call cpu_time(end_GF) - t_GF2 = end_GF2 - start_GF2 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF2,' seconds' + t_GF = end_GF - start_GF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF2 = ',t_GF,' seconds' write(*,*) end if @@ -888,7 +856,7 @@ program QuAcK if(doqsGF2) then - call cpu_time(start_GF2) + call cpu_time(start_GF) if(unrestricted) then @@ -903,10 +871,10 @@ program QuAcK end if - call cpu_time(end_GF2) + call cpu_time(end_GF) - t_GF2 = end_GF2 - start_GF2 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF2,' seconds' + t_GF = end_GF - start_GF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGF2 = ',t_GF,' seconds' write(*,*) end if @@ -917,7 +885,7 @@ program QuAcK if(doG0F3) then - call cpu_time(start_GF3) + call cpu_time(start_GF) if(unrestricted) then @@ -929,10 +897,10 @@ program QuAcK end if - call cpu_time(end_GF3) + call cpu_time(end_GF) - t_GF3 = end_GF3 - start_GF3 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF3,' seconds' + t_GF = end_GF - start_GF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds' write(*,*) end if @@ -943,7 +911,7 @@ program QuAcK if(doevGF3) then - call cpu_time(start_GF3) + call cpu_time(start_GF) if(unrestricted) then @@ -955,10 +923,10 @@ program QuAcK end if - call cpu_time(end_GF3) + call cpu_time(end_GF) - t_GF3 = end_GF3 - start_GF3 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF3,' seconds' + t_GF = end_GF - start_GF + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for GF3 = ',t_GF,' seconds' write(*,*) end if @@ -971,7 +939,7 @@ program QuAcK if(doG0W0) then - call cpu_time(start_G0W0) + call cpu_time(start_GW) if(unrestricted) then call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip, & @@ -987,16 +955,14 @@ program QuAcK else call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,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) -! call ehTM(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,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 end if - call cpu_time(end_G0W0) + call cpu_time(end_GW) - t_G0W0 = end_G0W0 - start_G0W0 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_G0W0,' seconds' + t_GW = end_GW - start_GW + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_GW,' seconds' write(*,*) end if @@ -1007,7 +973,7 @@ program QuAcK if(doevGW) then - call cpu_time(start_evGW) + call cpu_time(start_GW) if(unrestricted) then call evUGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & @@ -1021,10 +987,10 @@ program QuAcK 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) + call cpu_time(end_GW) - t_evGW = end_evGW - start_evGW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_evGW,' seconds' + t_GW = end_GW - start_GW + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGW = ',t_GW,' seconds' write(*,*) end if @@ -1035,7 +1001,7 @@ program QuAcK if(doqsGW) then - call wall_time(start_qsGW) + call wall_time(start_GW) if(unrestricted) then @@ -1052,10 +1018,10 @@ program QuAcK end if - call wall_time(end_qsGW) + call wall_time(end_GW) - t_qsGW = end_qsGW - start_qsGW - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_qsGW,' seconds' + t_GW = end_GW - start_GW + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds' write(*,*) end if @@ -1066,7 +1032,7 @@ program QuAcK if(doSRGqsGW) then - call wall_time(start_qsGW) + call wall_time(start_GW) if(unrestricted) then @@ -1080,10 +1046,10 @@ program QuAcK end if - call wall_time(end_qsGW) + call wall_time(end_GW) - t_qsGW = end_qsGW - start_qsGW - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_qsGW,' seconds' + t_GW = end_GW - start_GW + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GW,' seconds' write(*,*) end if @@ -1094,12 +1060,12 @@ program QuAcK if(doufG0W0) then - call cpu_time(start_ufGW) + call cpu_time(start_GW) call ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,TDA_W) - call cpu_time(end_ufGW) + call cpu_time(end_GW) - t_ufGW = end_ufGW - start_ufGW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_ufGW,' seconds' + t_GW = end_GW - start_GW + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufG0W0 = ',t_GW,' seconds' write(*,*) if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) @@ -1112,13 +1078,13 @@ program QuAcK if(doufGW) then - call cpu_time(start_ufGW) + call cpu_time(start_GW) call ufGW(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) ! call CCGW(maxSCF_CC,thresh_CC,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) - call cpu_time(end_ufGW) + call cpu_time(end_GW) - t_ufGW = end_ufGW - start_ufGW - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_ufGW,' seconds' + t_GW = end_GW - start_GW + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for ufGW = ',t_GW,' seconds' write(*,*) if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eG0W0) @@ -1133,7 +1099,7 @@ program QuAcK if(doG0T0pp) then - call cpu_time(start_G0T0) + call cpu_time(start_GT) if(unrestricted) then @@ -1152,10 +1118,10 @@ program QuAcK end if - call cpu_time(end_G0T0) + call cpu_time(end_GT) - t_G0T0 = end_G0T0 - start_G0T0 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds' + t_GT = end_GT - start_GT + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_GT,' seconds' write(*,*) end if @@ -1166,7 +1132,7 @@ program QuAcK if(doevGTpp) then - call cpu_time(start_evGT) + call cpu_time(start_GT) if(unrestricted) then @@ -1185,10 +1151,10 @@ program QuAcK end if - call cpu_time(end_evGT) + call cpu_time(end_GT) - t_evGT = end_evGT - start_evGT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_evGT,' seconds' + t_GT = end_GT - start_GT + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_GT,' seconds' write(*,*) end if @@ -1199,7 +1165,7 @@ program QuAcK if(doqsGTpp) then - call cpu_time(start_qsGT) + call cpu_time(start_GT) if(unrestricted) then @@ -1216,10 +1182,10 @@ program QuAcK end if - call cpu_time(end_qsGT) + call cpu_time(end_GT) - t_qsGT = end_qsGT - start_qsGT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGT = ',t_qsGT,' seconds' + t_GT = end_GT - start_GT + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for qsGT = ',t_GT,' seconds' write(*,*) end if @@ -1232,7 +1198,7 @@ program QuAcK if(doG0T0eh) then - call cpu_time(start_G0T0) + call cpu_time(start_GT) if(unrestricted) then @@ -1245,10 +1211,10 @@ program QuAcK end if - call cpu_time(end_G0T0) + call cpu_time(end_GT) - t_G0T0 = end_G0T0 - start_G0T0 - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds' + t_GT = end_GT - start_GT + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_GT,' seconds' write(*,*) end if @@ -1259,7 +1225,7 @@ program QuAcK if(doevGTeh) then - call cpu_time(start_evGT) + call cpu_time(start_GT) if(unrestricted) then else @@ -1268,10 +1234,36 @@ program QuAcK BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,ppBSE,singlet,triplet,linGT,eta_GT,regGT, & nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO,PHF,cHF,eHF,Vxc,eG0T0) end if - call cpu_time(end_evGT) + call cpu_time(end_GT) - t_evGT = end_evGT - start_evGT - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_evGT,' seconds' + t_GT = end_GT - start_GT + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for evGT = ',t_GT,' seconds' + write(*,*) + + end if + +!------------------------------------------------------------------------ +! Perform qsGTeh calculation +!------------------------------------------------------------------------ + + if(doqsGTeh) then + + call wall_time(start_GT) + + if(unrestricted) then + + else + + call qsGTeh(maxSCF_GT,thresh_GT,n_diis_GT,doACFDT,exchange_kernel,doXBS, & + BSE,BSE2,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta_GT,regGT,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) + + end if + + call wall_time(end_GT) + + t_GT = end_GT - start_GT + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for qsGW = ',t_GT,' seconds' write(*,*) end if @@ -1282,13 +1274,13 @@ program QuAcK if(doFCI) then - call cpu_time(start_FCI) + call cpu_time(start_CI) write(*,*) ' FCI is not yet implemented! Sorry.' ! call FCI(nBas,nC,nO,nV,nR,ERI_MO,eHF) - call cpu_time(end_FCI) + call cpu_time(end_CI) - t_FCI = end_FCI - start_FCI - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_FCI,' seconds' + t_CI = end_CI - start_CI + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_CI,' seconds' write(*,*) end if