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

RS in GW

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-04-16 17:02:01 +02:00
parent a0b8383971
commit 2b0244e010
19 changed files with 353 additions and 90 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
GoQCaml
View File

@ -1,5 +1,5 @@
#! /bin/bash
cd int
../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz
###../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz -m 0.5
###../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz
../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz -m 0.5

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.3
F 0. 0. 2

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. 1.4
H 0. 0. 2.5

51
input/basis
View File

@ -1,18 +1,39 @@
1 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 10
S 8
1 24350.0000000 0.0005020
2 3650.0000000 0.0038810
3 829.6000000 0.0199970
4 234.0000000 0.0784180
5 75.6100000 0.2296760
6 26.7300000 0.4327220
7 9.9270000 0.3506420
8 1.1020000 -0.0076450
S 8
1 24350.0000000 -0.0001180
2 3650.0000000 -0.0009150
3 829.6000000 -0.0047370
4 234.0000000 -0.0192330
5 75.6100000 -0.0603690
6 26.7300000 -0.1425080
7 9.9270000 -0.1777100
8 1.1020000 0.6058360
S 1
1 0.1220000 1.0000000
P 1
1 0.7270000 1.0000000
2 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 2.8360000 1.0000000
S 1
1 0.1220000 1.0000000
1 0.3782000 1.0000000
P 3
1 54.7000000 0.0171510
2 12.4300000 0.1076560
3 3.6790000 0.3216810
P 1
1 0.7270000 1.0000000
1 1.1430000 1.0000000
P 1
1 0.3300000 1.0000000
D 1
1 4.0140000 1.0000000
D 1
1 1.0960000 1.0000000
F 1
1 2.5440000 1.0000000

4
input/methods
View File

@ -11,8 +11,8 @@
# G0F2 evGF2 G0F3 evGF3
F F F F
# G0W0 evGW qsGW
F T F
T F F
# G0T0 evGT qsGT
F F F
T F F
# MCMP2
F

5
input/molecule
View File

@ -1,5 +1,4 @@
# nAt nEla nElb nCore nRyd
2 1 1 0 0
1 5 5 0 0
# Znuc x y z
H 0. 0. 0.
H 0. 0. 1.4
Ne 0.0 0.0 0.0

5
input/molecule.xyz
View File

@ -1,4 +1,3 @@
2
1
H 0.0000000000 0.0000000000 0.0000000000
H 0.0000000000 0.0000000000 0.7408481486
Ne 0.0000000000 0.0000000000 0.0000000000

2
input/options
View File

@ -5,7 +5,7 @@
# CC: maxSCF thresh DIIS n_diis
64 0.0000001 T 5
# CIS/TDHF/BSE: singlet triplet
T T
T F
# GF: maxSCF thresh DIIS n_diis lin renorm
256 0.00001 T 5 T 3
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta

51
input/weight
View File

@ -1,18 +1,39 @@
1 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 10
S 8
1 24350.0000000 0.0005020
2 3650.0000000 0.0038810
3 829.6000000 0.0199970
4 234.0000000 0.0784180
5 75.6100000 0.2296760
6 26.7300000 0.4327220
7 9.9270000 0.3506420
8 1.1020000 -0.0076450
S 8
1 24350.0000000 -0.0001180
2 3650.0000000 -0.0009150
3 829.6000000 -0.0047370
4 234.0000000 -0.0192330
5 75.6100000 -0.0603690
6 26.7300000 -0.1425080
7 9.9270000 -0.1777100
8 1.1020000 0.6058360
S 1
1 0.1220000 1.0000000
P 1
1 0.7270000 1.0000000
2 3
S 3
1 13.0100000 0.0196850
2 1.9620000 0.1379770
3 0.4446000 0.4781480
1 2.8360000 1.0000000
S 1
1 0.1220000 1.0000000
1 0.3782000 1.0000000
P 3
1 54.7000000 0.0171510
2 12.4300000 0.1076560
3 3.6790000 0.3216810
P 1
1 0.7270000 1.0000000
1 1.1430000 1.0000000
P 1
1 0.3300000 1.0000000
D 1
1 4.0140000 1.0000000
D 1
1 1.0960000 1.0000000
F 1
1 2.5440000 1.0000000

8
src/QuAcK/G0T0.f90
View File

@ -1,5 +1,5 @@
subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_manifold, &
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF)
linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF,eG0T0)
! Perform one-shot calculation with a T-matrix self-energy (G0T0)
@ -52,8 +52,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
double precision,allocatable :: SigT(:)
double precision,allocatable :: Z(:)
double precision,allocatable :: eG0T0(:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
@ -61,6 +59,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
! Output variables
double precision,intent(out) :: eG0T0(nBas)
! Hello world
write(*,*)
@ -88,7 +88,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nO,nVVt),rho2t(nBas,nV,nOOt), &
SigT(nBas),Z(nBas),eG0T0(nBas))
SigT(nBas),Z(nBas))
!----------------------------------------------
! alpha-beta block

9
src/QuAcK/G0W0.f90
View File

@ -34,6 +34,7 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, &
! Local variables
logical :: print_W = .false.
integer :: ispin
double precision :: EcRPA(nspin)
double precision :: EcBSE(nspin)
@ -104,7 +105,10 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, &
eGWlin(:) = eHF(:) + Z(:)*SigC(:)
if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
eGW(:) = eGWlin(:)
! Find all the roots of the QP equation if necessary
@ -121,7 +125,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, &
! Dump results
! call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
if(print_W) call print_excitation('RPA ',ispin,nS,Omega(:,ispin))
call print_G0W0(nBas,nO,eHF,ENuc,ERHF,SigC,Z,eGW,EcRPA(ispin),EcGM)
! Compute the RPA correlation energy

136
src/QuAcK/QuAcK.f90
View File

@ -14,6 +14,7 @@ program QuAcK
logical :: doG0W0,doevGW,doqsGW
logical :: doG0T0,doevGT,doqsGT
logical :: doMCMP2,doMinMCMP2
logical :: doGTGW = .false.
logical :: doBas
integer :: nNuc,nBas,nBasCABS
@ -46,8 +47,10 @@ program QuAcK
double precision,allocatable :: cTrial(:),gradient(:),hessian(:,:)
double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),H(:,:),X(:,:)
double precision,allocatable :: ERI_AO_basis(:,:,:,:)
double precision,allocatable :: ERI_MO_basis(:,:,:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
double precision,allocatable :: ERI_ERF_AO(:,:,:,:)
double precision,allocatable :: ERI_ERF_MO(:,:,:,:)
double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:)
double precision :: start_QuAcK ,end_QuAcK ,t_QuAcK
@ -202,7 +205,7 @@ program QuAcK
allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas),eG0T0(nBas),PHF(nBas,nBas,nspin), &
S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas), &
ERI_AO_basis(nBas,nBas,nBas,nBas),ERI_MO_basis(nBas,nBas,nBas,nBas))
ERI_AO(nBas,nBas,nBas,nBas),ERI_MO(nBas,nBas,nBas,nBas))
! Read integrals
@ -210,12 +213,12 @@ program QuAcK
if(doSph) then
call read_integrals_sph(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis)
call read_integrals_sph(nEl(:),nBas,S,T,V,Hc,ERI_AO)
else
call system('./GoQCaml')
call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis)
call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO)
end if
@ -237,7 +240,7 @@ program QuAcK
if(doRHF) then
call cpu_time(start_HF)
call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO_basis,X,ENuc,ERHF,eHF,cHF,PHF)
call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,ERHF,eHF,cHF,PHF)
call cpu_time(end_HF)
t_HF = end_HF - start_HF
@ -253,7 +256,7 @@ program QuAcK
if(doUHF) then
call cpu_time(start_HF)
call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO_basis,X,ENuc,EUHF,eHF,cHF,PHF)
call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,EUHF,eHF,cHF,PHF)
call cpu_time(end_HF)
t_HF = end_HF - start_HF
@ -270,7 +273,7 @@ program QuAcK
call cpu_time(start_MOM)
call MOM(maxSCF_HF,thresh_HF,n_diis_HF, &
nBas,nO,S,T,V,Hc,ERI_AO_basis,X,ENuc,ERHF,cHF,eHF,PHF)
nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,ERHF,cHF,eHF,PHF)
call cpu_time(end_MOM)
t_MOM = end_MOM - start_MOM
@ -285,7 +288,7 @@ program QuAcK
! Compute Hartree Hamiltonian in the MO basis
call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO_basis,H)
call Hartree_matrix_MO_basis(nBas,cHF,PHF,Hc,ERI_AO,H)
call cpu_time(start_AOtoMO)
@ -296,13 +299,13 @@ program QuAcK
if(doSph) then
ERI_MO_basis = ERI_AO_basis
ERI_MO(:,:,:,:) = ERI_AO(:,:,:,:)
print*,'!!! MO = AO !!!'
deallocate(ERI_AO_basis)
deallocate(ERI_AO)
else
call AOtoMO_integral_transform(nBas,cHF,ERI_AO_basis,ERI_MO_basis)
call AOtoMO_integral_transform(nBas,cHF,ERI_AO,ERI_MO)
end if
@ -319,7 +322,7 @@ program QuAcK
if(doMP2) then
call cpu_time(start_MP2)
call MP2(nBas,nC,nO,nV,nR,ERI_MO_basis,ENuc,ERHF,eHF,EcMP2)
call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,EcMP2)
call cpu_time(end_MP2)
t_MP2 = end_MP2 - start_MP2
@ -335,7 +338,7 @@ program QuAcK
if(doMP3) then
call cpu_time(start_MP3)
call MP3(nBas,nEl,ERI_MO_basis,eHF,ENuc,ERHF)
call MP3(nBas,nEl,ERI_MO,eHF,ENuc,ERHF)
call cpu_time(end_MP3)
t_MP3 = end_MP3 - start_MP3
@ -358,8 +361,8 @@ program QuAcK
! Read integrals
call read_F12_integrals(nBas,S,ERI_AO_basis,F12,Yuk,FC)
call MP2F12(nBas,nC,nO,nV,ERI_AO_basis,F12,Yuk,FC,ERHF,eHF,cHF)
call read_F12_integrals(nBas,S,ERI_AO,F12,Yuk,FC)
call MP2F12(nBas,nC,nO,nV,ERI_AO,F12,Yuk,FC,ERHF,eHF,cHF)
call cpu_time(end_MP2F12)
t_MP2F12 = end_MP2F12 - start_MP2F12
@ -376,7 +379,7 @@ program QuAcK
call cpu_time(start_CCD)
call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), &
ERI_MO_basis,ENuc,ERHF,eHF)
ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD
@ -395,7 +398,7 @@ program QuAcK
call cpu_time(start_CCSD)
call CCSD(maxSCF_CC,thresh_CC,n_diis_CC,doCCSDT,nBas,nC(1),nO(1),nV(1),nR(1), &
ERI_MO_basis,ENuc,ERHF,eHF)
ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCSD)
t_CCSD = end_CCSD - start_CCSD
@ -412,7 +415,7 @@ program QuAcK
call cpu_time(start_CCD)
call drCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), &
ERI_MO_basis,ENuc,ERHF,eHF)
ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD
@ -429,7 +432,7 @@ program QuAcK
call cpu_time(start_CCD)
call rCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), &
ERI_MO_basis,ENuc,ERHF,eHF)
ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD
@ -446,7 +449,7 @@ program QuAcK
call cpu_time(start_CCD)
call lCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), &
ERI_MO_basis,ENuc,ERHF,eHF)
ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD
@ -463,7 +466,7 @@ program QuAcK
call cpu_time(start_CCD)
call pCCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nC(1),nO(1),nV(1),nR(1), &
ERI_MO_basis,ENuc,ERHF,eHF)
ERI_MO,ENuc,ERHF,eHF)
call cpu_time(end_CCD)
t_CCD = end_CCD - start_CCD
@ -479,7 +482,7 @@ program QuAcK
if(doCIS) then
call cpu_time(start_CIS)
call CIS(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ERI_MO_basis,eHF)
call CIS(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF)
call cpu_time(end_CIS)
t_CIS = end_CIS - start_CIS
@ -496,7 +499,7 @@ program QuAcK
call cpu_time(start_RPA)
call RPA(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_RPA)
t_RPA = end_RPA - start_RPA
@ -513,7 +516,7 @@ program QuAcK
call cpu_time(start_RPAx)
call RPAx(doACFDT,exchange_kernel,singlet_manifold,triplet_manifold,eta, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO_basis,eHF)
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_RPAx)
t_RPAx = end_RPAx - start_RPAx
@ -530,7 +533,7 @@ program QuAcK
call cpu_time(start_ppRPA)
call ppRPA(singlet_manifold,triplet_manifold, &
nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF)
nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO,eHF)
call cpu_time(end_ppRPA)
t_ppRPA = end_ppRPA - start_ppRPA
@ -547,7 +550,7 @@ program QuAcK
call cpu_time(start_ADC)
call ADC(singlet_manifold,triplet_manifold,maxSCF_GF,thresh_GF,n_diis_GF, &
nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO_basis)
nBas,nC(1),nO(1),nV(1),nR(1),eHF,ERI_MO)
call cpu_time(end_ADC)
t_ADC = end_ADC - start_ADC
@ -563,7 +566,7 @@ program QuAcK
if(doG0F2) then
call cpu_time(start_GF2)
call G0F2(linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
call G0F2(linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF)
call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2
@ -579,7 +582,7 @@ program QuAcK
if(doevGF2) then
call cpu_time(start_GF2)
call evGF2(maxSCF_GF,thresh_GF,n_diis_GF,linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
call evGF2(maxSCF_GF,thresh_GF,n_diis_GF,linGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF)
call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2
@ -595,7 +598,7 @@ program QuAcK
if(doG0F3) then
call cpu_time(start_GF3)
call G0F3(renormGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
call G0F3(renormGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF)
call cpu_time(end_GF3)
t_GF3 = end_GF3 - start_GF3
@ -611,7 +614,7 @@ program QuAcK
if(doevGF3) then
call cpu_time(start_GF3)
call evGF3(maxSCF_GF,thresh_GF,n_diis_GF,renormGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO_basis,eHF)
call evGF3(maxSCF_GF,thresh_GF,n_diis_GF,renormGF,nBas,nC(1),nO(1),nV(1),nR(1),ERI_MO,eHF)
call cpu_time(end_GF3)
t_GF3 = end_GF3 - start_GF3
@ -631,7 +634,7 @@ program QuAcK
call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, &
singlet_manifold,triplet_manifold,linGW,eta, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0
@ -649,7 +652,7 @@ program QuAcK
call cpu_time(start_evGW)
call evGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, &
singlet_manifold,triplet_manifold,eta, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO_basis,PHF,cHF,eHF,eG0W0)
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_evGW)
t_evGW = end_evGW - start_evGW
@ -667,7 +670,7 @@ program QuAcK
call cpu_time(start_qsGW)
call qsGW(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA,G0W,GW0, &
singlet_manifold,triplet_manifold,eta, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO_basis,ERI_MO_basis,PHF,cHF,eHF)
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,S,X,T,V,Hc,ERI_AO,ERI_MO,PHF,cHF,eHF)
call cpu_time(end_qsGW)
t_qsGW = end_qsGW - start_qsGW
@ -687,8 +690,7 @@ program QuAcK
call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, &
singlet_manifold,triplet_manifold,linGW,eta, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF)
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0)
call cpu_time(end_G0T0)
t_G0T0 = end_G0T0 - start_G0T0
@ -705,7 +707,7 @@ program QuAcK
call cpu_time(start_evGT)
call evGT(maxSCF_GW,thresh_GW,n_diis_GW,doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_manifold, &
eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF,eG0T0)
eta,nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO,eHF,eG0T0)
call cpu_time(end_evGT)
t_evGT = end_evGT - start_evGT
@ -789,6 +791,64 @@ program QuAcK
end if
!------------------------------------------------------------------------
! Range-separeted GT/GW
!------------------------------------------------------------------------
if(doGTGW) then
! Read and transform long-range two-electron integrals
allocate(ERI_ERF_AO(nBas,nBas,nBas,nBas),ERI_ERF_MO(nBas,nBas,nBas,nBas))
call read_LR(nBas,ERI_ERF_AO)
call cpu_time(start_AOtoMO)
write(*,*)
write(*,*) 'AO to MO transformation for long-range ERIs... Please be patient'
write(*,*)
call AOtoMO_integral_transform(nBas,cHF,ERI_ERF_AO,ERI_ERF_MO)
call cpu_time(end_AOtoMO)
deallocate(ERI_ERF_AO)
t_AOtoMO = end_AOtoMO - start_AOtoMO
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds'
write(*,*)
! Long-range G0W0 calculation
call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA, &
singlet_manifold,triplet_manifold,linGW,eta, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0W0 = ',t_G0W0,' seconds'
write(*,*)
! Short-range G0T0 calculation
ERI_ERF_MO(:,:,:,:) = ERI_MO(:,:,:,:) - ERI_ERF_MO(:,:,:,:)
call cpu_time(start_G0T0)
call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, &
singlet_manifold,triplet_manifold,linGW,eta, &
nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_ERF_MO,eHF,eG0T0)
call cpu_time(end_G0T0)
t_G0T0 = end_G0T0 - start_G0T0
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds'
write(*,*)
call matout(nBas,1,(eG0W0+eG0T0-eHF(:,1))*HaToeV)
end if
!------------------------------------------------------------------------
! Basis set correction
!------------------------------------------------------------------------
@ -799,7 +859,7 @@ program QuAcK
call cpu_time(start_Bas)
call basis_correction(nBas,nO,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
ERI_MO_basis,eHF,cHF,PHF,eG0W0)
ERI_MO,eHF,cHF,PHF,eG0W0)
call cpu_time(end_Bas)
t_Bas = end_Bas - start_Bas

2
src/QuAcK/linear_response_C_pp.f90
View File

@ -25,7 +25,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp)
! Define the chemical potential
eF = e(nO) + e(nO+1)
eF = 0d0
! eF = 0d0
! Build C matrix for the singlet manifold

2
src/QuAcK/linear_response_D_pp.f90
View File

@ -25,7 +25,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp)
! Define the chemical potential
eF = e(nO) + e(nO+1)
eF = 0d0
! eF = 0d0
! Build the D matrix for the singlet manifold

6
src/QuAcK/linear_response_pp.f90
View File

@ -121,10 +121,16 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z)
! allocate(order(nOO+nVV))
! call quick_sort(Omega(:),order(:),nOO+nVV)
! call matout(nOO+nVV,1,Omega(:)*HaToeV)
! Split the various quantities in p-p and h-h parts
call sort_ppRPA(ortho_eigvec,nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:))
! call matout(32,1,(Omega1(:) - Omega1(1))*HaToeV)
! Compute the RPA correlation energy
EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )

3
src/QuAcK/ppRPA.f90
View File

@ -56,6 +56,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER
nOO = nO*(nO+1)/2
nVV = nV*(nV+1)/2
! print*,'nOO = ',nOO
! print*,'nVV = ',nVV
! Memory allocation
allocate(Omega1(nVV,nspin),X1(nVV,nVV,nspin),Y1(nOO,nVV,nspin), &

56
src/utils/read_LR.f90 Normal file
View File

@ -0,0 +1,56 @@
subroutine read_LR(nBas,G)
! Read the long-range two-electron integrals from files
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
! Local variables
integer :: mu,nu,la,si
double precision :: ERI
double precision :: lambda
! Output variables
double precision,intent(out) :: G(nBas,nBas,nBas,nBas)
! Open file with integrals
lambda = 1d0
print*, 'Scaling integrals by ',lambda
open(unit=11,file='int/ERI_lr.dat')
! Read two-electron integrals
G(:,:,:,:) = 0d0
do
read(11,*,end=11) mu,nu,la,si,ERI
ERI = lambda*ERI
! <12|34>
G(mu,nu,la,si) = ERI
! <32|14>
G(la,nu,mu,si) = ERI
! <14|32>
G(mu,si,la,nu) = ERI
! <34|12>
G(la,si,mu,nu) = ERI
! <41|23>
G(si,mu,nu,la) = ERI
! <23|41>
G(nu,la,si,mu) = ERI
! <21|43>
G(nu,mu,si,la) = ERI
! <43|21>
G(si,la,nu,mu) = ERI
enddo
11 close(unit=11)
end subroutine read_LR

6
src/utils/read_integrals.f90
View File

@ -19,7 +19,11 @@ subroutine read_integrals(nEl,nBas,S,T,V,Hc,G)
! Output variables
double precision,intent(out) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),G(nBas,nBas,nBas,nBas)
double precision,intent(out) :: S(nBas,nBas)
double precision,intent(out) :: T(nBas,nBas)
double precision,intent(out) :: V(nBas,nBas)
double precision,intent(out) :: Hc(nBas,nBas)
double precision,intent(out) :: G(nBas,nBas,nBas,nBas)
! Open file with integrals

89
src/utils/sort.f90 Normal file
View File

@ -0,0 +1,89 @@
subroutine quick_sort(x,iorder,isize)
implicit none
integer,intent(in) :: isize
double precision,intent(inout) :: x(isize)
integer,intent(inout) :: iorder(isize)
call rec_quicksort(x,iorder,isize,1,isize,1)
end subroutine quick_sort
recursive subroutine rec_quicksort(x,iorder,isize,first,last,level)
implicit none
integer, intent(in) :: isize, first, last, level
integer,intent(inout) :: iorder(isize)
double precision, intent(inout):: x(isize)
double precision :: c, tmp
integer :: itmp
integer :: i, j
c = x( shiftr(first+last,1) )
i = first
j = last
do
do while (x(i) < c)
i=i+1
end do
do while (c < x(j))
j=j-1
end do
if (i >= j) then
exit
endif
tmp = x(i)
x(i) = x(j)
x(j) = tmp
itmp = iorder(i)
iorder(i) = iorder(j)
iorder(j) = itmp
i=i+1
j=j-1
enddo
if ( ((i-first <= 10000).and.(last-j <= 10000)).or.(level<=0) ) then
if (first < i-1) then
call rec_quicksort(x, iorder, isize, first, i-1,level/2)
endif
if (j+1 < last) then
call rec_quicksort(x, iorder, isize, j+1, last,level/2)
endif
else
if (first < i-1) then
call rec_quicksort(x, iorder, isize, first, i-1,level/2)
endif
if (j+1 < last) then
call rec_quicksort(x, iorder, isize, j+1, last,level/2)
endif
endif
end subroutine rec_quicksort
subroutine set_order(x,iorder,isize,jsize)
implicit none
integer :: isize,jsize
double precision :: x(isize,jsize)
double precision,allocatable :: xtmp(:,:)
integer :: iorder(*)
integer :: i,j
allocate(xtmp(isize,jsize))
do i=1,isize
do j=1,jsize
xtmp(i,j) = x(i,iorder(j))
end do
end do
do i=1,isize
do j=1,jsize
x(i,j) = xtmp(i,j)
end do
end do
deallocate(xtmp)
end subroutine set_order