PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00
Code Releases Activity

Sph

Browse Source
This commit is contained in:
Pierre-Francois Loos 2019-04-02 22:58:24 +02:00
parent 8ad895ad85
commit 31e3707ce6
11 changed files with 157 additions and 119 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

18
GoSph
View File

@ -1,17 +1,17 @@
#! /bin/bash
if [ $# -ne 1 ]
if [ $# -ne 2 ]
then
echo "You need one argument [BasisSetSize] !!"
echo "You need two arguments [Number of electrons] [BasisSetSize] !!"
fi
if [ $# = 1 ]
if [ $# = 2 ]
then
cp examples/molecule.Sph input/molecule
cp examples/basis.Sph.Ylm"$1" input/basis
cp int/Sph_ERI_"$1".dat int/ERI.dat
cp int/Sph_Kin_"$1".dat int/Kin.dat
cp int/Sph_Nuc_"$1".dat int/Nuc.dat
cp int/Sph_Ov_"$1".dat int/Ov.dat
cp examples/molecule.Sph_"$1" input/molecule
cp examples/basis.Sph.Ylm"$2" input/basis
cp int/Sph_ERI_"$2".dat int/ERI.dat
cp int/Sph_Kin_"$2".dat int/Kin.dat
cp int/Sph_Nuc_"$2".dat int/Nuc.dat
cp int/Sph_Ov_"$2".dat int/Ov.dat
./bin/QuAcK
fi

4
examples/molecule.sph
View File

@ -1,4 +0,0 @@
# nAt nEla nElb nCore nRyd
1 16 16 0 0
# Znuc x y z
X 0.0 0.0 0.0

26
input/basis
View File

@ -1,19 +1,9 @@
1 9
1 3
S 3 1.00
38.3600000 0.0238090
5.7700000 0.1548910
1.2400000 0.4699870
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
S 1 1.00
1.0000000 1.0000000
0.2976000 1.0000000
P 1 1.00
1.2750000 1.0000000

4
input/molecule
View File

@ -1,4 +1,4 @@
# nAt nEla nElb nCore nRyd
1 16 16 0 0
1 1 1 0 0
# Znuc x y z
X 0.0 0.0 0.0
He 0.0 0.0 0.0

2
input/options
View File

@ -9,6 +9,6 @@
# GF: maxSCF thresh DIIS n_diis renormalization
64 0.00001 T 5 3
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize
64 0.00001 T 5 F F F F F F F
64 0.00001 T 5 F F T F F F F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

34
input/weight
View File

@ -1,29 +1,9 @@
1 6
S 8 1.00
2940.0000000 0.0006800
441.2000000 0.0052360
100.5000000 0.0266060
28.4300000 0.0999930
9.1690000 0.2697020
3.1960000 0.4514690
1.1590000 0.2950740
0.1811000 0.0125870
S 8 1.00
2940.0000000 -0.0001230
441.2000000 -0.0009660
100.5000000 -0.0048310
28.4300000 -0.0193140
9.1690000 -0.0532800
3.1960000 -0.1207230
1.1590000 -0.1334350
0.1811000 0.5307670
1 3
S 3 1.00
38.3600000 0.0238090
5.7700000 0.1548910
1.2400000 0.4699870
S 1 1.00
0.0589000 1.0000000
P 3 1.00
3.6190000 0.0291110
0.7110000 0.1693650
0.1951000 0.5134580
0.2976000 1.0000000
P 1 1.00
0.0601800 1.0000000
D 1 1.00
0.2380000 1.0000000
1.2750000 1.0000000

42
src/QuAcK/G0W0.f90
View File

@ -1,5 +1,5 @@
subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,P,cHF,eHF,eG0W0)
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
! Perform G0W0 calculation
@ -8,19 +8,34 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
! Input variables
logical,intent(in) :: COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold
logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
logical,intent(in) :: BSE
logical,intent(in) :: TDA
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ENuc,ERHF
double precision,intent(in) :: cHF(nBas,nBas),eHF(nBas),Hc(nBas,nBas),P(nBas,nBas)
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),ERI_MO_basis(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO_basis(nBas,nBas,nBas,nBas)
! Local variables
logical :: dRPA
integer :: ispin
double precision :: EcRPA,EcGM
double precision,allocatable :: H(:,:),SigmaC(:),Z(:)
double precision,allocatable :: Omega(:,:),XpY(:,:,:),rho(:,:,:,:),rhox(:,:,:,:)
double precision :: EcRPA
double precision :: EcGM
double precision,allocatable :: H(:,:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
! Output variables
@ -49,13 +64,12 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
! Memory allocation
allocate(H(nBas,nBas),SigmaC(nBas),Z(nBas), &
Omega(nS,nspin),XpY(nS,nS,nspin), &
allocate(H(nBas,nBas),SigC(nBas),Z(nBas),Omega(nS,nspin),XpY(nS,nS,nspin), &
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin))
! Compute Hartree Hamiltonian in the MO basis
call Hartree_matrix_MO_basis(nBas,cHF,P,Hc,ERI_AO_basis,H)
call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO_basis,H)
! Compute linear response
@ -69,7 +83,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigmaC)
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
! COHSEX static approximation
@ -85,12 +99,12 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
! Solve the quasi-particle equation
eG0W0(:) = eHF(:) + Z(:)*SigmaC(:)
eG0W0(:) = eHF(:) + Z(:)*SigC(:)
! Dump results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigmaC,Z,eG0W0,EcRPA,EcGM)
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eG0W0,EcRPA,EcGM)
! Plot stuff
@ -118,7 +132,7 @@ subroutine G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
ispin = 2
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eHF,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,cHF,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eG0W0,ERI_MO_basis, &
rho(:,:,:,1),EcRPA,Omega(:,ispin),XpY(:,:,ispin))

4
src/QuAcK/QuAcK.f90
View File

@ -405,7 +405,7 @@ program QuAcK
call cpu_time(start_G0W0)
call G0W0(COHSEX,SOSEX,BSE,TDA,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0
@ -440,7 +440,7 @@ program QuAcK
call cpu_time(start_qsGW)
call qsGW(maxSCF_GW,thresh_GW,n_diis_GW, &
COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold, &
nBas,nC,nO,nV,nR,nS,ENuc,S,X,T,V,Hc,ERI_AO_basis,PHF,cHF,eHF)
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,S,X,T,V,Hc,ERI_AO_basis,PHF,cHF,eHF)
call cpu_time(end_qsGW)
t_qsGW = end_qsGW - start_qsGW

76
src/QuAcK/evGW.f90
View File

@ -8,24 +8,52 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Input variables
integer,intent(in) :: maxSCF,max_diis
double precision,intent(in) :: thresh,ENuc,ERHF
logical,intent(in) :: COHSEX,SOSEX,BSE,TDA,G0W,GW0
logical,intent(in) :: singlet_manifold,triplet_manifold,linearize
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
double precision,intent(in) :: ENuc
double precision,intent(in) :: ERHF
logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
logical,intent(in) :: BSE
logical,intent(in) :: TDA
logical,intent(in) :: G0W
logical,intent(in) :: GW0
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
logical,intent(in) :: linearize
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: cHF(nBas,nBas),eHF(nBas),eG0W0(nBas),Hc(nBas,nBas),PHF(nBas,nBas)
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),ERI_MO_basis(nBas,nBas,nBas,nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: eG0W0(nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_MO_basis(nBas,nBas,nBas,nBas)
! Local variables
logical :: dRPA,linear_mixing
integer :: ispin,nSCF,n_diis
logical :: dRPA
logical :: linear_mixing
integer :: ispin
integer :: nSCF
integer :: n_diis
double precision :: rcond
double precision :: Conv,EcRPA,EcGM,lambda
double precision,allocatable :: error_diis(:,:),e_diis(:,:)
double precision,allocatable :: eGW(:),eOld(:),Z(:)
double precision,allocatable :: H(:,:),SigmaC(:)
double precision,allocatable :: Omega(:,:),XpY(:,:,:),rho(:,:,:,:),rhox(:,:,:,:)
double precision :: Conv
double precision :: EcRPA
double precision :: EcGM
double precision :: lambda
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: e_diis(:,:)
double precision,allocatable :: eGW(:)
double precision,allocatable :: eOld(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: H(:,:)
double precision,allocatable :: SigC(:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
! Hello world
@ -51,10 +79,8 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Memory allocation
allocate(eGW(nBas),eOld(nBas),Z(nBas), &
H(nBas,nBas),SigmaC(nBas), &
Omega(nS,nspin),XpY(nS,nS,nspin), &
rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
allocate(eGW(nBas),eOld(nBas),Z(nBas),H(nBas,nBas),SigC(nBas),Omega(nS,nspin), &
XpY(nS,nS,nspin),rho(nBas,nBas,nS,nspin),rhox(nBas,nBas,nS,nspin), &
error_diis(nBas,max_diis),e_diis(nBas,max_diis))
! Initialization
@ -90,22 +116,22 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Compute correlation part of the self-energy
call excitation_density(nBas,nC,nO,nR,nS,cHF,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,cHF,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
! Correlation self-energy
if(G0W) then
call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eHF, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigmaC)
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eHF,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
else
call self_energy_correlation_diag(COHSEX,SOSEX,nBas,nC,nO,nV,nR,nS,eGW, &
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigmaC)
Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),EcGM,SigC)
call renormalization_factor(SOSEX,nBas,nC,nO,nV,nR,nS,eGW,Omega(:,ispin),rho(:,:,:,ispin),rhox(:,:,:,ispin),Z)
endif
@ -114,11 +140,11 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
if(linearize) then
eGW(:) = eHF(:) + Z(:)*SigmaC(:)
eGW(:) = eHF(:) + Z(:)*SigC(:)
else
eGW(:) = eHF(:) + SigmaC(:)
eGW(:) = eHF(:) + SigC(:)
endif
@ -129,7 +155,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Print results
call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigmaC,Z,eGW,EcRPA)
call print_evGW(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA)
! Linear mixing or DIIS extrapolation
@ -199,7 +225,7 @@ subroutine evGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
ispin = 2
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,cHF,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,eGW,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))

3
src/QuAcK/excitation_density.f90
View File

@ -1,4 +1,4 @@
subroutine excitation_density(nBas,nC,nO,nR,nS,c,ERI,XpY,rho)
subroutine excitation_density(nBas,nC,nO,nR,nS,ERI,XpY,rho)
! Compute excitation densities
@ -7,7 +7,6 @@ subroutine excitation_density(nBas,nC,nO,nR,nS,c,ERI,XpY,rho)
! Input variables
integer,intent(in) :: nBas,nC,nO,nR,nS
double precision,intent(in) :: c(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: XpY(nS,nS)

63
src/QuAcK/qsGW.f90
View File

@ -8,28 +8,61 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Input variables
integer,intent(in) :: maxSCF,max_diis
integer,intent(in) :: maxSCF
integer,intent(in) :: max_diis
double precision,intent(in) :: thresh
logical,intent(in) :: COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_manifold,triplet_manifold
logical,intent(in) :: COHSEX
logical,intent(in) :: SOSEX
logical,intent(in) :: BSE
logical,intent(in) :: TDA
logical,intent(in) :: G0W
logical,intent(in) :: GW0
logical,intent(in) :: singlet_manifold
logical,intent(in) :: triplet_manifold
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
double precision,intent(in) :: ENuc
double precision,intent(in) :: PHF(nBas,nBas),cHF(nBas,nBas),eHF(nBas)
double precision,intent(in) :: S(nBas,nBas),T(nBas,nBAs),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas)
double precision,intent(in) :: eHF(nBas)
double precision,intent(in) :: cHF(nBas,nBas)
double precision,intent(in) :: PHF(nBas,nBas)
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBAs)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas)
! Local variables
logical :: dRPA
integer :: nSCF,nBasSq,ispin,n_diis
double precision :: EcRPA,EcGM,Conv
integer :: nSCF
integer :: nBasSq
integer :: ispin
integer :: n_diis
double precision :: EcRPA
double precision :: EcGM
double precision :: Conv
double precision :: rcond
double precision,external :: trace_matrix
double precision,allocatable :: error_diis(:,:),F_diis(:,:)
double precision,allocatable :: Omega(:,:),XpY(:,:,:),rho(:,:,:,:),rhox(:,:,:,:)
double precision,allocatable :: c(:,:),cp(:,:),e(:),P(:,:)
double precision,allocatable :: F(:,:),Fp(:,:),J(:,:),K(:,:)
double precision,allocatable :: SigC(:,:),SigCp(:,:),SigCm(:,:),Z(:)
double precision,allocatable :: error(:,:),ERI_MO_basis(:,:,:,:)
double precision,allocatable :: error_diis(:,:)
double precision,allocatable :: F_diis(:,:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: rho(:,:,:,:)
double precision,allocatable :: rhox(:,:,:,:)
double precision,allocatable :: c(:,:)
double precision,allocatable :: cp(:,:)
double precision,allocatable :: e(:)
double precision,allocatable :: P(:,:)
double precision,allocatable :: F(:,:)
double precision,allocatable :: Fp(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: K(:,:)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: SigCp(:,:)
double precision,allocatable :: SigCm(:,:)
double precision,allocatable :: Z(:)
double precision,allocatable :: ERI_MO_basis(:,:,:,:)
double precision,allocatable :: error(:,:)
! Hello world
@ -101,8 +134,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
! Compute correlation part of the self-energy
call excitation_density(nBas,nC,nO,nR,nS,c,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,c,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
if(SOSEX) call excitation_density_SOSEX(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rhox(:,:,:,ispin))
if(G0W) then
@ -214,7 +247,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,COHSEX,SOSEX,BSE,TDA,G0W,GW0,singlet_mani
ispin = 2
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,c,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call excitation_density(nBas,nC,nO,nR,nS,ERI_MO_basis,XpY(:,:,ispin),rho(:,:,:,ispin))
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI_MO_basis, &
rho(:,:,:,ispin),EcRPA,Omega(:,ispin),XpY(:,:,ispin))