diff --git a/README.md b/README.md
index 5ed8fa7..61be55b 100644
--- a/README.md
+++ b/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
+
+
+
+
+
+
+---
-**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
+
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).
+
+---
diff --git a/input/methods.default b/input/methods.default
index bc212bd..5ed2560 100644
--- a/input/methods.default
+++ b/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
diff --git a/src/GF/RG0F2.f90 b/src/GF/RG0F2.f90
index afbbb8f..63b5453 100644
--- a/src/GF/RG0F2.f90
+++ b/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
diff --git a/src/GF/RGF.f90 b/src/GF/RGF.f90
index 8665c02..07acb38 100644
--- a/src/GF/RGF.f90
+++ b/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, &
- singlet,triplet,linearize,eta,regularize,nBas,nOrb,nC,nO,nV,nR,nS, &
+ call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, &
+ 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, &
- dBSE,dTDA,singlet,triplet,eta,regularize,nNuc,ZNuc, &
- rNuc,ENuc,nBas,nOrb,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc, &
+ call qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA, &
+ 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)
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)
diff --git a/src/GF/RGF2_QP_graph.f90 b/src/GF/RGF2_QP_graph.f90
index d4d49bd..6250668 100644
--- a/src/GF/RGF2_QP_graph.f90
+++ b/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
- logical,intent(in) :: reg
+ double precision,intent(in) :: eta
+ 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
- SigC = RGF2_SRG_SigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eOld,ERI)
- dSigC = RGF2_SRG_dSigC(flow,p,w,eta,nBas,nC,nO,nV,nR,eOld,ERI)
+ if(doSRG) then
+ SigC = RGF2_SRG_SigC(p,w,flow,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)
diff --git a/src/GF/RGF2_SRG_SigC.f90 b/src/GF/RGF2_SRG_SigC.f90
index c3809d9..0880523 100644
--- a/src/GF/RGF2_SRG_SigC.f90
+++ b/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
diff --git a/src/GF/RGF2_SRG_dSigC.f90 b/src/GF/RGF2_SRG_dSigC.f90
index 7c4bf25..ae0b9fc 100644
--- a/src/GF/RGF2_SRG_dSigC.f90
+++ b/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)&
- *(eps**2 - eta**2)/(eps**2 + eta**2)**2
+ 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
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)&
- *(eps**2 - eta**2)/(eps**2 + eta**2)**2
+ 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
end do
end do
diff --git a/src/GF/RGF2_SRG_self_energy.f90 b/src/GF/RGF2_SRG_self_energy.f90
index 2623934..9aed51f 100644
--- a/src/GF/RGF2_SRG_self_energy.f90
+++ b/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)&
- /(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
+ SigC(p,q) = SigC(p,q) + num*(eps_p + eps_q)/(eps_p**2 + eps_q**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)&
- /(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
+ SigC(p,q) = SigC(p,q) + num*(eps_p + eps_q)/(eps_p**2 + eps_q**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
diff --git a/src/GF/RGF2_SRG_self_energy_diag.f90 b/src/GF/RGF2_SRG_self_energy_diag.f90
index cfe67f9..d58f5a2 100644
--- a/src/GF/RGF2_SRG_self_energy_diag.f90
+++ b/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)
- Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+ SigC(p) = SigC(p) + num/eps
+ 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)
- Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2
+ SigC(p) = SigC(p) + num/eps
+ 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
diff --git a/src/GF/RGF2_self_energy.f90 b/src/GF/RGF2_self_energy.f90
index ec09c0b..b1be4da 100644
--- a/src/GF/RGF2_self_energy.f90
+++ b/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
diff --git a/src/GF/RGF2_self_energy_diag.f90 b/src/GF/RGF2_self_energy_diag.f90
index 71822ab..977dcf3 100644
--- a/src/GF/RGF2_self_energy_diag.f90
+++ b/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
diff --git a/src/GF/evRGF2.f90 b/src/GF/evRGF2.f90
index 2489db8..60d8d7e 100644
--- a/src/GF/evRGF2.f90
+++ b/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
diff --git a/src/GF/print_qsRGF2.f90 b/src/GF/print_qsRGF2.f90
index feb590f..908e9a1 100644
--- a/src/GF/print_qsRGF2.f90
+++ b/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(*,*)
diff --git a/src/GF/qsRGF2.f90 b/src/GF/qsRGF2.f90
index 63999d2..499a3e9 100644
--- a/src/GF/qsRGF2.f90
+++ b/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
- n_diis = 0
- ispin = 1
- Conv = 1d0
- P(:,:) = PHF(:,:)
- eOld(:) = eHF(:)
- eGF(:) = eHF(:)
- c(:,:) = cHF(:,:)
- F_diis(:,:) = 0d0
- error_diis(:,:) = 0d0
- rcond = 0d0
- flow = 500d0
+ nSCF = -1
+ n_diis = 0
+ ispin = 1
+ Conv = 1d0
+ P(:,:) = PHF(:,:)
+ eGF(:) = eHF(:)
+ c(:,:) = cHF(:,:)
+ F_diis(:,:) = 0d0
+ err_diis(:,:) = 0d0
+ rcond = 0d0
!------------------------------------------------------------------------
! 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 MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1))
+ call AOtoMO(nBas,nOrb,c,F,Fp)
+ 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
-
- 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
- !------------------------------------------------------------------------
+ Conv = maxval(abs(err))
! 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))
-
- ! Correlation energy
-
- call RMP2(.false., regularize, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, EqsGF2, eGF, Ec)
+ Ex = 0.25d0*trace_matrix(nBas,matmul(P,K))
! 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 print_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, &
- c, SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF2, dipole)
+ call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole)
+ call print_qsRGF2(nBas,nOrb,nO,nSCF,Conv,thresh,eHF,eGF,&
+ 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, &
- nV, nR, nS, ERI_MO, dipole_int_MO, eGF, EcBSE)
+ call RGF2_phBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,&
+ 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, &
- nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE)
+ call RGF2_ppBSE(TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,&
+ nC,nO,nV,nR,ERI_MO,dipole_int_MO,eGF,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
diff --git a/src/GW/RGW_SRG_self_energy.f90 b/src/GW/RGW_SRG_self_energy.f90
index 53b4521..52209a4 100644
--- a/src/GW/RGW_SRG_self_energy.f90
+++ b/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 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
- end do
- end do
+ do a=nO+1,nOrb-nR
+ 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
+
+ end do
+ end do
end do
Z(:) = 1d0/(1d0 - Z(:))
diff --git a/src/GW/evRGW.f90 b/src/GW/evRGW.f90
index 90677ce..3967dfd 100644
--- a/src/GW/evRGW.f90
+++ b/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
diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90
index 7adf22b..9f81545 100644
--- a/src/GW/qsRGW.f90
+++ b/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 MOtoAO(nBas,nOrb,S(1,1),c(1,1),Fp(1,1),F(1,1))
+ call AOtoMO(nBas,nOrb,c,F,Fp)
+ 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)