starting cleaning SRG

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
 

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.
+  logical                       :: print_W = .false.
-  double precision,allocatable  :: error_diis(:,:)
+  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,allocatable  :: err(:,:)
-
-  double precision,parameter    :: flow = 500d0
 
 ! 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,7 +149,7 @@ 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
@@ -165,7 +163,7 @@ subroutine SRG_qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS
   eOld(:)       = eHF(:)
   c(:,:)        = cHF(:,:)
   F_diis(:,:)   = 0d0
-  error_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(:)
 
     !------------------------------------------------------------------------
@@ -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 Fock and exchange build", tt
- print *, "Wall Time for AO to MO", tao
+! print *, "Wall Time for AO to MO", tao
- print *, "Wall Time for LR", tlr
+! print *, "Wall Time for LR", tlr
- print *, "Wall Time for excitation density", tex
+! print *, "Wall Time for excitation density", tex
- print *, "Wall Time for SRG", tsrg
+! print *, "Wall Time for SRG", tsrg
- print *, "Wall time MO to AO Sigma", tmo
+! 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
 

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
    

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)
+! call wall_time(t2)
- print *, "first loop", (t2-t1)
+! 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)
+! call wall_time(t2)
-  print *, "second loop", (t2-t1)
+! 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

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)
+! call wall_time(t2)
- print *, "first loop", (t2-t1)
+! 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)
+! call wall_time(t2)
-  print *, "second loop", (t2-t1)
+! 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

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