Merge branch 'master' of github.com:pfloos/QuAcK

README.md

@@ -1,8 +1,19 @@
-# QuAcK: an open-source software for emerging quantum electronic structure methods
+
+# 🦆 QuAcK: an open-source software for emerging quantum electronic structure methods
+
+![License: GPL v3](https://img.shields.io/badge/License-GPLv3-blue.svg)
+![Fortran 90](https://img.shields.io/badge/language-Fortran%2090-yellow)
+![Stars](https://img.shields.io/github/stars/pfloos/QuAcK?style=social)
+![Forks](https://img.shields.io/github/forks/pfloos/QuAcK?style=social)
+
+---
 
 <img src="logo/logo_quack.png"  width="250">
 
-**Contributors:**
+---
+
+## 👥 Contributors
+
 - [Pierre-Francois Loos](https://pfloos.github.io/WEB_LOOS)
 - [Anthony Scemama](https://scemama.github.io)
 - [Enzo Monino](https://enzomonino.github.io)
@@ -10,8 +21,19 @@
 - [Abdallah Ammar](https://scholar.google.com/citations?user=y437T5sAAAAJ&hl=en)
 - [Mauricio Rodriguez-Mayorga](https://scholar.google.com/citations?user=OLGOgQgAAAAJ&hl=es)
 - [Loris Burth](https://github.com/lburth)
-  
+
-# What is it?
+---
+
+## 🚀 Features
+
+- **Rapid Prototyping:** Ideal for testing and developing new quantum chemistry methods.
+- **Modular Design:** Easily integrate with other tools and libraries.
+- **Educational Tool:** Serves as an excellent entry point for researchers familiar with electronic structure theory.
+- **Integration with PySCF:** Utilizes [PySCF](https://github.com/pyscf/pyscf) for computing one- and two-electron integrals.
+
+---
+
+## What is it?
 
 **QuAcK** is a lightweight electronic structure program written in `Fortran 90`, developed at the [Laboratoire de Chimie et Physique Quantiques (LCPQ)](https://www.lcpq.ups-tlse.fr) in Toulouse, France. Originally designed as a platform for rapid prototyping of new ideas in quantum chemistry, QuAcK serves as a flexible and accessible environment for testing novel methods before integrating them more efficiently into larger-scale projects such as the [Quantum Package](https://quantumpackage.github.io/qp2/).
 
@@ -19,11 +41,12 @@ Thanks to its compact and transparent codebase, QuAcK is particularly well-suite
 
 For beginners in the field or those with less programming experience, we recommend starting with [qcmath](https://github.com/LCPQ/qcmath/), a symbolic and numerical quantum chemistry toolkit built in [Mathematica](https://www.wolfram.com/mathematica/). qcmath is specifically designed to help newcomers explore and develop ideas without the complexity of full-fledged numerical implementations.
 
-QuAcK is under active and ongoing development, which means that bugs, inconsistencies, and incomplete features are to be expected. It is a tool made *by* experts *for* experts—users are expected to understand what they’re doing and to remain cautious when interpreting results. The code may allow questionable inputs or behavior *on purpose*, to encourage flexibility during prototyping—so always double-check your results and assumptions.
+> ⚠️ **Note:** QuAcK is under active and ongoing development, which means that bugs, inconsistencies, and incomplete features are to be expected. It is a tool made *by* experts *for* experts—users are expected to understand what they’re doing and to remain cautious when interpreting results. The code may allow questionable inputs or behavior *on purpose*, to encourage flexibility during prototyping—so always double-check your results and assumptions. In short: use QuAcK at your own risk—but also to your advantage, if you're ready to experiment and explore.
 
-In short: use QuAcK at your own risk—but also to your advantage, if you're ready to experiment and explore.
+---
+
+## 🛠️ Installation
 
-# Installation guide
 The QuAcK software can be downloaded on GitHub as a Git repository
 ```
 git clone https://github.com/pfloos/QuAcK.git
@@ -41,7 +64,12 @@ pip install pyscf
 PySCF is used for the computation of one- and two-electron integrals (mainly) which are dumped in files and read by QuAcK.
 Therefore, it is very easy to use other software to compute the integrals or to add other types of integrals.
 
-# Quick start
+Then, go to the `src` directory and compile
+```
+cd src; make
+```
+
+## ⚡ Quick Start
 
 ```
 ~ 💩 % cd $QUACK_ROOT
@@ -145,6 +173,24 @@ QuAcK 💩 % python PyDuck.py -x water -b cc-pvdz
 
 QuAcK runs calculations in the `QUACK_ROOT` directory which is quite unusual but it also use the `--working_dir` option to run calculations elsewhere.
 
+---
+
+## 📄 License
+
+QuAcK is licensed under the [GNU General Public License v3.0](https://www.gnu.org/licenses/gpl-3.0.en.html).
+
+---
+
+## 📫 Contact
+
+For questions or contributions, please open an issue or submit a pull request on the [GitHub repository](https://github.com/pfloos/QuAcK).
+
+--- 
+
+## 💰 Funding
+
 <img src="https://lcpq.github.io/PTEROSOR/img/ERC.png" width="200" />
 
 QuAcK is supported by the [PTEROSOR](https://lcpq.github.io/PTEROSOR/) project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 863481).
+
+---

input/methods.default

@@ -13,7 +13,7 @@
 # G0F2 evGF2 qsGF2 ufGF2 G0F3 evGF3
   F    F     F     F    F    F
 # G0W0 evGW qsGW ufG0W0 ufGW 
-  T    F    F    F      F    
+  F    F    F    F      F    
 # G0T0pp evGTpp qsGTpp ufG0T0pp
   F      F      F      F
 # G0T0eh evGTeh qsGTeh

src/GF/RG0F2.f90

@@ -1,4 +1,4 @@
-subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,regularize, &
+subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,eta,doSRG, &
                  nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform a one-shot second-order Green function calculation
@@ -19,7 +19,7 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
   logical,intent(in)            :: triplet
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
+  logical,intent(in)            :: doSRG
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nOrb
@@ -52,7 +52,16 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
   write(*,*)'*******************************'
   write(*,*)
 
-  flow =  500d0
+! SRG regularization
+
+  flow = 500d0
+
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized G0F2 scheme ***'
+    write(*,*)
+
+  end if
 
 ! Memory allocation
 
@@ -60,15 +69,17 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
 
 ! Frequency-dependent second-order contribution
 
-  if(regularize) then 
+  if(doSRG) then 
 
-    call RGF2_SRG_self_energy_diag(flow,eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z)
+    call RGF2_SRG_self_energy_diag(flow,nOrb,nC,nO,nV,nR,eHF,ERI,Ec,SigC,Z)
 
   else
 
-    call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,SigC,Z)
+    call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eHF,ERI,Ec,SigC,Z)
 
   end if
+
+  print*,Ec
   
   eGFlin(:) = eHF(:) + Z(:)*SigC(:)
 
@@ -83,13 +94,12 @@ subroutine RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize,
     write(*,*) ' *** Quasiparticle energies obtained by root search *** '
     write(*,*)
 
-    call RGF2_QP_graph(flow,regularize,eta,nOrb,nC,nO,nV,nR,eHF,ERI,eGFlin,eHF,eGF,Z)
+    call RGF2_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,eHF,ERI,eGFlin,eHF,eGF,Z)
 
   end if
 
   ! Print results
 
-  call RMP2(.false.,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
   call print_RG0F2(nOrb,nO,eHF,SigC,eGF,Z,ENuc,ERHF,Ec)
 
 ! Perform BSE@GF2 calculation

src/GF/RGF.f90

@@ -1,7 +1,7 @@
 subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,          &
                renorm,maxSCF,                                                           &
                thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet,linearize, &
-               eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,        &
+               eta,doSRG,nNuc,ZNuc,rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,             &
                S,X,T,V,Hc,ERI_AO,ERI_MO,CAP,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
 ! Green's function module
@@ -34,7 +34,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,
   logical,intent(in)            :: triplet
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
+  logical,intent(in)            :: doSRG
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)
@@ -76,7 +76,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,
 
     call wall_time(start_GF)
     call RG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, &
-               linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,    &
+               linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,         &
                ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
     call wall_time(end_GF)
 
@@ -93,8 +93,8 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,
   if(doevGF2) then
 
     call wall_time(start_GF)
-    call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,       & 
+    call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,  & 
-                singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS, &
+                singlet,triplet,linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS, &
                 ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
     call wall_time(end_GF)
 
@@ -111,9 +111,9 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,
   if(doqsGF2) then 
 
     call wall_time(start_GF)
-    call qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,    &
+    call qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
-                dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc,   &
+                dBSE,dTDA,singlet,triplet,eta,doSRG,nNuc,ZNuc,      &
-                rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,   & 
+                rNuc,ENuc,nBas,nOrb,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 wall_time(end_GF)
 
@@ -179,7 +179,7 @@ subroutine RGF(dotest,doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,docG0F2,
 
     call wall_time(start_GF)
     call cRG0F2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,singlet,triplet, &
-               linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,    &
+               linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,          &
                ENuc,ERHF,ERI_MO,CAP,dipole_int_MO,eHF)
     call wall_time(end_GF)
 

src/GF/RGF2_QP_graph.f90

@@ -1,4 +1,4 @@
-subroutine RGF2_QP_graph(flow,reg,eta,nBas,nC,nO,nV,nR,eHF,ERI,eGFlin,eOld,eGF,Z)
+subroutine RGF2_QP_graph(doSRG,eta,flow,nBas,nC,nO,nV,nR,eHF,ERI,eGFlin,eOld,eGF,Z)
 
 ! Compute the graphical solution of the GF2 QP equation
 
@@ -7,8 +7,9 @@ subroutine RGF2_QP_graph(flow,reg,eta,nBas,nC,nO,nV,nR,eHF,ERI,eGFlin,eOld,eGF,Z
 
 ! Input variables
 
-  double precision,intent(in)   :: eta,flow
+  double precision,intent(in)   :: eta
-  logical,intent(in)            :: reg
+  double precision,intent(in)   :: flow
+  logical,intent(in)            :: doSRG
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
@@ -51,9 +52,9 @@ subroutine RGF2_QP_graph(flow,reg,eta,nBas,nC,nO,nV,nR,eHF,ERI,eGFlin,eOld,eGF,Z
     
       nIt = nIt + 1
 
-      if(reg) then
+      if(doSRG) then
-        SigC  = RGF2_SRG_SigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eOld,ERI)
+        SigC  = RGF2_SRG_SigC(p,w,flow,nBas,nC,nO,nV,nR,eOld,ERI)
-        dSigC = RGF2_SRG_dSigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eOld,ERI)
+        dSigC = RGF2_SRG_dSigC(p,w,flow,nBas,nC,nO,nV,nR,eOld,ERI)
       else
         SigC  = RGF2_SigC(p,w,eta,nBas,nC,nO,nV,nR,eOld,ERI)
         dSigC = RGF2_dSigC(p,w,eta,nBas,nC,nO,nV,nR,eOld,ERI)

src/GF/RGF2_SRG_SigC.f90

@@ -1,4 +1,4 @@
-double precision function RGF2_SRG_SigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI)
+double precision function RGF2_SRG_SigC(p,w,flow,nBas,nC,nO,nV,nR,eHF,ERI)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -8,9 +8,8 @@ double precision function RGF2_SRG_SigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI)
 ! Input variables
 
   integer,intent(in)            :: p
-  double precision,intent(in)   :: flow
   double precision,intent(in)   :: w
-  double precision,intent(in)   :: eta
+  double precision,intent(in)   :: flow
   integer,intent(in)            :: nBas,nC,nO,nV,nR
   double precision,intent(in)   :: eHF(nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
@@ -31,7 +30,7 @@ double precision function RGF2_SRG_SigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI)
 
         eps = w + eHF(a) - eHF(i) - eHF(j)
         kappa = 1d0 - exp(-2d0*eps**2*s)
-        RGF2_SRG_SigC = RGF2_SRG_SigC + kappa*(2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)*eps/(eps**2 + eta**2)
+        RGF2_SRG_SigC = RGF2_SRG_SigC + kappa*(2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)/eps
 
         end do
      end do
@@ -45,7 +44,7 @@ double precision function RGF2_SRG_SigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI)
 
         eps = w + eHF(i) - eHF(a) - eHF(b)
         kappa = 1d0 - exp(-2d0*eps**2*s)
-        RGF2_SRG_SigC = RGF2_SRG_SigC + kappa*(2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)*eps/(eps**2 + eta**2)
+        RGF2_SRG_SigC = RGF2_SRG_SigC + kappa*(2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)/eps
 
       end do
     end do

src/GF/RGF2_SRG_dSigC.f90

@@ -1,4 +1,4 @@
-double precision function RGF2_SRG_dSigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI)
+double precision function RGF2_SRG_dSigC(p,w,flow,nBas,nC,nO,nV,nR,eHF,ERI)
 
 ! Compute diagonal of the correlation part of the self-energy
 
@@ -10,7 +10,6 @@ double precision function RGF2_SRG_dSigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI)
   integer,intent(in)            :: p
   double precision,intent(in)   :: w
   double precision,intent(in)   :: flow
-  double precision,intent(in)   :: eta
   integer,intent(in)            :: nBas,nC,nO,nV,nR
   double precision,intent(in)   :: eHF(nBas)
   double precision,intent(in)   :: ERI(nBas,nBas,nBas,nBas)
@@ -33,8 +32,7 @@ double precision function RGF2_SRG_dSigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI)
 
         eps = w + eHF(a) - eHF(i) - eHF(j)
         kappa = 1d0 - exp(-2d0*eps**2*s)
-        RGF2_SRG_dSigC = RGF2_SRG_dSigC - kappa*(2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)&
+        RGF2_SRG_dSigC = RGF2_SRG_dSigC - kappa*(2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)/eps**2
-                *(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
       end do
     end do
@@ -48,8 +46,7 @@ double precision function RGF2_SRG_dSigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eHF,ERI)
 
         eps = w + eHF(i) - eHF(a) - eHF(b)
         kappa = 1d0 - exp(-2d0*eps**2*s)
-        RGF2_SRG_dSigC = RGF2_SRG_dSigC - kappa*(2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)&
+        RGF2_SRG_dSigC = RGF2_SRG_dSigC - kappa*(2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)/eps**2
-                *(eps**2 - eta**2)/(eps**2 + eta**2)**2
 
       end do
     end do

src/GF/RGF2_SRG_self_energy.f90

@@ -1,4 +1,4 @@
-subroutine RGF2_SRG_self_energy(flow, eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
+subroutine RGF2_SRG_self_energy(flow,nBas,nC,nO,nV,nR,e,ERI,Ec,SigC,Z)
 
 ! Compute GF2 self-energy and its renormalization factor
 
@@ -7,7 +7,6 @@ subroutine RGF2_SRG_self_energy(flow, eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
 ! Input variables
 
-  double precision,intent(in)   :: eta
   double precision,intent(in)   :: flow
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
@@ -22,30 +21,28 @@ subroutine RGF2_SRG_self_energy(flow, eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
   integer                       :: i,j,a,b
   integer                       :: p,q
   double precision              :: eps_p,eps_q
-  double precision              :: num
+  double precision              :: num,eps
 
   double precision              :: s
   double precision              :: kappa
 
 ! Output variables
 
+  double precision,intent(out)  :: Ec
   double precision,intent(out)  :: SigC(nBas,nBas)
   double precision,intent(out)  :: Z(nBas)
 
-! Initialize 
-
-  SigC(:,:) = 0d0
-  Z(:)      = 0d0
-
 !-----------------------------------------!
 ! Parameters for regularized calculations !
 !-----------------------------------------!
 
   s = flow
 
-!----------------------------------------------------!
+!-------------------------!
-! Compute GF2 self-energy and renormalization factor !
+! Compute GF2 self-energy !
-!----------------------------------------------------!
+!-------------------------!
+
+  SigC(:,:) = 0d0
 
   do p=nC+1,nBas-nR
     do q=nC+1,nBas-nR
@@ -55,12 +52,10 @@ subroutine RGF2_SRG_self_energy(flow, eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
             eps_p = e(p) + e(a) - e(i) - e(j)
             eps_q = e(q) + e(a) - e(i) - e(j)
-            kappa = 1d0 - exp(-s*(eps_p**2 + eps_q**2 + 2*eta**2))
+            kappa = 1d0 - exp(-s*(eps_p**2 + eps_q**2))
             num = kappa*(2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(q,a,i,j)
 
-            SigC(p,q) = SigC(p,q) + num*(eps_p + eps_q)&
+            SigC(p,q) = SigC(p,q) + num*(eps_p + eps_q)/(eps_p**2 + eps_q**2)
-                         /(eps_p**2 + eps_q**2 + 2*eta**2)
-            if(p == q) Z(p) = Z(p) - num*(eps_p**2 - eta**2)/(eps_p**2 + eta**2)**2
 
           end do
         end do
@@ -76,12 +71,10 @@ subroutine RGF2_SRG_self_energy(flow, eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
             eps_p = e(p) + e(i) - e(a) - e(b)
             eps_q = e(q) + e(i) - e(a) - e(b)
-            kappa = 1d0 - exp(-s*(eps_p**2 + eps_q**2 + 2*eta**2))
+            kappa = 1d0 - exp(-s*(eps_p**2 + eps_q**2))
             num = kappa*(2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(q,i,a,b)
 
-            SigC(p,q) = SigC(p,q) + num*(eps_p + eps_q)&
+            SigC(p,q) = SigC(p,q) + num*(eps_p + eps_q)/(eps_p**2 + eps_q**2)
-                        /(eps_p**2 + eps_q**2 + 2*eta**2)
-            if(p == q) Z(p) = Z(p) - num*(eps_p**2 - eta**2)/(eps_p**2 + eta**2)**2
 
           end do
         end do
@@ -89,6 +82,66 @@ subroutine RGF2_SRG_self_energy(flow, eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
     end do
   end do
 
+!------------------------------------!
+! Compute GF2 renormalization factor !
+!------------------------------------!
+
+  Z(:) = 0d0
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+      do j=nC+1,nO
+        do a=nO+1,nBas-nR
+
+          eps_p = e(p) + e(a) - e(i) - e(j)
+          kappa = 1d0 - exp(-2d0*s*eps_p**2)
+          num = kappa*(2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
+
+          Z(p) = Z(p) - num/eps_p**2
+
+        end do
+      end do
+    end do
+  end do
+
+  do p=nC+1,nBas-nR
+    do i=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
+
+          eps_p = e(p) + e(i) - e(a) - e(b)
+          kappa = 1d0 - exp(-2d0*s*eps_p**2)
+          num = kappa*(2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
+
+          Z(p) = Z(p) - num/eps_p**2
+
+        end do
+      end do
+    end do
+  end do
+
   Z(:) = 1d0/(1d0 - Z(:))
 
+!----------------------------!
+! Compute correlation energy !
+!----------------------------!
+
+  Ec = 0d0
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
+
+          eps = e(i) + e(j) - e(a) - e(b)
+          kappa = 1d0 - exp(-2d0*s*eps**2)
+          num = kappa*(2d0*ERI(i,j,a,b) - ERI(i,j,b,a))*ERI(i,j,a,b)
+
+          Ec = Ec + num/eps
+
+        end do
+      end do
+    end do
+  end do
+
 end subroutine

src/GF/RGF2_SRG_self_energy_diag.f90

@@ -1,4 +1,4 @@
-subroutine RGF2_SRG_self_energy_diag(flow,eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
+subroutine RGF2_SRG_self_energy_diag(flow,nBas,nC,nO,nV,nR,e,ERI,Ec,SigC,Z)
 
 ! Compute diagonal part of the GF2 self-energy and its renormalization factor
 
@@ -7,7 +7,7 @@ subroutine RGF2_SRG_self_energy_diag(flow,eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
 ! Input variables
 
-  double precision,intent(in)   :: eta,flow
+  double precision,intent(in)   :: flow
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nC
   integer,intent(in)            :: nO
@@ -28,6 +28,7 @@ subroutine RGF2_SRG_self_energy_diag(flow,eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
 ! Output variables
 
+  double precision,intent(out)  :: Ec
   double precision,intent(out)  :: SigC(nBas)
   double precision,intent(out)  :: Z(nBas)
 
@@ -55,8 +56,8 @@ subroutine RGF2_SRG_self_energy_diag(flow,eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
           kappa = 1d0 - exp(-2d0*eps**2*s)
           num = kappa*(2d0*ERI(p,a,i,j) - ERI(p,a,j,i))*ERI(p,a,i,j)
 
-          SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2)
+          SigC(p) = SigC(p) + num/eps
-          Z(p)    = Z(p)    - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+          Z(p)    = Z(p)    - num/eps**2
 
         end do
       end do
@@ -72,8 +73,8 @@ subroutine RGF2_SRG_self_energy_diag(flow,eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
           kappa = 1d0 - exp(-2d0*eps**2*s)
           num = kappa*(2d0*ERI(p,i,a,b) - ERI(p,i,b,a))*ERI(p,i,a,b)
 
-          SigC(p) = SigC(p) + num*eps/(eps**2 + eta**2)
+          SigC(p) = SigC(p) + num/eps
-          Z(p)    = Z(p)    - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+          Z(p)    = Z(p)    - num/eps**2
 
         end do
       end do
@@ -81,4 +82,27 @@ subroutine RGF2_SRG_self_energy_diag(flow,eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
   end do
 
   Z(:) = 1d0/(1d0 - Z(:))
+
+!----------------------------!  
+! Compute correlation energy !  
+!----------------------------!  
+
+  Ec = 0d0
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
+
+          eps = e(i) + e(j) - e(a) - e(b)
+          kappa = 1d0 - exp(-2d0*s*eps**2)
+          num = kappa*(2d0*ERI(i,j,a,b) - ERI(i,j,b,a))*ERI(i,j,a,b)
+
+          Ec = Ec + num/eps
+
+        end do
+      end do
+    end do
+  end do
+
 end subroutine 

src/GF/RGF2_self_energy.f90

@@ -1,4 +1,4 @@
-subroutine RGF2_self_energy(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
+subroutine RGF2_self_energy(eta,nBas,nC,nO,nV,nR,e,ERI,Ec,SigC,Z)
 
 ! Compute GF2 self-energy and its renormalization factor
 
@@ -25,6 +25,7 @@ subroutine RGF2_self_energy(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
 ! Output variables
 
+  double precision,intent(out)  :: Ec
   double precision,intent(out)  :: SigC(nBas,nBas)
   double precision,intent(out)  :: Z(nBas)
 
@@ -73,4 +74,23 @@ subroutine RGF2_self_energy(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
   Z(:) = 1d0/(1d0 - Z(:))
 
+! Compute correlation energy
+
+  Ec = 0d0
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
+
+          eps = e(i) + e(j) - e(a) - e(b)
+          num = (2d0*ERI(i,j,a,b) - ERI(i,j,b,a))*ERI(i,j,a,b)
+
+          Ec = Ec + num/eps
+
+        end do
+      end do
+    end do
+  end do
+  
 end subroutine 

src/GF/RGF2_self_energy_diag.f90

@@ -1,4 +1,4 @@
-subroutine RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
+subroutine RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,Ec,SigC,Z)
 
 ! Compute diagonal part of the GF2 self-energy and its renormalization factor
 
@@ -25,6 +25,7 @@ subroutine RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
 
 ! Output variables
 
+  double precision,intent(out)  :: Ec
   double precision,intent(out)  :: SigC(nBas)
   double precision,intent(out)  :: Z(nBas)
 
@@ -66,6 +67,26 @@ subroutine RGF2_self_energy_diag(eta,nBas,nC,nO,nV,nR,e,ERI,SigC,Z)
       end do
     end do
   end do
+
   Z(:) = 1d0/(1d0 - Z(:))
 
+! Compute correlation energy
+
+  Ec = 0d0
+
+  do i=nC+1,nO
+    do j=nC+1,nO
+      do a=nO+1,nBas-nR
+        do b=nO+1,nBas-nR
+
+          eps = e(i) + e(j) - e(a) - e(b)
+          num = (2d0*ERI(i,j,a,b) - ERI(i,j,b,a))*ERI(i,j,a,b)
+
+          Ec = Ec + num/eps
+
+        end do
+      end do
+    end do
+  end do
+
 end subroutine 

src/GF/evRGF2.f90

@@ -1,5 +1,5 @@
 subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,singlet,triplet, &
-                 linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
+                 linearize,eta,doSRG,nBas,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF)
 
 ! Perform eigenvalue self-consistent second-order Green function calculation
 
@@ -22,7 +22,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
   logical,intent(in)            :: triplet
   logical,intent(in)            :: linearize
   double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
+  logical,intent(in)            :: doSRG
 
   integer,intent(in)            :: nBas
   integer,intent(in)            :: nOrb
@@ -55,13 +55,23 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
 
 ! Hello world
 
-
   write(*,*)
   write(*,*)'********************************'
   write(*,*)'* Restricted evGF2 Calculation *'
   write(*,*)'********************************'
   write(*,*)
 
+! SRG regularization
+
+  flow = 500d0
+
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized evGF2 scheme ***'
+    write(*,*)
+
+  end if
+
 ! Memory allocation
 
   allocate(SigC(nOrb), Z(nOrb), eGF(nOrb), eOld(nOrb), error_diis(nOrb,max_diis), e_diis(nOrb,max_diis))
@@ -76,7 +86,6 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
   eGF(:)         = eHF(:)
   eOld(:)         = eHF(:)
   rcond           = 0d0
-  flow            = 500d0
 !------------------------------------------------------------------------
 ! Main SCF loop
 !------------------------------------------------------------------------
@@ -85,13 +94,13 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
 
     ! Frequency-dependent second-order contribution
 
-    if(regularize) then 
+    if(doSRG) then 
 
-      call RGF2_SRG_self_energy_diag(flow,eta,nOrb,nC,nO,nV,nR,eGF,ERI,SigC,Z)
+      call RGF2_SRG_self_energy_diag(flow,nOrb,nC,nO,nV,nR,eGF,ERI,Ec,SigC,Z)
 
     else
 
-      call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eGF,ERI,SigC,Z)
+      call RGF2_self_energy_diag(eta,nOrb,nC,nO,nV,nR,eGF,ERI,Ec,SigC,Z)
 
     end if
 
@@ -104,7 +113,7 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
       write(*,*) ' *** Quasiparticle energies obtained by root search *** '
       write(*,*)
  
-      call RGF2_QP_graph(flow,regularize,eta,nOrb,nC,nO,nV,nR,eHF,ERI,eOld,eOld,eGF,Z)
+      call RGF2_QP_graph(doSRG,eta,flow,nOrb,nC,nO,nV,nR,eHF,ERI,eOld,eOld,eGF,Z)
  
     end if
 
@@ -112,7 +121,6 @@ subroutine evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis,si
 
     ! Print results
 
-    call RMP2(.false.,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eGF,Ec)
     call print_evRGF2(nOrb,nO,nSCF,Conv,eHF,SigC,Z,eGF,ENuc,ERHF,Ec)
 
     ! DIIS extrapolation

src/GF/print_qsRGF2.f90

@@ -97,10 +97,10 @@ subroutine print_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, c, &
     write(*,'(A32,1X,F16.10,A3)') ' Nuclear   repulsion: ',ENuc,' au'
     write(*,'(A32,1X,F16.10,A3)') ' qsGF2        energy: ',ENuc + EqsGF,' au'
     write(*,'(A50)')           '---------------------------------------'
-    write(*,'(A35)')           ' Dipole moment (Debye)    '
+    write(*,'(A32)')           ' 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(*,'(A50)')           '---------------------------------------'
     write(*,*)
  
     write(*,'(A50)') '---------------------------------------'
@@ -109,7 +109,7 @@ subroutine print_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, c, &
     call matout(nBas, nOrb, c)
     write(*,*)
     write(*,'(A50)') '---------------------------------------'
-    write(*,'(A32)') ' qsGF2 MO energies'
+    write(*,'(A32)') ' qsGF2 MO energies    '
     write(*,'(A50)') '---------------------------------------'
     call matout(nOrb, 1, eGF)
     write(*,*)

src/GF/qsRGF2.f90

@@ -1,5 +1,5 @@
 subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
-                  dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc, &
+                  dBSE,dTDA,singlet,triplet,eta,doSRG,nNuc,ZNuc,      &
                   rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, & 
                   ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF)
 
@@ -23,7 +23,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
   logical,intent(in)            :: singlet
   logical,intent(in)            :: triplet
   double precision,intent(in)   :: eta
-  logical,intent(in)            :: regularize
+  logical,intent(in)            :: doSRG
 
   integer,intent(in)            :: nNuc
   double precision,intent(in)   :: ZNuc(nNuc)
@@ -64,12 +64,11 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
   double precision              :: Ec
   double precision              :: EcBSE(nspin)
 
-  double precision,allocatable  :: error_diis(:,:)
+  double precision,allocatable  :: err_diis(:,:)
   double precision,allocatable  :: F_diis(:,:)
   double precision,allocatable  :: c(:,:)
   double precision,allocatable  :: cp(:,:)
   double precision,allocatable  :: eGF(:)
-  double precision,allocatable  :: eOld(:)
   double precision,allocatable  :: P(:,:)
   double precision,allocatable  :: F(:,:)
   double precision,allocatable  :: Fp(:,:)
@@ -78,7 +77,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
   double precision,allocatable  :: SigC(:,:)
   double precision,allocatable  :: SigCp(:,:)
   double precision,allocatable  :: Z(:)
-  double precision,allocatable  :: error(:,:)
+  double precision,allocatable  :: err(:,:)
 
 ! Hello world
 
@@ -89,6 +88,17 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
   write(*,*)'********************************'
   write(*,*)
 
+! SRG regularization
+
+  flow = 500d0
+
+  if(doSRG) then
+
+    write(*,*) '*** SRG regularized qsGF2 scheme ***'
+    write(*,*)
+
+  end if
+
 ! Warning 
 
   write(*,*) '!! ERIs in MO basis will be overwritten in qsGF2 !!'
@@ -108,7 +118,6 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
 ! Memory allocation
 
   allocate(eGF(nOrb))
-  allocate(eOld(nOrb))
 
   allocate(c(nBas,nOrb))
 
@@ -119,36 +128,34 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
   allocate(F(nBas,nBas))
   allocate(J(nBas,nBas))
   allocate(K(nBas,nBas))
-  allocate(error(nBas,nBas))
+  allocate(err(nBas,nBas))
 
   allocate(Z(nOrb))
   allocate(SigC(nOrb,nOrb))
 
   allocate(SigCp(nBas,nBas))
 
-  allocate(error_diis(nBas_Sq,max_diis))
+  allocate(err_diis(nBas_Sq,max_diis))
   allocate(F_diis(nBas_Sq,max_diis))
 
 ! Initialization
   
-  nSCF            = -1
+  nSCF          = -1
-  n_diis          = 0
+  n_diis        = 0
-  ispin           = 1
+  ispin         = 1
-  Conv            = 1d0
+  Conv          = 1d0
-  P(:,:)          = PHF(:,:)
+  P(:,:)        = PHF(:,:)
-  eOld(:)         = eHF(:)
+  eGF(:)        = eHF(:)
-  eGF(:)          = eHF(:)
+  c(:,:)        = cHF(:,:)
-  c(:,:)          = cHF(:,:)
+  F_diis(:,:)   = 0d0
-  F_diis(:,:)     = 0d0
+  err_diis(:,:) = 0d0
-  error_diis(:,:) = 0d0
+  rcond         = 0d0
-  rcond           = 0d0
-  flow = 500d0
 
 !------------------------------------------------------------------------
 ! Main loop
 !------------------------------------------------------------------------
 
-  do while(Conv > thresh .and. nSCF <= maxSCF)
+  do while(Conv > thresh .and. nSCF < maxSCF)
 
     ! Increment
 
@@ -156,115 +163,104 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
 
     ! Buid Hartree matrix
 
-    call Hartree_matrix_AO_basis(nBas, P, ERI_AO, J)
+    call Hartree_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)
+    call exchange_matrix_AO_basis(nBas,P,ERI_AO,K)
 
     ! AO to MO transformation of two-electron integrals
 
-    call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO)
+    call AOtoMO_ERI_RHF(nBas,nOrb,c,ERI_AO,ERI_MO)
 
     ! Compute self-energy and renormalization factor
 
-    if(regularize) then
+    if(doSRG) then
 
-      call RGF2_SRG_self_energy(flow,eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z)
+      call RGF2_SRG_self_energy(flow,nOrb,nC,nO,nV,nR,eGF,ERI_MO,Ec,SigC,Z)
 
     else
 
-      call RGF2_self_energy(eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z)
+      call RGF2_self_energy(eta,nOrb,nC,nO,nV,nR,eGF,ERI_MO,Ec,SigC,Z)
 
     end if
 
+    print*,Ec
+
     ! Make correlation self-energy Hermitian and transform it back to AO basis
    
     SigC = 0.5d0*(SigC + transpose(SigC))
 
-    call MOtoAO(nBas, nOrb, S, c, SigC, SigCp)
+    call MOtoAO(nBas,nOrb,S,c,SigC,SigCp)
  
     ! Solve the quasi-particle equation
 
     F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:)
     if(nBas .ne. nOrb) then
-      call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1))
+      call AOtoMO(nBas,nOrb,c,F,Fp)
-      call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1))
+      call MOtoAO(nBas,nOrb,S,c,Fp,F)
     endif
 
     ! 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)
 
-    ! DIIS extrapolation 
+    Conv = maxval(abs(err))
-
-    n_diis = min(n_diis+1, max_diis)
-    if(abs(rcond) > 1d-7) then
-      call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F)
-    else
-      n_diis = 0
-    end if
-
-    ! Diagonalize Hamiltonian in AO basis
-
-    if(nBas .eq. nOrb) then
-      Fp = matmul(transpose(X), matmul(F, X))
-      cp(:,:) = Fp(:,:)
-      call diagonalize_matrix(nOrb, cp, eGF)
-      c = matmul(X, cp)
-    else
-      Fp = matmul(transpose(c), matmul(F, c))
-      cp(:,:) = Fp(:,:)
-      call diagonalize_matrix(nOrb, cp, eGF)
-      c = matmul(c, cp)
-    endif
-
-
-    ! 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(eGF - eOld))
-    eOld(:) = eGF(:)
-
-    !------------------------------------------------------------------------
-    !   Compute total energy
-    !------------------------------------------------------------------------
 
     ! Kinetic energy
 
-    ET = trace_matrix(nBas, matmul(P, T))
+    ET = trace_matrix(nBas,matmul(P,T))
 
     ! Potential energy
 
-    EV = trace_matrix(nBas, matmul(P, V))
+    EV = trace_matrix(nBas,matmul(P,V))
 
     ! Hartree energy
 
-    EJ = 0.5d0*trace_matrix(nBas, matmul(P, J))
+    EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
 
     ! Exchange energy
 
-    Ex = 0.25d0*trace_matrix(nBas, matmul(P, K))
+    Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
-
-    ! Correlation energy
-
-    call RMP2(.false., regularize, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, EqsGF2, eGF, Ec)
 
     ! Total energy
 
     EqsGF2 = ET + EV + EJ + Ex + Ec
 
+    ! 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,err_diis,F_diis,err,F)
+
+    end if
+
+    ! Diagonalize Hamiltonian in AO basis
+
+    if(nBas == nOrb) then
+      Fp = matmul(transpose(X),matmul(F,X))
+      cp(:,:) = Fp(:,:)
+      call diagonalize_matrix(nOrb,cp,eGF)
+      c = matmul(X,cp)
+    else
+      Fp = matmul(transpose(c),matmul(F,c))
+      cp(:,:) = Fp(:,:)
+      call diagonalize_matrix(nOrb,cp,eGF)
+      c = matmul(c,cp)
+    endif
+
+    call AOtoMO(nBas,nOrb,c,SigCp,SigC)
+
+    ! 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 dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
-    call print_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, &
+    call print_qsRGF2(nBas,nOrb,nO,nSCF,Conv,thresh,eHF,eGF,&
-                      c, SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF2, dipole)
+                      c,SigC,Z,ENuc,ET,EV,EJ,Ex,Ec,EqsGF2,dipole)
 
   end do
 !------------------------------------------------------------------------
@@ -281,21 +277,21 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
     write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
     write(*,*)
 
-    deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
+    deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,err,err_diis,F_diis)
     stop
 
   end if
 
 ! Deallocate memory
 
-  deallocate(c, cp, P, F, Fp, J, K, SigC, SigCp, Z, error, error_diis, F_diis)
+  deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,Z,err,err_diis,F_diis)
 
 ! Perform phBSE@GF2 calculation
 
   if(dophBSE) then
 
-    call RGF2_phBSE(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, nC, nO, &
+    call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,&
-                    nV, nR, nS, ERI_MO, dipole_int_MO, eGF, EcBSE)
+                    nV,nR,nS,ERI_MO,dipole_int_MO,eGF,EcBSE)
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'
@@ -313,8 +309,8 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,  &
 
   if(doppBSE) then
 
-    call RGF2_ppBSE(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, &
+    call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,&
-                    nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE)
+                    nC,nO,nV,nR,ERI_MO,dipole_int_MO,eGF,EcBSE)
 
     write(*,*)
     write(*,*)'-------------------------------------------------------------------------------'

src/GW/RGW_SRG_self_energy.f90

@@ -117,12 +117,14 @@ subroutine RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,e,Om,rho,EcGM,SigC,
   ! Virtual part of the renormlization factor
 
   do p=nC+1,nOrb-nR
-     do a=nO+1,nOrb-nR
+    do a=nO+1,nOrb-nR
-        do m=1,nS
+      do m=1,nS
-           Dpam = e(p) - e(a) - Om(m)
+
-           Z(p) = Z(p)  - 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2
+        Dpam = e(p) - e(a) - Om(m)
-        end do
+        Z(p) = Z(p)  - 2d0*rho(p,a,m)**2*(1d0-exp(-2d0*s*Dpam*Dpam))/Dpam**2
-     end do
+
+      end do
+    end do
   end do
 
   Z(:) = 1d0/(1d0 - Z(:))

src/GW/evRGW.f90

@@ -85,7 +85,7 @@ subroutine evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
 ! SRG regularization
 
-  flow = 100d0
+  flow = 500d0
 
   if(doSRG) then
 

src/GW/qsRGW.f90

@@ -175,7 +175,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
   c(:,:)        = cHF(:,:)
   F_diis(:,:)   = 0d0
   err_diis(:,:) = 0d0
-  rcond          = 0d0
+  rcond         = 0d0
 
 !------------------------------------------------------------------------
 ! Main loop
@@ -211,9 +211,13 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
     call RGW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho)
 
     if(doSRG) then 
+
       call RGW_SRG_self_energy(flow,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+
     else
+
       call RGW_self_energy(eta,nBas,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z)
+
     end if
 
     ! Make correlation self-energy Hermitian and transform it back to AO basis
@@ -226,8 +230,8 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
     F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:)
     if(nBas .ne. nOrb) then
-      call AOtoMO(nBas,nOrb,c(1,1),F(1,1),Fp(1,1))
+      call AOtoMO(nBas,nOrb,c,F,Fp)
-      call MOtoAO(nBas,nOrb,S(1,1),c(1,1),Fp(1,1),F(1,1))
+      call MOtoAO(nBas,nOrb,S,c,Fp,F)
     endif
 
     ! Compute commutator and convergence criteria
@@ -267,7 +271,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop
 
     ! Diagonalize Hamiltonian in AO basis
 
-    if(nBas .eq. nOrb) then
+    if(nBas == nOrb) then
       Fp = matmul(transpose(X),matmul(F,X))
       cp(:,:) = Fp(:,:)
       call diagonalize_matrix(nOrb,cp,eGW)