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

UG0F2

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-03-08 17:00:05 +01:00
parent 2d826d0f1e
commit 72470287f7
12 changed files with 284 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
input/methods
View File

@ -6,12 +6,12 @@
F F F F F F F F
# drCCD rCCD lCCD pCCD # drCCD rCCD lCCD pCCD
F F F F F F F F
# CIS* CIS(D) CID CISD # CIS* CIS(D) CID CISD FCI
F F F F F F F F F
# RPA* RPAx* ppRPA # RPA* RPAx* ppRPA
F F F F F F
# G0F2 evGF2 qsGF2* G0F3 evGF3 # G0F2* evGF2 qsGF2* G0F3 evGF3
F F T F F T F F F F
# G0W0* evGW* qsGW* # G0W0* evGW* qsGW*
F F F F F F
# G0T0 evGT qsGT # G0T0 evGT qsGT

2
input/options
View File

@ -7,7 +7,7 @@
# spin: TDA singlet triplet spin_conserved spin_flip # spin: TDA singlet triplet spin_conserved spin_flip
F T T T T F T T T T
# GF: maxSCF thresh DIIS n_diis lin eta renorm # GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.001 3 256 0.00001 T 5 T 0.001 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
256 0.00001 T 5 T 0.0 F F F F F 256 0.00001 T 5 T 0.0 F F F F F
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS

2
mol/butadiene.xyz
View File

@ -1,5 +1,5 @@
10 10
Butadiene,^1A_g,CC3,aug-cc-pVTZ
C 0.60673471 0.00000000 0.39936380 C 0.60673471 0.00000000 0.39936380
C -0.60673471 0.00000000 -0.39936380 C -0.60673471 0.00000000 -0.39936380
C 1.84223863 0.00000000 -0.11897388 C 1.84223863 0.00000000 -0.11897388

31
src/CI/FCI.f90 Normal file
View File

@ -0,0 +1,31 @@
subroutine FCI(nBas,nC,nO,nV,nR,ERI,e)
! Perform a full configuration interaction calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nC,nO,nV,nR
double precision,intent(in) :: e(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
! Hello world
write(*,*)
write(*,*)'**********************************'
write(*,*)'| Full Configuration Interaction |'
write(*,*)'**********************************'
write(*,*)
! Form FCI vector
! Form FCI matrix
! Diagonalize FCI matrix
end subroutine FCI

122
src/GF/UG0F2.f90 Normal file
View File

@ -0,0 +1,122 @@
subroutine UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linearize,eta,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
S,ERI,ERI_aaaa,ERI_aabb,ERI_bbbb,dipole_int_aa,dipole_int_bb,eHF)
! Perform unrestricted G0W0 calculation
implicit none
include 'parameters.h'
include 'quadrature.h'
! Input variables
logical,intent(in) :: BSE
logical,intent(in) :: TDA
logical,intent(in) :: dBSE
logical,intent(in) :: dTDA
logical,intent(in) :: evDyn
logical,intent(in) :: spin_conserved
logical,intent(in) :: spin_flip
logical,intent(in) :: linearize
double precision,intent(in) :: eta
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
integer,intent(in) :: nS(nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart)
double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart)
double precision,intent(in) :: eHF(nBas,nspin)
! Local variables
integer :: is
integer :: ispin
double precision :: Ec(nsp)
double precision :: EcBSE(nspin)
double precision,allocatable :: SigC(:,:)
double precision,allocatable :: Z(:,:)
integer :: nS_aa,nS_bb,nS_sc
double precision,allocatable :: eGF2lin(:,:)
double precision,allocatable :: eGF2(:,:)
! Output variables
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| One-shot G0F2 calculation |'
write(*,*)'| *** Unrestricted version *** |'
write(*,*)'************************************************'
write(*,*)
! TDA
if(TDA) then
write(*,*) 'Tamm-Dancoff approximation activated!'
write(*,*)
end if
! Memory allocation
nS_aa = nS(1)
nS_bb = nS(2)
nS_sc = nS_aa + nS_bb
allocate(SigC(nBas,nspin),Z(nBas,nspin),eGF2(nBas,nspin),eGF2lin(nBas,nspin))
!---------------------!
! Compute self-energy !
!---------------------!
call unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aaaa,ERI_aabb,ERI_bbbb,eHF,eGF2,SigC,Z)
!-----------------------------------!
! Solve the quasi-particle equation !
!-----------------------------------!
eGF2lin(:,:) = eHF(:,:) + Z(:,:)*SigC(:,:)
if(linearize) then
write(*,*) ' *** Quasiparticle energies obtained by linearization *** '
write(*,*)
eGF2(:,:) = eGF2lin(:,:)
else
! Find graphical solution of the QP equation
print*,'!!! Graphical solution NYI for UG0F2 !!!'
stop
end if
! Compute MP2 correlation energy
call UMP2(nBas,nC,nO,nV,nR,ERI_aaaa,ERI_aabb,ERI_bbbb,ENuc,EUHF,eHF,Ec)
! Dump results
call print_UG0F2(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec)
! Perform BSE calculation
if(BSE) then
print*,'!!! BSE2 NYI for UG0F2 !!!'
end if
end subroutine UG0F2

73
src/GF/print_UG0F2.f90 Normal file
View File

@ -0,0 +1,73 @@
subroutine print_UG0F2(nBas,nO,eHF,ENuc,EUHF,SigC,Z,eGF2,Ec)
! Print one-electron energies and other stuff for G0W0
implicit none
include 'parameters.h'
integer,intent(in) :: nBas
integer,intent(in) :: nO(nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: EUHF
double precision,intent(in) :: Ec(nsp)
double precision,intent(in) :: eHF(nBas,nspin)
double precision,intent(in) :: SigC(nBas,nspin)
double precision,intent(in) :: Z(nBas,nspin)
double precision,intent(in) :: eGF2(nBas,nspin)
integer :: p
integer :: ispin
double precision :: HOMO(nspin)
double precision :: LUMO(nspin)
double precision :: Gap(nspin)
! HOMO and LUMO
do ispin=1,nspin
if(nO(ispin) > 0) then
HOMO(ispin) = eGF2(nO(ispin),ispin)
LUMO(ispin) = eGF2(nO(ispin)+1,ispin)
Gap(ispin) = LUMO(ispin) - HOMO(ispin)
else
HOMO(ispin) = 0d0
LUMO(ispin) = eGF2(1,ispin)
Gap(ispin) = 0d0
end if
end do
! Dump results
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,*)' Unrestricted one-shot G0F2 calculation (eV)'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(A1,A3,A1,A30,A1,A30,A1,A30,A1,A30,A1)') &
'|',' ','|','e_HF ','|','Sig_c ','|','Z ','|','e_QP ','|'
write(*,'(A1,A3,A1,2A15,A1,2A15,A1,2A15,A1,2A15,A1)') &
'|','#','|','up ','dw ','|','up ','dw ','|','up ','dw ','|','up ','dw ','|'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
do p=1,nBas
write(*,'(A1,I3,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1,2F15.6,A1)') &
'|',p,'|',eHF(p,1)*HaToeV,eHF(p,2)*HaToeV,'|',SigC(p,1)*HaToeV,SigC(p,2)*HaToeV,'|', &
Z(p,1),Z(p,2),'|',eGF2(p,1)*HaToeV,eGF2(p,2)*HaToeV,'|'
enddo
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') 'UG0F2 HOMO energy:',maxval(HOMO(:))*HaToeV,' eV'
write(*,'(2X,A30,F15.6,A3)') 'UG0F2 LUMO energy:',minval(LUMO(:))*HaToeV,' eV'
write(*,'(2X,A30,F15.6,A3)') 'UG0F2 HOMO-LUMO gap :',(minval(LUMO(:))-maxval(HOMO(:)))*HaToeV,' eV'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,'(2X,A30,F15.6,A3)') ' UG0F2 total energy :',ENuc + EUHF + sum(Ec(:)),' au'
write(*,'(2X,A30,F15.6,A3)') ' UG0F2 correlation energy:',sum(Ec(:)),' au'
write(*,*)'-------------------------------------------------------------------------------&
-------------------------------------------------'
write(*,*)
end subroutine print_UG0F2

8
src/GF/qsUGF2.f90
View File

@ -306,4 +306,12 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,
deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,error,error_diis,F_diis) deallocate(cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,error,error_diis,F_diis)
! Perform BSE calculation
if(BSE) then
print*,'!!! BSE2 NYI for qsUGF2 !!!'
end if
end subroutine qsUGF2 end subroutine qsUGF2

2
src/GF/unrestricted_self_energy_GF2.f90
View File

@ -35,6 +35,8 @@ subroutine unrestricted_self_energy_GF2(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,ERI_b
! Compute self-energy | ! Compute self-energy |
!---------------------! !---------------------!
SigC(:,:,:) = 0d0
!----------------! !----------------!
! Spin-up sector ! Spin-up sector
!----------------! !----------------!

2
src/GF/unrestricted_self_energy_GF2_diag.f90
View File

@ -35,6 +35,8 @@ subroutine unrestricted_self_energy_GF2_diag(nBas,nC,nO,nV,nR,eta,ERI_aa,ERI_ab,
! Compute self-energy | ! Compute self-energy |
!---------------------! !---------------------!
SigC(:,:) = 0d0
!----------------! !----------------!
! Spin-up sector ! Spin-up sector
!----------------! !----------------!

2
src/MBPT/qsGW.f90
View File

@ -207,6 +207,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,SOSE
F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + SigCp(:,:)
call matout(nBas,nBAs,SigCp)
! Compute commutator and convergence criteria ! Compute commutator and convergence criteria
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)

36
src/QuAcK/QuAcK.f90
View File

@ -11,7 +11,7 @@ program QuAcK
logical :: doMP2,doMP3,doMP2F12 logical :: doMP2,doMP3,doMP2F12
logical :: doCCD,doDCD,doCCSD,doCCSDT logical :: doCCD,doDCD,doCCSD,doCCSDT
logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD
logical :: doCIS,doCIS_D,doCID,doCISD logical :: doCIS,doCIS_D,doCID,doCISD,doFCI
logical :: doRPA,doRPAx,doppRPA logical :: doRPA,doRPAx,doppRPA
logical :: doADC logical :: doADC
logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
@ -89,6 +89,7 @@ program QuAcK
double precision :: start_CIS ,end_CIS ,t_CIS double precision :: start_CIS ,end_CIS ,t_CIS
double precision :: start_CID ,end_CID ,t_CID double precision :: start_CID ,end_CID ,t_CID
double precision :: start_CISD ,end_CISD ,t_CISD double precision :: start_CISD ,end_CISD ,t_CISD
double precision :: start_FCI ,end_FCI ,t_FCI
double precision :: start_RPA ,end_RPA ,t_RPA double precision :: start_RPA ,end_RPA ,t_RPA
double precision :: start_RPAx ,end_RPAx ,t_RPAx double precision :: start_RPAx ,end_RPAx ,t_RPAx
double precision :: start_ppRPA ,end_ppRPA ,t_ppRPA double precision :: start_ppRPA ,end_ppRPA ,t_ppRPA
@ -162,7 +163,7 @@ program QuAcK
doMP2,doMP3,doMP2F12, & doMP2,doMP3,doMP2F12, &
doCCD,doDCD,doCCSD,doCCSDT, & doCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_lCCD,do_pCCD, & do_drCCD,do_rCCD,do_lCCD,do_pCCD, &
doCIS,doCIS_D,doCID,doCISD, & doCIS,doCIS_D,doCID,doCISD,doFCI, &
doRPA,doRPAx,doppRPA, & doRPA,doRPAx,doppRPA, &
doG0F2,doevGF2,doqsGF2, & doG0F2,doevGF2,doqsGF2, &
doG0F3,doevGF3, & doG0F3,doevGF3, &
@ -812,8 +813,19 @@ program QuAcK
if(doG0F2) then if(doG0F2) then
call cpu_time(start_GF2) call cpu_time(start_GF2)
call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, &
eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) if(unrestricted) then
call UG0F2(BSE,TDA,dBSE,dTDA,evDyn,spin_conserved,spin_flip,linGF,eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,EUHF, &
S,ERI_AO,ERI_MO_aaaa,ERI_MO_aabb,ERI_MO_bbbb,dipole_int_aa,dipole_int_bb,eHF)
else
call G0F2(BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet,linGF, &
eta_GF,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF)
end if
call cpu_time(end_GF2) call cpu_time(end_GF2)
t_GF2 = end_GF2 - start_GF2 t_GF2 = end_GF2 - start_GF2
@ -1176,6 +1188,22 @@ program QuAcK
end if end if
!------------------------------------------------------------------------
! Compute FCI
!------------------------------------------------------------------------
if(doFCI) then
call cpu_time(start_FCI)
call FCI(nBas,nC,nO,nV,nR,ERI_MO,eHF)
call cpu_time(end_FCI)
t_FCI = end_FCI - start_FCI
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for FCI = ',t_FCI,' seconds'
write(*,*)
end if
!------------------------------------------------------------------------ !------------------------------------------------------------------------
! End of QuAcK ! End of QuAcK
!------------------------------------------------------------------------ !------------------------------------------------------------------------

8
src/QuAcK/read_methods.f90
View File

@ -2,7 +2,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doMP2,doMP3,doMP2F12, & doMP2,doMP3,doMP2F12, &
doCCD,doDCD,doCCSD,doCCSDT, & doCCD,doDCD,doCCSD,doCCSDT, &
do_drCCD,do_rCCD,do_lCCD,do_pCCD, & do_drCCD,do_rCCD,do_lCCD,do_pCCD, &
doCIS,doCIS_D,doCID,doCISD, & doCIS,doCIS_D,doCID,doCISD,doFCI, &
doRPA,doRPAx,doppRPA, & doRPA,doRPAx,doppRPA, &
doG0F2,doevGF2,doqsGF2, & doG0F2,doevGF2,doqsGF2, &
doG0F3,doevGF3, & doG0F3,doevGF3, &
@ -20,7 +20,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
logical,intent(out) :: doMP2,doMP3,doMP2F12 logical,intent(out) :: doMP2,doMP3,doMP2F12
logical,intent(out) :: doCCD,doDCD,doCCSD,doCCSDT logical,intent(out) :: doCCD,doDCD,doCCSD,doCCSDT
logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD
logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD,doFCI
logical,intent(out) :: doRPA,doRPAx,doppRPA logical,intent(out) :: doRPA,doRPAx,doppRPA
logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3 logical,intent(out) :: doG0F2,doevGF2,doqsGF2,doG0F3,doevGF3
logical,intent(out) :: doG0W0,doevGW,doqsGW logical,intent(out) :: doG0W0,doevGW,doqsGW
@ -60,6 +60,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
doCIS_D = .false. doCIS_D = .false.
doCID = .false. doCID = .false.
doCISD = .false. doCISD = .false.
doFCI = .false.
doRPA = .false. doRPA = .false.
doRPAx = .false. doRPAx = .false.
@ -118,11 +119,12 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, &
! Read excited state methods ! Read excited state methods
read(1,*) read(1,*)
read(1,*) answer1,answer2,answer3,answer4 read(1,*) answer1,answer2,answer3,answer4,answer5
if(answer1 == 'T') doCIS = .true. if(answer1 == 'T') doCIS = .true.
if(answer2 == 'T') doCIS_D = .true. if(answer2 == 'T') doCIS_D = .true.
if(answer3 == 'T') doCID = .true. if(answer3 == 'T') doCID = .true.
if(answer4 == 'T') doCISD = .true. if(answer4 == 'T') doCISD = .true.
if(answer5 == 'T') doFCI = .true.
if(doCIS_D) doCIS = .true. if(doCIS_D) doCIS = .true.
read(1,*) read(1,*)