From d051302107116ca93d1bdc99ec43cdd2e0c4ac63 Mon Sep 17 00:00:00 2001 From: Loris Burth Date: Tue, 15 Apr 2025 10:20:26 +0200 Subject: [PATCH] complex linearized evGW@cRHF --- src/GW/complex_RGW.f90 | 22 +- src/GW/complex_cRG0W0.f90 | 9 +- src/GW/complex_evRGW.f90 | 292 +++++++++++++++++++++++ src/GW/print_complex_evRGW.f90 | 67 ++++++ src/QuAcK/QuAcK.f90 | 3 +- src/QuAcK/RQuAcK.f90 | 2 +- src/utils/complex_DIIS_extrapolation.f90 | 2 +- 7 files changed, 385 insertions(+), 12 deletions(-) create mode 100644 src/GW/complex_evRGW.f90 create mode 100644 src/GW/print_complex_evRGW.f90 diff --git a/src/GW/complex_RGW.f90 b/src/GW/complex_RGW.f90 index 4e18d8a..40b1225 100644 --- a/src/GW/complex_RGW.f90 +++ b/src/GW/complex_RGW.f90 @@ -1,4 +1,4 @@ -subroutine complex_RGW(dotest,docG0W0,maxSCF,thresh,max_diis,doACFDT, & +subroutine complex_RGW(dotest,docG0W0,doevGW,maxSCF,thresh,max_diis,doACFDT, & exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & linearize,eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF, & S,X,T,V,Hc,ERI_AO,ERI_MO,CAP_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) @@ -12,7 +12,7 @@ subroutine complex_RGW(dotest,docG0W0,maxSCF,thresh,max_diis,doACFDT, logical,intent(in) :: dotest - logical,intent(in) :: docG0W0 + logical,intent(in) :: docG0W0,doevGW integer,intent(in) :: maxSCF integer,intent(in) :: max_diis @@ -82,4 +82,22 @@ subroutine complex_RGW(dotest,docG0W0,maxSCF,thresh,max_diis,doACFDT, write(*,*) end if + +!------------------------------------------------------------------------ +! Perform evGW calculation +!------------------------------------------------------------------------ + + if(doevGW) then + + call wall_time(start_GW) + call complex_evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & + singlet,triplet,linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + call wall_time(end_GW) + + t_GW = end_GW - start_GW + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for evGW = ',t_GW,' seconds' + write(*,*) + + end if + end subroutine diff --git a/src/GW/complex_cRG0W0.f90 b/src/GW/complex_cRG0W0.f90 index 54aaa6f..e625f5f 100644 --- a/src/GW/complex_cRG0W0.f90 +++ b/src/GW/complex_cRG0W0.f90 @@ -59,8 +59,6 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, double precision,allocatable :: Re_Z(:) double precision,allocatable :: Im_Z(:) complex*16,allocatable :: Om(:) - double precision,allocatable :: Re_Om(:) - double precision,allocatable :: Im_Om(:) complex*16,allocatable :: XpY(:,:) complex*16,allocatable :: XmY(:,:) complex*16,allocatable :: rho(:,:,:) @@ -106,7 +104,7 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, ! Memory allocation allocate(Aph(nS,nS),Bph(nS,nS),Re_SigC(nOrb),Im_SigC(nOrb),Re_Z(nOrb),Im_Z(nOrb),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS), & - Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),Re_eHF(nOrb),Im_eHF(nOrb),Re_Om(nS),Im_Om(nS)) + Re_eGW(nOrb),Im_eGW(nOrb),Re_eGWlin(nOrb),Im_eGWlin(nOrb),Re_eHF(nOrb),Im_eHF(nOrb),) Re_eHF(:) = real(eHF(:)) Im_eHF(:) = aimag(eHF(:)) !-------------------! @@ -117,8 +115,6 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, if(.not.TDA_W) call complex_phRLR_B(isp_W,dRPA_W,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) call complex_phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) - Re_Om(:) = real(Om(:)) - Im_Om(:) = aimag(Om(:)) !if(print_W) call print_excitation_energies('phRPA@RHF','singlet',nS,Om) @@ -152,7 +148,6 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, write(*,*) ' *** Quasiparticle energies obtained by root search *** ' write(*,*) - write(*,*) "ROOT SEARCH NOT IMPLEMENTED YET" call complex_RGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS, & Re_eHF,Im_eHF,Om,rho,Re_eGWlin,Im_eGWlin,Re_eHF,Im_eHF, & Re_eGW,Im_eGW,Re_Z,Im_Z) @@ -160,7 +155,7 @@ subroutine complex_cRG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2, ! Plot self-energy, renormalization factor, and spectral function ! -! if(plot_self) call RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) + if(plot_self) call complex_RGW_plot_self_energy(nOrb,eta,nC,nO,nV,nR,nS,eHF,eHF,Om,rho) ! !! Cumulant expansion ! diff --git a/src/GW/complex_evRGW.f90 b/src/GW/complex_evRGW.f90 new file mode 100644 index 0000000..8405553 --- /dev/null +++ b/src/GW/complex_evRGW.f90 @@ -0,0 +1,292 @@ +subroutine complex_evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & + singlet,triplet,linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + +! Perform self-consistent eigenvalue-only GW calculation + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dotest + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + logical,intent(in) :: doACFDT + logical,intent(in) :: exchange_kernel + logical,intent(in) :: doXBS + logical,intent(in) :: dophBSE + logical,intent(in) :: dophBSE2 + logical,intent(in) :: TDA_W + logical,intent(in) :: TDA + logical,intent(in) :: dBSE + logical,intent(in) :: dTDA + logical,intent(in) :: doppBSE + logical,intent(in) :: singlet + logical,intent(in) :: triplet + logical,intent(in) :: linearize + double precision,intent(in) :: eta + logical,intent(in) :: doSRG + + integer,intent(in) :: nBas + integer,intent(in) :: nOrb + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + complex*16,intent(in) :: eHF(nOrb) + complex*16,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + complex*16,intent(in) :: dipole_int(nOrb,nOrb,ncart) + +! Local variables + + logical :: dRPA = .true. + integer :: ispin + integer :: nSCF + integer :: n_diis + double precision :: flow + double precision :: rcond + double precision :: Conv + complex*16 :: EcRPA + complex*16 :: EcBSE(nspin) + complex*16 :: EcGM + double precision :: alpha + complex*16,allocatable :: Aph(:,:) + complex*16,allocatable :: Bph(:,:) + complex*16,allocatable :: error_diis(:,:) + complex*16,allocatable :: e_diis(:,:) + complex*16,allocatable :: eGW(:) + complex*16,allocatable :: eOld(:) + double precision,allocatable :: Re_eGW(:) + double precision,allocatable :: Im_eGW(:) + double precision,allocatable :: Re_eOld(:) + double precision,allocatable :: Im_eOld(:) + double precision,allocatable :: Re_eHF(:) + double precision,allocatable :: Im_eHF(:) + double precision,allocatable :: Re_Z(:) + double precision,allocatable :: Im_Z(:) + double precision,allocatable :: Re_SigC(:) + double precision,allocatable :: Im_SigC(:) + complex*16,allocatable :: Om(:) + complex*16,allocatable :: XpY(:,:) + complex*16,allocatable :: XmY(:,:) + complex*16,allocatable :: rho(:,:,:) + +! Hello world + + write(*,*) + write(*,*)'*******************************' + write(*,*)'* Restricted evGW Calculation *' + write(*,*)'*******************************' + write(*,*) + +! TDA for W + + if(TDA_W) then + write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' + write(*,*) + end if + +! SRG regularization + + flow = 100d0 + + if(doSRG) then + + write(*,*) '*** SRG regularized evGW scheme ***' + write(*,*) + + end if + +! Memory allocation + + allocate(Aph(nS,nS),Bph(nS,nS),Re_eGW(nOrb),Im_eGW(nOrb),eGW(nOrb),Re_eOld(nOrb),Im_eOld(nOrb),& + eOld(nOrb),Re_Z(nOrb),Im_Z(nOrb),Re_SigC(nOrb),Im_SigC(nOrb), & + Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nOrb,nOrb,nS),error_diis(nOrb,max_diis),e_diis(nOrb,max_diis),& + Re_eHF(nOrb),Im_eHF(nOrb)) + +! Initialization + + nSCF = 0 + ispin = 1 + n_diis = 0 + Conv = 1d0 + e_diis(:,:) = 0d0 + error_diis(:,:) = 0d0 + Re_eHF(:) = real(eHF(:)) + Im_eHF(:) = aimag(eHF(:)) + eGW(:) = eHF(:) + Re_eGW(:) = Re_eHF(:) + Im_eGW(:) = Im_eHF(:) + eOld(:) = eGW(:) + Re_eOld(:) = Re_eGW(:) + Im_eOld(:) = Im_eGW(:) + Re_Z(:) = 1d0 + Im_Z(:) = 0d0 + rcond = 0d0 + +!------------------------------------------------------------------------ +! Main loop +!------------------------------------------------------------------------ + + do while(Conv > thresh .and. nSCF <= maxSCF) + + ! Compute screening + + call complex_phRLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI,Aph) + if(.not.TDA_W) call complex_phRLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI,Bph) + + call complex_phRLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + + ! Compute spectral weights + + call complex_RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI,XpY,rho) + + ! Compute correlation part of the self-energy + + call complex_RGW_self_energy_diag(eta,nBas,nOrb,nC,nO,nV,nR,nS,Re_eGW,Im_eGW,Om,rho,& + EcGM,Re_SigC,Im_SigC,Re_Z,Im_Z) + + ! Solve the quasi-particle equation + + if(linearize) then + + write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' + write(*,*) + + Re_eGW(:) = Re_eHF(:) + Re_SigC(:) + Im_eGW(:) = Im_eHF(:) + Im_SigC(:) + eGW = cmplx(Re_eGW + Im_eGW,kind=8) + else + + write(*,*) ' *** Quasiparticle energies obtained by root search *** ' + write(*,*) + + call complex_RGW_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,nS,Re_eHF,Im_eHF,Om,& + rho,Re_eOld,Im_eOld,Re_eGW,Im_eGW,Re_Z,Im_Z) + + eGW = cmplx(Re_eGW + Im_eGW,kind=8) + end if + + ! Convergence criteria + + Conv = maxval(abs(eGW - eOld)) + + ! Print results + + call print_complex_evRGW(nBas,nO,nSCF,Conv,Re_eHF,Im_eHF,ENuc,ERHF,Re_SigC,Im_SigC,& + Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM) + ! DIIS extrapolation + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call complex_DIIS_extrapolation(rcond,nOrb,nOrb,n_diis,error_diis,e_diis,eGW-eOld,eGW) + + end if + + ! Save quasiparticles energy for next cycle + Re_eGW(:) = real(eGW(:)) + Im_eGW(:) = aimag(eGW(:)) + eOld(:) = eGW(:) + Re_eOld(:) = real(eOld(:)) + Im_eOld(:) = aimag(eOld(:)) + + ! Increment + + nSCF = nSCF + 1 + + end do +!------------------------------------------------------------------------ +! End main loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF+1) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + end if + +!--------------------! +! Cumulant expansion ! +!--------------------! +! +!! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z) +! +!! Deallocate memory +! +! deallocate(Aph,Bph,eOld,Z,SigC,Om,XpY,XmY,rho,error_diis,e_diis) +! +!! Perform BSE calculation +! +! if(dophBSE) then +! +! call RGW_phBSE(dophBSE2,exchange_kernel,TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta, & +! nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eGW,eGW,EcBSE) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@evGW@RHF correlation energy (singlet) = ',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@evGW@RHF correlation energy (triplet) = ',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@evGW@RHF correlation energy = ',sum(EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@BSE@evGW@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +!! Compute the BSE correlation energy via the adiabatic connection +! +! if(doACFDT) then +! +! call RGW_phACFDT(exchange_kernel,doXBS,TDA_W,TDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,eGW,eGW,EcBSE) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@evGW@RHF correlation energy (singlet) = ',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@evGW@RHF correlation energy (triplet) = ',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@evGW@RHF correlation energy = ',sum(EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@evGW@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +! end if +! +! end if +! +! if(doppBSE) then +! +! call RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGW,EcBSE) +! +! write(*,*) +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGW@RHF correlation energy (singlet) = ',EcBSE(1),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGW@RHF correlation energy (triplet) = ',EcBSE(2),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGW@RHF correlation energy = ',sum(EcBSE),' au' +! write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@evGW@RHF total energy = ',ENuc + ERHF + sum(EcBSE),' au' +! write(*,*)'-------------------------------------------------------------------------------' +! write(*,*) +! +! end if +! +!! Testing zone +! +! if(dotest) then +! +! call dump_test_value('R','evGW correlation energy',EcRPA) +! call dump_test_value('R','evGW HOMO energy',eGW(nO)) +! call dump_test_value('R','evGW LUMO energy',eGW(nO+1)) +! +! end if +! +end subroutine diff --git a/src/GW/print_complex_evRGW.f90 b/src/GW/print_complex_evRGW.f90 new file mode 100644 index 0000000..e7f0c74 --- /dev/null +++ b/src/GW/print_complex_evRGW.f90 @@ -0,0 +1,67 @@ +subroutine print_complex_evRGW(nBas,nO,nSCF,Conv,Re_eHF,Im_eHF,ENuc,ERHF,Re_SigC,Im_SigC,Re_Z,Im_Z,Re_eGW,Im_eGW,EcRPA,EcGM) + +! Print one-electron energies and other stuff for G0W0 + + implicit none + include 'parameters.h' + + integer,intent(in) :: nBas,nO,nSCF + double precision,intent(in) :: Conv + double precision,intent(in) :: ENuc + complex*16,intent(in) :: ERHF + complex*16,intent(in) :: EcRPA + complex*16,intent(in) :: EcGM + double precision,intent(in) :: Re_eHF(nBas) + double precision,intent(in) :: Im_eHF(nBas) + double precision,intent(in) :: Re_SigC(nBas) + double precision,intent(in) :: Im_SigC(nBas) + double precision,intent(in) :: Re_Z(nBas) + double precision,intent(in) :: Im_Z(nBas) + double precision,intent(in) :: Re_eGW(nBas) + double precision,intent(in) :: Im_eGW(nBas) + + integer :: p,index_homo,index_lumo + double precision :: Re_eHOMO,Re_eLUMO,Im_eHOMO,Im_eLUMO,Re_Gap,Im_Gap + +! HOMO and LUMO + + index_homo = maxloc(Re_eGW(1:nO),1) + Re_eHOMO = Re_eGW(index_homo) + Im_eHOMO = Im_eGW(index_homo) + index_lumo = minloc(Re_eGW(nO+1:nBas),1) + nO + Re_eLUMO = Re_eGW(index_lumo) + Im_eLUMO = Im_eGW(index_lumo) + Re_Gap = Re_eLUMO-Re_eHOMO + Im_Gap = Im_eLUMO-Im_eHOMO + +! Dump results + write(*,*)'-------------------------------------------------------------------------------' + if(nSCF < 10) then + write(*,'(1X,A20,I1,A1,I1,A17)')' Self-consistent evG',nSCF,'W',nSCF,'@cRHF calculation' + elseif(nSCF < 100) then + write(*,'(1X,A20,I2,A1,I2,A17)')' Self-consistent evG',nSCF,'W',nSCF,'@cRHF calculation' + else + write(*,'(1X,A20,I3,A1,I3,A17)')' Self-consistent evG',nSCF,'W',nSCF,'@cRHF calculation' + end if + 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)') & + '|','#','|','Re(e_HF) (eV)','|','Re(Sig_GW) (eV)','|','Re(Z)','|','Re(e_GW) (eV)','|' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & + '|','#','|','Im(e_HF) (eV)','|','Im(Sig_GW) (eV)','|','Im(Z)','|','Im(e_GW) (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,'|',Re_eHF(p)*HaToeV,'|',Re_SigC(p)*HaToeV,'|',Re_Z(p),'|',Re_eGW(p)*HaToeV,'|' + 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,'|',Im_eHF(p)*HaToeV,'|',Im_SigC(p)*HaToeV,'|',Im_Z(p),'|',Im_eGW(p)*HaToeV,'|' + if(p==nO) then + write(*,*)'-------------------------------------------------------------------------------' + end if + write(*,*)'-------------------------------------------------------------------------------' + end do + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A10,I3)') 'Iteration ',nSCF + write(*,'(2X,A14,F15.5)')'Convergence = ',Conv + write(*,*) +end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 8d99d4d..c93a27f 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -131,7 +131,8 @@ program QuAcK docG0W0,docG0F2, & doRtest,doUtest,doGtest) doCAP = docG0W0 .or. docG0F2 .or. docRHF ! Add different cases if they need CAP - + docG0W0 = docG0W0 .or. (doG0W0 .and. docRHF) + docG0F2 = docG0F2 .or. (doG0F2 .and. docRHF) !--------------------------! ! Read options for methods ! !--------------------------! diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 85a8333..092dbce 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -418,7 +418,7 @@ doGF = doG0F2 .or. doevGF2 .or. doqsGF2 .or. doufG0F02 .or. doG0F3 .or. doevGF3 if(doGW .and. docRHF) then call wall_time(start_GW) - call complex_RGW(dotest,docG0W0,maxSCF_GW,thresh_GW,max_diis_GW, & + call complex_RGW(dotest,docG0W0,doevGW,maxSCF_GW,thresh_GW,max_diis_GW, & doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,doppBSE,TDA_W,TDA,dBSE,dTDA,singlet,triplet, & lin_GW,eta_GW,reg_GW,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T, & V,Hc,ERI_AO,complex_ERI_MO,complex_CAP_MO,dipole_int_AO,complex_dipole_int_MO,complex_PHF,complex_cHF,complex_eHF) diff --git a/src/utils/complex_DIIS_extrapolation.f90 b/src/utils/complex_DIIS_extrapolation.f90 index 37ccc9b..cc039c7 100644 --- a/src/utils/complex_DIIS_extrapolation.f90 +++ b/src/utils/complex_DIIS_extrapolation.f90 @@ -49,7 +49,7 @@ subroutine complex_DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_ call complex_linear_solve(n_diis+1,A,b,w,rcond) ! Check condition number if (rcond < 1.0d-10) then - write(*,*) "!!! DIIS system ill-conditioned: rcond = ", rcond, " - Skipping DIIS update !!!" + !write(*,*) "!!! DIIS system ill-conditioned: rcond = ", rcond, " - Skipping DIIS update !!!" else ! Perform extrapolation only if the system is well-conditioned e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis)))