PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-03 20:53:53 +01:00
Code Releases Activity

fix ACFDT for Tmatrix

Browse Source
This commit is contained in:
Pierre-Francois Loos 2022-02-03 10:05:58 +01:00
parent 03f41ef20c
commit 21c1ceb975
13 changed files with 68 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -3,7 +3,7 @@
# MP2* MP3 MP2-F12
F F F
# CCD pCCD DCD CCSD CCSD(T)
F F F F F
F F F T F
# drCCD rCCD crCCD lCCD
F F F F
# CIS* CIS(D) CID CISD FCI

6
input/options
View File

@ -1,5 +1,5 @@
# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability
256 0.0000001 T 5 2 1 T T F
256 0.0000001 T 5 2 1 T 1.0 F
# MP:
# CC: maxSCF thresh DIIS n_diis
@ -13,8 +13,8 @@
# GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg
256 0.00001 T 5 T 0.0 F F
# ACFDT: AC Kx XBS
F F F
F T T
# BSE: BSE dBSE dTDA evDyn
T T T F
T F T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

2
mol/h2.xyz
View File

@ -1,4 +1,4 @@
2
H 0. 0. 0.
H 0. 0. 2.000000
H 0. 0. 0.741

8
src/GT/G0T0.f90
View File

@ -238,11 +238,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
write(*,*)'-------------------------------------------------------------------------------'
write(*,*)
! Free memory
deallocate(Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s, &
Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t)
! Compute the BSE correlation energy via the adiabatic connection
if(doACFDT) then
@ -260,7 +255,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing
end if
call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI_MO,eHF,eG0T0,EcAC)
nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, &
Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eHF,eG0T0,EcAC)
if(exchange_kernel) then

3
src/GT/evGT.f90
View File

@ -265,7 +265,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, &
end if
call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI_MO,eGT,eGT,EcAC)
nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, &
Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC)
if(exchange_kernel) then

3
src/GT/qsGT.f90
View File

@ -389,7 +389,8 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T
end if
call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI_MO,eGT,eGT,EcAC)
nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, &
Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

4
src/HF/RHF.f90
View File

@ -12,7 +12,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
double precision,intent(in) :: thresh
logical,intent(in) :: level_shift
double precision,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nO
@ -135,7 +135,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc
! Level-shifting
if(level_shift .and. Conv > thresh) call level_shifting(nBas,nO,S,c,F)
if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,c,F)
! Diagonalize Fock matrix

6
src/HF/UHF.f90
View File

@ -12,7 +12,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
logical,intent(in) :: level_shift
double precision,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
@ -187,10 +187,10 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,
! Level-shifting
if(level_shift .and. Conv > thresh) then
if(level_shift > 0d0 .and. Conv > thresh) then
do ispin=1,nspin
call level_shifting(nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin))
call level_shifting(level_shift,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin))
end do
end if

4
src/QuAcK/QuAcK.f90
View File

@ -111,8 +111,8 @@ program QuAcK
double precision :: start_Bas ,end_Bas ,t_Bas
integer :: maxSCF_HF,n_diis_HF
double precision :: thresh_HF
logical :: DIIS_HF,guess_type,ortho_type,mix,level_shift
double precision :: thresh_HF,level_shift
logical :: DIIS_HF,guess_type,ortho_type,mix
integer :: maxSCF_CC,n_diis_CC
double precision :: thresh_CC

9
src/QuAcK/read_options.f90
View File

@ -22,7 +22,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
integer,intent(out) :: guess_type
integer,intent(out) :: ortho_type
logical,intent(out) :: mix
logical,intent(out) :: level_shift
double precision,intent(out) :: level_shift
logical,intent(out) :: dostab
integer,intent(out) :: maxSCF_CC
@ -101,16 +101,15 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t
guess_type = 1
ortho_type = 1
mix = .false.
level_shift = .false.
level_shift = 0d0
dostab = .false.
read(1,*)
read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,answer3,answer4
read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,level_shift,answer3
if(answer1 == 'T') DIIS_HF = .true.
if(answer2 == 'T') mix = .true.
if(answer3 == 'T') level_shift = .true.
if(answer4 == 'T') dostab = .true.
if(answer3 == 'T') dostab = .true.
if(.not.DIIS_HF) n_diis_HF = 1

90
src/RPA/ACFDT_Tmatrix.f90
View File

@ -1,5 +1,6 @@
subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, &
ERI,eT,eGT,EcAC)
nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, &
Omega2t,X2t,Y2t,rho1t,rho2t,ERI,eT,eGT,EcAC)
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix
@ -26,6 +27,28 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
integer,intent(in) :: nR
integer,intent(in) :: nS
integer,intent(in) :: nOOs
integer,intent(in) :: nOOt
integer,intent(in) :: nVVs
integer,intent(in) :: nVVt
double precision,intent(in) :: Omega1s(nVVs)
double precision,intent(in) :: X1s(nVVs,nVVs)
double precision,intent(in) :: Y1s(nOOs,nVVs)
double precision,intent(in) :: Omega2s(nOOs)
double precision,intent(in) :: X2s(nVVs,nOOs)
double precision,intent(in) :: Y2s(nOOs,nOOs)
double precision,intent(in) :: rho1s(nBas,nBas,nVVs)
double precision,intent(in) :: rho2s(nBas,nBas,nOOs)
double precision,intent(in) :: Omega1t(nVVt)
double precision,intent(in) :: X1t(nVVt,nVVt)
double precision,intent(in) :: Y1t(nOOt,nVVt)
double precision,intent(in) :: Omega2t(nOOt)
double precision,intent(in) :: X2t(nVVt,nOOt)
double precision,intent(in) :: Y2t(nOOt,nOOt)
double precision,intent(in) :: rho1t(nBas,nBas,nVVt)
double precision,intent(in) :: rho2t(nBas,nBas,nOOt)
double precision,intent(in) :: eT(nBas)
double precision,intent(in) :: eGT(nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
@ -39,46 +62,23 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
double precision :: lambda
double precision,allocatable :: Ec(:,:)
integer :: nOOs,nOOt
integer :: nVVs,nVVt
double precision :: EcRPA(nspin)
double precision,allocatable :: TA(:,:)
double precision,allocatable :: TB(:,:)
double precision,allocatable :: TAs(:,:)
double precision,allocatable :: TBs(:,:)
double precision,allocatable :: TAt(:,:)
double precision,allocatable :: TBt(:,:)
double precision,allocatable :: Omega(:,:)
double precision,allocatable :: XpY(:,:,:)
double precision,allocatable :: XmY(:,:,:)
double precision,allocatable :: Omega1s(:),Omega1t(:)
double precision,allocatable :: X1s(:,:),X1t(:,:)
double precision,allocatable :: Y1s(:,:),Y1t(:,:)
double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:)
double precision,allocatable :: Omega2s(:),Omega2t(:)
double precision,allocatable :: X2s(:,:),X2t(:,:)
double precision,allocatable :: Y2s(:,:),Y2t(:,:)
double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:)
! Output variables
double precision,intent(out) :: EcAC(nspin)
! Useful quantities
nOOs = nO*nO
nVVs = nV*nV
nOOt = nO*(nO-1)/2
nVVt = nV*(nV-1)/2
! Memory allocation
allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), &
Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), &
rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), &
Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), &
Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), &
rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt))
allocate(TA(nS,nS),TB(nS,nS),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
allocate(TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS), &
Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin))
allocate(Ec(nAC,nspin))
! Antisymmetrized kernel version
@ -113,11 +113,6 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
lambda = rAC(iAC)
! Initialize T matrix
TA(:,:) = 0d0
TB(:,:) = 0d0
if(doXBS) then
isp_T = 1
@ -128,8 +123,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TAs)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TBs)
isp_T = 2
iblock = 4
@ -139,12 +134,12 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TAt)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TBt)
end if
call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, &
call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt+TAs,TBt+TBs, &
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))
@ -183,11 +178,6 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
lambda = rAC(iAC)
! Initialize T matrix
TA(:,:) = 0d0
TB(:,:) = 0d0
if(doXBS) then
isp_T = 1
@ -198,8 +188,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TAs)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TBs)
isp_T = 2
iblock = 4
@ -209,12 +199,12 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple
call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB)
call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TAt)
if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TBt)
end if
call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, &
call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt-TAs,TBt-TBs, &
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))

6
src/eDFT/UKS.f90
View File

@ -20,7 +20,7 @@ subroutine UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,t
integer,intent(in) :: max_diis
integer,intent(in) :: guess_type
logical,intent(in) :: mix
logical,intent(in) :: level_shift
double precision,intent(in) :: level_shift
double precision,intent(in) :: thresh
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
@ -275,10 +275,10 @@ subroutine UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,t
! Level-shifting
if(level_shift .and. Conv > thresh) then
if(level_shift > 0d0 .and. Conv > thresh) then
do ispin=1,nspin
call level_shifting(nBas,maxval(nO(ispin,:)),S,c,F(:,:,ispin))
call level_shifting(level_shift,nBas,maxval(nO(ispin,:)),S,c,F(:,:,ispin))
end do
end if

6
src/utils/level_shifting.f90
View File

@ -1,4 +1,4 @@
subroutine level_shifting(nBas,nO,S,c,F)
subroutine level_shifting(level_shift,nBas,nO,S,c,F)
! Perform level-shifting on the Fock matrix
@ -6,6 +6,7 @@ subroutine level_shifting(nBas,nO,S,c,F)
! Input variables
double precision,intent(in) :: level_shift
integer,intent(in) :: nBas
integer,intent(in) :: nO
double precision,intent(in) :: S(nBas,nBas)
@ -15,7 +16,6 @@ subroutine level_shifting(nBas,nO,S,c,F)
double precision,allocatable :: F_MO(:,:)
double precision,allocatable :: Sc(:,:)
double precision,parameter :: LS = 0.1d0
integer :: a
@ -28,7 +28,7 @@ subroutine level_shifting(nBas,nO,S,c,F)
F_MO(:,:) = matmul(transpose(c),matmul(F,c))
do a=nO+1,nBas
F_MO(a,a) = F_MO(a,a) + LS
F_MO(a,a) = F_MO(a,a) + level_shift
end do
Sc(:,:) = matmul(S,c)