gdirect; gthresh,energy=1.d-8,orbital=1.d-8,grid=1.d-8
symmetry,nosym
orient,noorient
geometry={
Ne,,0.0,0.0,0.0
Ar,,0.0,0.0,6.5}
basis={
set,orbital; default,avtz         !for orbitals
set,jkfit;   default,avtz/jkfit   !for JK integrals
set,mp2fit;  default,avtz/mp2fit  !for E2disp/E2exch-disp
set,dflhf;   default,avtz/jkfit   !for LHF
}

!=========delta(HF) contribution for higher order interaction terms====
ca=2101.2; cb=2102.2 !sapt files

!dimer
{df-hf,basis=jkfit,locorb=0}
edm=energy

!monomer A
dummy,ar
{df-hf,basis=jkfit,locorb=0; save,$ca}
ema=energy
{sapt;monomerA}

!monomer B
dummy,ne
{df-hf,basis=jkfit,locorb=0; save,$cb}
emb=energy
{sapt;monomerB}

!interaction contributions
{sapt,Sinf=0,SAPT_LEVEL=2;intermol,ca=$ca,cb=$cb,icpks=1,fitlevel=3
dfit,basis_coul=jkfit,basis_exch=jkfit,cfit_scf=3}

!calculate high-order terms by subtracting 1st+2nd order energies
eint_hf=(edm-ema-emb)*1000 mH
delta_hf=eint_hf-e1pol-e1ex-e2ind-e2exind

!=========DFT-SAPT at second order intermol. perturbation theory====
ca=2103.2; cb=2104.2 !sapt files;

!shifts for asymptotic correction to xc potential
eps_homo_pbe0_ar=-0.440936      !HOMO(Ar)/PBE0 functional
eps_homo_pbe0_ne=-0.589207      !HOMO(Ne)/PBE0
ip_ar=0.5792                    !exp. ionisation potential
ip_ne=0.7925                    !exp. ionisation potential
shift_ar=ip_ar+eps_homo_pbe0_ar !shift for bulk xc potential (Ar)
shift_ne=ip_ne+eps_homo_pbe0_ne !shift for bulk xc potential (Ne)

!monomer A, perform LPBE0AC calculation
dummy,ar
{df-ks,pbex,pw91c,lhf; dftfac,0.75,1.0,0.25; asymp,shift_ne; save,$ca}
{sapt;monomerA}

!monomer B, perform LPBE0AC calculation
dummy,ne
{df-ks,pbex,pw91c,lhf; dftfac,0.75,1.0,0.25; start,atdens; asymp,shift_ar; save,$cb}
{sapt;monomerB}

!interaction contributions
{sapt,Sinf=0,SAPT_LEVEL=3;intermol,ca=$ca,cb=$cb,icpks=0,fitlevel=3,nlexfac=0.0
dfit,basis_coul=jkfit,basis_exch=jkfit,cfit_scf=3}

!add high-order approximation to obtain the total interaction energy
eint_dftsapt=e12tot+delta_hf