LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 12:23:50 +01:00
Code Releases Activity

eta in linear response

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-01-23 21:22:41 +01:00
parent 76aa00962d
commit 9e3c3f9c89
29 changed files with 221 additions and 56756 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.BF
View File

@ -2,4 +2,4 @@
2 7 7 0 0
# Znuc x y z
B 0. 0. 0.
F 0. 0. 3.4
F 0. 0. 2.41

4
examples/molecule.CO
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd
2 7 7 4 0
2 7 7 0 0
# Znuc x y z
C 0. 0. 0.
O 0. 0. 3.4
O 0. 0. 2.154

2
examples/molecule.F2
View File

@ -2,4 +2,4 @@
2 9 9 0 0
# Znuc x y z
F 0. 0. 0.
F 0. 0. 3.4
F 0. 0. 3

2
examples/molecule.H2
View File

@ -2,4 +2,4 @@
2 1 1 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 3
H 0. 0. 1.4

4
examples/molecule.HCl
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd
2 9 9 10 0
2 9 9 0 0
# Znuc x y z
H 0. 0. 0.
Cl 0. 0. 3.4
Cl 0. 0. 2.42

2
examples/molecule.LiF
View File

@ -2,4 +2,4 @@
2 6 6 0 0
# Znuc x y z
Li 0. 0. 0.
F 0. 0. 3.4
F 0. 0. 3.003

2
examples/molecule.LiH
View File

@ -2,4 +2,4 @@
2 2 2 0 0
# Znuc x y z
Li 0. 0. 0.
H 0. 0. 3.9
H 0. 0. 3.05

4
examples/molecule.N2
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd
2 7 7 4 0
2 7 7 0 0
# Znuc x y z
N 0. 0. 0.
N 0. 0. 3.4
N 0. 0. 2.067

109
input/basis
View File

@ -1,75 +1,58 @@
1 10
1 6
S 8
1 11420.0000000 0.0005230
2 1712.0000000 0.0040450
3 389.3000000 0.0207750
4 110.0000000 0.0807270
5 35.5700000 0.2330740
6 12.5400000 0.4335010
7 4.6440000 0.3474720
8 0.5118000 -0.0085080
1 14710.0000000 0.0007210
2 2207.0000000 0.0055530
3 502.8000000 0.0282670
4 142.6000000 0.1064440
5 46.4700000 0.2868140
6 16.7000000 0.4486410
7 6.3560000 0.2647610
8 1.3160000 0.0153330
S 8
1 11420.0000000 -0.0001150
2 1712.0000000 -0.0008950
3 389.3000000 -0.0046240
4 110.0000000 -0.0185280
5 35.5700000 -0.0573390
6 12.5400000 -0.1320760
7 4.6440000 -0.1725100
8 0.5118000 0.5999440
1 14710.0000000 -0.0001650
2 2207.0000000 -0.0013080
3 502.8000000 -0.0064950
4 142.6000000 -0.0266910
5 46.4700000 -0.0736900
6 16.7000000 -0.1707760
7 6.3560000 -0.1123270
8 1.3160000 0.5628140
S 1
1 1.2930000 1.0000000
S 1
1 0.1787000 1.0000000
1 0.3897000 1.0000000
P 3
1 26.6300000 0.0146700
2 5.9480000 0.0917640
3 1.7420000 0.2986830
1 22.6700000 0.0448780
2 4.9770000 0.2357180
3 1.3470000 0.5085210
P 1
1 0.5550000 1.0000000
P 1
1 0.1725000 1.0000000
1 0.3471000 1.0000000
D 1
1 1.6540000 1.0000000
D 1
1 0.4690000 1.0000000
F 1
1 1.0930000 1.0000000
2 10
1 1.6400000 1.0000000
2 6
S 8
1 11420.0000000 0.0005230
2 1712.0000000 0.0040450
3 389.3000000 0.0207750
4 110.0000000 0.0807270
5 35.5700000 0.2330740
6 12.5400000 0.4335010
7 4.6440000 0.3474720
8 0.5118000 -0.0085080
1 14710.0000000 0.0007210
2 2207.0000000 0.0055530
3 502.8000000 0.0282670
4 142.6000000 0.1064440
5 46.4700000 0.2868140
6 16.7000000 0.4486410
7 6.3560000 0.2647610
8 1.3160000 0.0153330
S 8
1 11420.0000000 -0.0001150
2 1712.0000000 -0.0008950
3 389.3000000 -0.0046240
4 110.0000000 -0.0185280
5 35.5700000 -0.0573390
6 12.5400000 -0.1320760
7 4.6440000 -0.1725100
8 0.5118000 0.5999440
1 14710.0000000 -0.0001650
2 2207.0000000 -0.0013080
3 502.8000000 -0.0064950
4 142.6000000 -0.0266910
5 46.4700000 -0.0736900
6 16.7000000 -0.1707760
7 6.3560000 -0.1123270
8 1.3160000 0.5628140
S 1
1 1.2930000 1.0000000
S 1
1 0.1787000 1.0000000
1 0.3897000 1.0000000
P 3
1 26.6300000 0.0146700
2 5.9480000 0.0917640
3 1.7420000 0.2986830
1 22.6700000 0.0448780
2 4.9770000 0.2357180
3 1.3470000 0.5085210
P 1
1 0.5550000 1.0000000
P 1
1 0.1725000 1.0000000
1 0.3471000 1.0000000
D 1
1 1.6540000 1.0000000
D 1
1 0.4690000 1.0000000
F 1
1 1.0930000 1.0000000
1 1.6400000 1.0000000

6
input/molecule
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd
2 7 7 4 0
2 9 9 0 0
# Znuc x y z
N 0. 0. 0.
N 0. 0. 3.4
F 0. 0. 0.
F 0. 0. 3

4
input/molecule.xyz
View File

@ -1,4 +1,4 @@
2
N 0.0000000000 0.0000000000 0.0000000000
N 0.0000000000 0.0000000000 1.7992026466
F 0.0000000000 0.0000000000 0.0000000000
F 0.0000000000 0.0000000000 1.5875317470

2
input/options
View File

@ -1,5 +1,5 @@
# RHF: maxSCF thresh DIIS n_diis guess_type ortho_type
64 0.0000001 T 5 2 1
64 0.0000001 T 5 1 1
# MP:
# CC: maxSCF thresh DIIS n_diis

109
input/weight
View File

@ -1,75 +1,58 @@
1 10
1 6
S 8
1 11420.0000000 0.0005230
2 1712.0000000 0.0040450
3 389.3000000 0.0207750
4 110.0000000 0.0807270
5 35.5700000 0.2330740
6 12.5400000 0.4335010
7 4.6440000 0.3474720
8 0.5118000 -0.0085080
1 14710.0000000 0.0007210
2 2207.0000000 0.0055530
3 502.8000000 0.0282670
4 142.6000000 0.1064440
5 46.4700000 0.2868140
6 16.7000000 0.4486410
7 6.3560000 0.2647610
8 1.3160000 0.0153330
S 8
1 11420.0000000 -0.0001150
2 1712.0000000 -0.0008950
3 389.3000000 -0.0046240
4 110.0000000 -0.0185280
5 35.5700000 -0.0573390
6 12.5400000 -0.1320760
7 4.6440000 -0.1725100
8 0.5118000 0.5999440
1 14710.0000000 -0.0001650
2 2207.0000000 -0.0013080
3 502.8000000 -0.0064950
4 142.6000000 -0.0266910
5 46.4700000 -0.0736900
6 16.7000000 -0.1707760
7 6.3560000 -0.1123270
8 1.3160000 0.5628140
S 1
1 1.2930000 1.0000000
S 1
1 0.1787000 1.0000000
1 0.3897000 1.0000000
P 3
1 26.6300000 0.0146700
2 5.9480000 0.0917640
3 1.7420000 0.2986830
1 22.6700000 0.0448780
2 4.9770000 0.2357180
3 1.3470000 0.5085210
P 1
1 0.5550000 1.0000000
P 1
1 0.1725000 1.0000000
1 0.3471000 1.0000000
D 1
1 1.6540000 1.0000000
D 1
1 0.4690000 1.0000000
F 1
1 1.0930000 1.0000000
2 10
1 1.6400000 1.0000000
2 6
S 8
1 11420.0000000 0.0005230
2 1712.0000000 0.0040450
3 389.3000000 0.0207750
4 110.0000000 0.0807270
5 35.5700000 0.2330740
6 12.5400000 0.4335010
7 4.6440000 0.3474720
8 0.5118000 -0.0085080
1 14710.0000000 0.0007210
2 2207.0000000 0.0055530
3 502.8000000 0.0282670
4 142.6000000 0.1064440
5 46.4700000 0.2868140
6 16.7000000 0.4486410
7 6.3560000 0.2647610
8 1.3160000 0.0153330
S 8
1 11420.0000000 -0.0001150
2 1712.0000000 -0.0008950
3 389.3000000 -0.0046240
4 110.0000000 -0.0185280
5 35.5700000 -0.0573390
6 12.5400000 -0.1320760
7 4.6440000 -0.1725100
8 0.5118000 0.5999440
1 14710.0000000 -0.0001650
2 2207.0000000 -0.0013080
3 502.8000000 -0.0064950
4 142.6000000 -0.0266910
5 46.4700000 -0.0736900
6 16.7000000 -0.1707760
7 6.3560000 -0.1123270
8 1.3160000 0.5628140
S 1
1 1.2930000 1.0000000
S 1
1 0.1787000 1.0000000
1 0.3897000 1.0000000
P 3
1 26.6300000 0.0146700
2 5.9480000 0.0917640
3 1.7420000 0.2986830
1 22.6700000 0.0448780
2 4.9770000 0.2357180
3 1.3470000 0.5085210
P 1
1 0.5550000 1.0000000
P 1
1 0.1725000 1.0000000
1 0.3471000 1.0000000
D 1
1 1.6540000 1.0000000
D 1
1 0.4690000 1.0000000
F 1
1 1.0930000 1.0000000
1 1.6400000 1.0000000

1
int/.gitignore vendored
View File

@ -1 +0,0 @@
*.dat

56513
int/boxints.txt
View File

File diff suppressed because it is too large Load Diff

8
scan_H2.sh
View File

@ -1,10 +1,10 @@
#! /bin/bash
MOL="H2"
BASIS="VTZ"
R_START=0.5
R_END=3.0
DR=0.05
BASIS="cc-pvqz"
R_START=1.397
R_END=1.400
DR=0.001
for R in $(seq $R_START $DR $R_END)
do

6
scan_LiF.sh
View File

@ -2,9 +2,9 @@
MOL="LiF"
BASIS="cc-pvdz"
R_START=1.5
R_END=3.5
DR=0.1
R_START=3.001
R_END=3.003
DR=0.001
for R in $(seq $R_START $DR $R_END)
do

6
scan_LiH.sh
View File

@ -2,9 +2,9 @@
MOL="LiH"
BASIS="cc-pvtz"
R_START=1.5
R_END=4.0
DR=0.1
R_START=3.00
R_END=3.05
DR=0.001
for R in $(seq $R_START $DR $R_END)
do

8
scan_N2.sh
View File

@ -1,10 +1,10 @@
#! /bin/bash
MOL="N2"
BASIS="cc-pvqz"
R_START=1.5
R_END=3.5
DR=0.1
BASIS="cc-pvtz"
R_START=2.060
R_END=2.080
DR=0.001
for R in $(seq $R_START $DR $R_END)
do

25
src/QuAcK/ACFDT.f90
View File

@ -1,4 +1,4 @@
subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_manifold, &
subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
@ -17,6 +17,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
@ -74,13 +75,13 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man
if(doXBS) then
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
end if
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
@ -119,20 +120,20 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_man
lambda = rAC(iAC)
if(doXBS) then
! if(doXBS) then
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
! call linear_response(ispin,dRPA,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
! rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
! call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
end if
! end if
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
! call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
! rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
! call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
! write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
end do

11
src/QuAcK/Bethe_Salpeter.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
! Compute the Bethe-Salpeter excitation energies
@ -12,6 +12,7 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
double precision,intent(in) :: eta
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eW(nBas)
double precision,intent(in) :: eGW(nBas)
@ -38,11 +39,11 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
ispin = 1
EcBSE(ispin) = 0d0
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,.true.,TDA,.true.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))
@ -55,11 +56,11 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
ispin = 2
EcBSE(ispin) = 0d0
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,.true.,TDA,.true.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('BSE ',ispin,nS,Omega(:,ispin))

7
src/QuAcK/Bethe_Salpeter_A_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A_lr)
subroutine Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response
@ -8,6 +8,7 @@ subroutine Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A_lr
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nS)
@ -16,6 +17,7 @@ subroutine Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A_lr
! Local variables
double precision :: chi
double precision :: eps
integer :: i,j,a,b,ia,jb,kc
! Output variables
@ -33,7 +35,8 @@ subroutine Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A_lr
chi = 0d0
do kc=1,nS
chi = chi + rho(i,j,kc)*rho(a,b,kc)/Omega(kc)
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,j,kc)*rho(a,b,kc)*Omega(kc)/eps
enddo
A_lr(ia,jb) = A_lr(ia,jb) - lambda*ERI(i,a,j,b) + 4d0*lambda*chi

7
src/QuAcK/Bethe_Salpeter_B_matrix.f90
View File

@ -1,4 +1,4 @@
subroutine Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B_lr)
subroutine Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B_lr)
! Compute the extra term for Bethe-Salpeter equation for linear response
@ -8,6 +8,7 @@ subroutine Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B_lr
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: eta
double precision,intent(in) :: lambda
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: Omega(nS)
@ -16,6 +17,7 @@ subroutine Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B_lr
! Local variables
double precision :: chi
double precision :: eps
integer :: i,j,a,b,ia,jb,kc
! Output variables
@ -33,7 +35,8 @@ subroutine Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B_lr
chi = 0d0
do kc=1,nS
chi = chi + rho(i,b,kc)*rho(a,j,kc)/Omega(kc)
eps = Omega(kc)**2 + eta**2
chi = chi + rho(i,b,kc)*rho(a,j,kc)*Omega(kc)/eps
enddo
B_lr(ia,jb) = B_lr(ia,jb) - lambda*ERI(i,a,b,j) + 4d0*lambda*chi

32
src/QuAcK/G0W0.f90
View File

@ -78,7 +78,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
! Compute linear response
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
! Compute correlation part of the self-energy
@ -106,15 +106,15 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
! Compute the RPA correlation energy
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@G0W0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@G0W0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -126,7 +126,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
if(BSE) then
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eHF,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
if(exchange_kernel) then
@ -138,10 +138,10 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -161,7 +161,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho,EcAC)
if(exchange_kernel) then
@ -173,10 +173,10 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,singlet_manif
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

6
src/QuAcK/QuAcK.f90
View File

@ -442,7 +442,8 @@ program QuAcK
if(doRPA) then
call cpu_time(start_RPA)
call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
call cpu_time(end_RPA)
t_RPA = end_RPA - start_RPA
@ -458,7 +459,8 @@ program QuAcK
if(doRPAx) then
call cpu_time(start_RPAx)
call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
call cpu_time(end_RPAx)
t_RPAx = end_RPAx - start_RPAx

26
src/QuAcK/RPA.f90
View File

@ -1,4 +1,5 @@
subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e)
subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,e)
! Perform a direct random phase approximation calculation
@ -12,6 +13,7 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC
logical,intent(in) :: exchange_kernel
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC
integer,intent(in) :: nO
@ -57,7 +59,7 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC
ispin = 1
call linear_response(ispin,.true.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
@ -69,7 +71,7 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC
ispin = 2
call linear_response(ispin,.true.,.false.,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
call linear_response(ispin,.true.,.false.,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,e,ERI,rho, &
EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
@ -84,10 +86,10 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'Tr@RPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F15.6)') 'Tr@RPA correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPA correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -100,7 +102,7 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC
write(*,*) '------------------------------------------------------'
write(*,*)
call ACFDT(exchange_kernel,.false.,.true.,.false.,.false.,singlet_manifold,triplet_manifold, &
call ACFDT(exchange_kernel,.false.,.true.,.false.,.false.,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho,EcAC)
if(exchange_kernel) then
@ -112,10 +114,10 @@ subroutine RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,nBas,nC
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'AC@RPA correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F15.6)') 'AC@RPA correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@RPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@RPA correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@RPA total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

30
src/QuAcK/evGW.f90
View File

@ -115,7 +115,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif
@ -217,10 +217,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@evGW total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@evGW total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -228,15 +228,15 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
if(BSE) then
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -256,15 +256,15 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI,eGW,Omega,XpY,XmY,rho,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@evGW total energy =',ENuc + ERHF + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)

7
src/QuAcK/linear_response.f90
View File

@ -1,4 +1,4 @@
subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho,EcRPA,Omega,XpY,XmY)
subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,rho,EcRPA,Omega,XpY,XmY)
! Compute linear response
@ -8,6 +8,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
! Input variables
logical,intent(in) :: dRPA,TDA,BSE
double precision,intent(in) :: eta
integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: lambda
double precision,intent(in) :: e(nBas)
@ -40,7 +41,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
! Build A and B matrices
call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A)
if(BSE) call Bethe_Salpeter_A_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A)
if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A)
! Tamm-Dancoff approximation
@ -48,7 +49,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,r
if(.not. TDA) then
call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B)
if(BSE) call Bethe_Salpeter_B_matrix(nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B)
if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B)
endif

30
src/QuAcK/qsGW.f90
View File

@ -141,7 +141,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
if(.not. GW0 .or. nSCF == 0) then
call linear_response(ispin,.true.,TDA,.false.,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
endif
@ -248,10 +248,10 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F15.6)') 'Tr@RPA@qsGW total energy =',ENuc + EqsGW + EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy (singlet) =',EcRPA(1)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy (triplet) =',EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW correlation energy =',EcRPA(1) + EcRPA(2)
write(*,'(2X,A50,F20.10)') 'Tr@RPA@qsGW total energy =',ENuc + EqsGW + EcRPA(1) + EcRPA(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -259,15 +259,15 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
if(BSE) then
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold, &
call Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F15.6)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (singlet) =',EcBSE(1)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy (triplet) =',EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW correlation energy =',EcBSE(1) + EcBSE(2)
write(*,'(2X,A50,F20.10)') 'Tr@BSE@qsGW total energy =',ENuc + EqsGW + EcBSE(1) + EcBSE(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
@ -287,15 +287,15 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,
end if
call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold, &
call ACFDT(exchange_kernel,doXBS,.true.,TDA,BSE,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eGW,Omega,XpY,XmY,rho,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'
write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F15.6)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (singlet) =',EcAC(1)
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy (triplet) =',EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW correlation energy =',EcAC(1) + EcAC(2)
write(*,'(2X,A50,F20.10)') 'AC@BSE@qsGW total energy =',ENuc + EqsGW + EcAC(1) + EcAC(2)
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)