From eb6f35a7993f221c2ddd15a08d082bcc4ac4be75 Mon Sep 17 00:00:00 2001 From: pfloos Date: Wed, 11 Sep 2024 11:39:03 +0200 Subject: [PATCH] starting cleaning SRG --- src/GW/GW_regularization.f90 | 2 +- src/GW/SRG_qsRGW.f90 | 67 ++++++++++++++++++------------------ src/GW/SRG_qsUGW.f90 | 2 +- src/GW/SRG_self_energy.f90 | 30 +++++++++------- src/GW/USRG_self_energy.f90 | 32 +++++++++-------- src/GW/qsRGW.f90 | 2 -- 6 files changed, 71 insertions(+), 64 deletions(-) diff --git a/src/GW/GW_regularization.f90 b/src/GW/GW_regularization.f90 index ed77ae4..33dead2 100644 --- a/src/GW/GW_regularization.f90 +++ b/src/GW/GW_regularization.f90 @@ -28,7 +28,7 @@ subroutine GW_regularization(nBas,nC,nO,nV,nR,nS,e,Om,rho) ! SRG flow parameter - s = 100d0 + s = 500d0 ! Regularize excitation densities diff --git a/src/GW/SRG_qsRGW.f90 b/src/GW/SRG_qsRGW.f90 index ef52143..42ee7b7 100644 --- a/src/GW/SRG_qsRGW.f90 +++ b/src/GW/SRG_qsRGW.f90 @@ -76,8 +76,8 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS double precision :: dipole(ncart) logical :: dRPA_W = .true. - logical :: print_W = .true. - double precision,allocatable :: error_diis(:,:) + logical :: print_W = .false. + double precision,allocatable :: err_diis(:,:) double precision,allocatable :: F_diis(:,:) double precision,allocatable :: Aph(:,:) double precision,allocatable :: Bph(:,:) @@ -97,9 +97,7 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS double precision,allocatable :: SigC(:,:) double precision,allocatable :: SigCp(:,:) double precision,allocatable :: Z(:) - double precision,allocatable :: error(:,:) - - double precision,parameter :: flow = 500d0 + double precision,allocatable :: err(:,:) ! Hello world @@ -141,7 +139,7 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS allocate(F(nBas,nBas)) allocate(J(nBas,nBas)) allocate(K(nBas,nBas)) - allocate(error(nBas,nBas)) + allocate(err(nBas,nBas)) allocate(SigCp(nBas,nBas)) allocate(Aph(nS,nS)) @@ -151,21 +149,21 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS allocate(XmY(nS,nS)) allocate(rho(nOrb,nOrb,nS)) - allocate(error_diis(nBas_Sq,max_diis)) + allocate(err_diis(nBas_Sq,max_diis)) allocate(F_diis(nBas_Sq,max_diis)) ! Initialization - nSCF = -1 - n_diis = 0 - ispin = 1 - Conv = 1d0 - P(:,:) = PHF(:,:) - eGW(:) = eHF(:) - eOld(:) = eHF(:) - c(:,:) = cHF(:,:) - F_diis(:,:) = 0d0 - error_diis(:,:) = 0d0 + nSCF = -1 + n_diis = 0 + ispin = 1 + Conv = 1d0 + P(:,:) = PHF(:,:) + eGW(:) = eHF(:) + eOld(:) = eHF(:) + c(:,:) = cHF(:,:) + F_diis(:,:) = 0d0 + err_diis(:,:) = 0d0 rcond = 0d0 !------------------------------------------------------------------------ @@ -215,7 +213,7 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS tlr = tlr + tlr2 -tlr1 - if(print_W) call print_excitation_energies('phRPA@RGW','singlet',nS,Om) + if(print_W) call print_excitation_energies('phRPA@qsGW','singlet',nS,Om) ! Compute correlation part of the self-energy @@ -227,7 +225,8 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS tex=tex+tex2-tex1 call wall_time(tsrg1) - call SRG_self_energy(flow,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call SRG_self_energy(nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call wall_time(tsrg2) @@ -245,14 +244,16 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS ! Compute commutator and convergence criteria - error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) + err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) + + if(nSCF > 1) Conv = maxval(abs(err)) ! DIIS extrapolation if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F) end if @@ -271,7 +272,7 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS ! Save quasiparticles energy for next cycle - Conv = maxval(abs(error)) + Conv = maxval(abs(err)) eOld(:) = eGW(:) !------------------------------------------------------------------------ @@ -301,8 +302,8 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS ! Print results call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, & - SigC, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGW, dipole) + call print_qsRGW(nBas,nOrb,nO,nSCF,Conv,thresh,eHF,eGW,c, & + SigC,Z,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,dipole) end do !------------------------------------------------------------------------ @@ -319,26 +320,26 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) - deallocate(c, cp, P, F, Fp, J, K, SigC, Z, Om, XpY, XmY, rho, error, error_diis, F_diis) + deallocate(c,cp,P,F,Fp,J,K,SigC,Z,Om,XpY,XmY,rho,err,err_diis,F_diis) stop end if - print *, "Wall time for Fock and exchange build", tt - print *, "Wall Time for AO to MO", tao - print *, "Wall Time for LR", tlr - print *, "Wall Time for excitation density", tex - print *, "Wall Time for SRG", tsrg - print *, "Wall time MO to AO Sigma", tmo +! print *, "Wall time for Fock and exchange build", tt +! print *, "Wall Time for AO to MO", tao +! print *, "Wall Time for LR", tlr +! print *, "Wall Time for excitation density", tex +! print *, "Wall Time for SRG", tsrg +! print *, "Wall time MO to AO Sigma", tmo ! Cumulant expansion - call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z) +! call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z) ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,Z,Om,XpY,XmY,rho,error,error_diis,F_diis) + deallocate(c,cp,P,F,Fp,J,K,SigC,Z,Om,XpY,XmY,rho,err,err_diis,F_diis) ! Perform BSE calculation diff --git a/src/GW/SRG_qsUGW.f90 b/src/GW/SRG_qsUGW.f90 index 6939afb..825eebd 100644 --- a/src/GW/SRG_qsUGW.f90 +++ b/src/GW/SRG_qsUGW.f90 @@ -223,7 +223,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS end do end if - call USRG_self_energy(eta,nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,EcGM,SigC,Z) + call USRG_self_energy(nBas,nC,nO,nV,nR,nSt,eGW,Om,rho,EcGM,SigC,Z) ! Make correlation self-energy Hermitian and transform it back to AO basis diff --git a/src/GW/SRG_self_energy.f90 b/src/GW/SRG_self_energy.f90 index 04e27b4..2b585c3 100644 --- a/src/GW/SRG_self_energy.f90 +++ b/src/GW/SRG_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) +subroutine SRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Compute correlation part of the self-energy @@ -7,7 +7,6 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Input variables - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC integer,intent(in) :: nO @@ -25,6 +24,7 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) integer :: m double precision :: Dpim,Dqim,Dpam,Dqam,Diam double precision :: t1,t2 + double precision :: s ! Output variables @@ -32,6 +32,10 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) double precision,intent(out) :: SigC(nBas,nBas) double precision,intent(out) :: Z(nBas) +! SRG flow parameter + + s = 500d0 + ! Initialize SigC(:,:) = 0d0 @@ -45,7 +49,7 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) call wall_time(t1) !$OMP PARALLEL & - !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) & + !$OMP SHARED(SigC,rho,s,nS,nC,nO,nBas,nR,e,Om) & !$OMP PRIVATE(m,i,q,p,Dpim,Dqim) & !$OMP DEFAULT(NONE) !$OMP DO @@ -55,7 +59,7 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) do i=nC+1,nO Dpim = e(p) - e(i) + Om(m) Dqim = e(q) - e(i) + Om(m) - SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,m)*rho(q,i,m)*(1d0-dexp(-eta*Dpim*Dpim)*dexp(-eta*Dqim*Dqim)) & + SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,m)*rho(q,i,m)*(1d0-dexp(-s*Dpim*Dpim)*dexp(-s*Dqim*Dqim)) & *(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim) end do end do @@ -64,14 +68,14 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) !$OMP END DO !$OMP END PARALLEL - call wall_time(t2) - print *, "first loop", (t2-t1) +! call wall_time(t2) +! print *, "first loop", (t2-t1) ! Virtual part of the correlation self-energy call wall_time(t1) !$OMP PARALLEL & - !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nR,nBas,e,Om) & + !$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nBas,e,Om) & !$OMP PRIVATE(m,a,q,p,Dpam,Dqam) & !$OMP DEFAULT(NONE) !$OMP DO @@ -81,7 +85,7 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) do a=nO+1,nBas-nR Dpam = e(p) - e(a) - Om(m) Dqam = e(q) - e(a) - Om(m) - SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,m)*rho(q,a,m)*(1d0-exp(-eta*Dpam*Dpam)*exp(-eta*Dqam*Dqam)) & + SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,m)*rho(q,a,m)*(1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam)) & *(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam) end do end do @@ -90,8 +94,8 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) !$OMP END DO !$OMP END PARALLEL - call wall_time(t2) - print *, "second loop", (t2-t1) +! call wall_time(t2) +! print *, "second loop", (t2-t1) ! Initialize @@ -102,7 +106,7 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) do i=nC+1,nO do m=1,nS Dpim = e(p) - e(i) + Om(m) - Z(p) = Z(p) - 2d0*rho(p,i,m)**2*(1d0-dexp(-2d0*eta*Dpim*Dpim))/Dpim**2 + Z(p) = Z(p) - 2d0*rho(p,i,m)**2*(1d0-dexp(-2d0*s*Dpim*Dpim))/Dpim**2 end do end do end do @@ -113,7 +117,7 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) do a=nO+1,nBas-nR do m=1,nS Dpam = e(p) - e(a) - Om(m) - Z(p) = Z(p) - 2d0*rho(p,a,m)**2*(1d0-dexp(-2d0*eta*Dpam*Dpam))/Dpam**2 + Z(p) = Z(p) - 2d0*rho(p,a,m)**2*(1d0-dexp(-2d0*s*Dpam*Dpam))/Dpam**2 end do end do end do @@ -129,7 +133,7 @@ subroutine SRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) do a=nO+1,nBas-nR do m=1,nS Diam = e(a) - e(i) + Om(m) - EcGM = EcGM - 4d0*rho(a,i,m)*rho(a,i,m)*(1d0-exp(-2d0*eta*Diam*Diam))/Diam + EcGM = EcGM - 4d0*rho(a,i,m)*rho(a,i,m)*(1d0-exp(-2d0*s*Diam*Diam))/Diam end do end do end do diff --git a/src/GW/USRG_self_energy.f90 b/src/GW/USRG_self_energy.f90 index 4dbd0fe..18820b0 100644 --- a/src/GW/USRG_self_energy.f90 +++ b/src/GW/USRG_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) +subroutine USRG_self_energy(nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Compute correlation part of the self-energy @@ -7,7 +7,6 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Input variables - double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC(nspin) integer,intent(in) :: nO(nspin) @@ -26,6 +25,7 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) integer :: m double precision :: Dpim,Dqim,Dpam,Dqam,Diam double precision :: t1,t2 + double precision :: s ! Output variables @@ -33,6 +33,10 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) double precision,intent(out) :: SigC(nBas,nBas,nspin) double precision,intent(out) :: Z(nBas,nspin) +! SRG flow parameter + + s = 500d0 + ! Initialize SigC(:,:,:) = 0d0 @@ -43,10 +47,10 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) ! Occupied part of the correlation self-energy - call wall_time(t1) +! call wall_time(t1) !$OMP PARALLEL & - !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nBas,nR,e,Om) & + !$OMP SHARED(SigC,rho,s,nS,nC,nO,nBas,nR,e,Om) & !$OMP PRIVATE(ispin,m,i,q,p,Dpim,Dqim) & !$OMP DEFAULT(NONE) !$OMP DO @@ -58,7 +62,7 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) Dpim = e(p,ispin) - e(i,ispin) + Om(m) Dqim = e(q,ispin) - e(i,ispin) + Om(m) SigC(p,q,ispin) = SigC(p,q,ispin) & - + rho(p,i,m,ispin)*rho(q,i,m,ispin)*(1d0-dexp(-eta*Dpim*Dpim)*dexp(-eta*Dqim*Dqim)) & + + rho(p,i,m,ispin)*rho(q,i,m,ispin)*(1d0-dexp(-s*Dpim*Dpim)*dexp(-s*Dqim*Dqim)) & *(Dpim + Dqim)/(Dpim*Dpim + Dqim*Dqim) end do end do @@ -68,14 +72,14 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) !$OMP END DO !$OMP END PARALLEL - call wall_time(t2) - print *, "first loop", (t2-t1) +! call wall_time(t2) +! print *, "first loop", (t2-t1) ! Virtual part of the correlation self-energy call wall_time(t1) !$OMP PARALLEL & - !$OMP SHARED(SigC,rho,eta,nS,nC,nO,nR,nBas,e,Om) & + !$OMP SHARED(SigC,rho,s,nS,nC,nO,nR,nBas,e,Om) & !$OMP PRIVATE(ispin,m,a,q,p,Dpam,Dqam) & !$OMP DEFAULT(NONE) !$OMP DO @@ -87,7 +91,7 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) Dpam = e(p,ispin) - e(a,ispin) - Om(m) Dqam = e(q,ispin) - e(a,ispin) - Om(m) SigC(p,q,ispin) = SigC(p,q,ispin) & - + rho(p,a,m,ispin)*rho(q,a,m,ispin)*(1d0-exp(-eta*Dpam*Dpam)*exp(-eta*Dqam*Dqam)) & + + rho(p,a,m,ispin)*rho(q,a,m,ispin)*(1d0-exp(-s*Dpam*Dpam)*exp(-s*Dqam*Dqam)) & *(Dpam + Dqam)/(Dpam*Dpam + Dqam*Dqam) end do end do @@ -97,8 +101,8 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) !$OMP END DO !$OMP END PARALLEL - call wall_time(t2) - print *, "second loop", (t2-t1) +! call wall_time(t2) +! print *, "second loop", (t2-t1) ! Initialize @@ -110,7 +114,7 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) do i=nC(ispin)+1,nO(ispin) do m=1,nS Dpim = e(p,ispin) - e(i,ispin) + Om(m) - Z(p,ispin) = Z(p,ispin) - rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*eta*Dpim*Dpim))/Dpim**2 + Z(p,ispin) = Z(p,ispin) - rho(p,i,m,ispin)**2*(1d0-dexp(-2d0*s*Dpim*Dpim))/Dpim**2 end do end do end do @@ -123,7 +127,7 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) do a=nO(ispin)+1,nBas-nR(ispin) do m=1,nS Dpam = e(p,ispin) - e(a,ispin) - Om(m) - Z(p,ispin) = Z(p,ispin) - rho(p,a,m,ispin)**2*(1d0-dexp(-2d0*eta*Dpam*Dpam))/Dpam**2 + Z(p,ispin) = Z(p,ispin) - rho(p,a,m,ispin)**2*(1d0-dexp(-2d0*s*Dpam*Dpam))/Dpam**2 end do end do end do @@ -141,7 +145,7 @@ subroutine USRG_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,Z) do a=nO(ispin)+1,nBas-nR(ispin) do m=1,nS Diam = e(a,ispin) - e(i,ispin) + Om(m) - EcGM = EcGM - rho(a,i,m,ispin)*rho(a,i,m,ispin)*(1d0-exp(-2d0*eta*Diam*Diam))/Diam + EcGM = EcGM - rho(a,i,m,ispin)*rho(a,i,m,ispin)*(1d0-exp(-2d0*s*Diam*Diam))/Diam end do end do end do diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index 5596225..61fc434 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -195,7 +195,6 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om) - ! Compute correlation part of the self-energy call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho) @@ -266,7 +265,6 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop c = matmul(c,cp) endif - call AOtoMO(nBas,nOrb,c,SigCp,SigC) ! Density matrix