====== Symmetry-adapted intermolecular perturbation theory ====== ===== Introduction ===== The SAPT (symmetry-adapted intermolecular perturbation theory) program calculates the total interaction energy between closed-shell molecules as a sum of individual first and second order interaction terms, namely electrostatic $E_{\text{pol}}^{(1)}$, induction $E_{\text{ind}}^{(2)}$ and dispersion $E_{\text{disp}}^{(2)}$ accompanied by their respective exchange counterparts ($E_{\text{exch}}^{(1)}$, $E_{\text{exch-ind}}^{(2)}$ and $E_{\text{exch-disp}}^{(2)}$). The latter ones arise due to electron exchange between the monomers when the molecules are close to each other and are sometimes denoted as Pauli repulsion. Since all above terms are accessible through density matrices and static and dynamic density-density response functions of the monomers, in principle (see section [[symmetry-adapted intermolecular perturbation theory#high order terms|high order terms]]) no calculation of the dimer wave function is required. Therefore SAPT is free from the basis set superposition error which occurs in the supermolecular approach. References: **General Symmetry-adapted perturbation theory and many-body SAPT:**\\ $[1]$ B. Jeziorski, R. Moszynski and K. Szalewicz, Chem. Rev. **94**, 1887. (1994).\\ **DFT-SAPT:**\\ $[2]$ G. Jansen and A. Heßelmann, J. Phys. Chem. A **105**, 646 (2001).\\ $[3]$ A. Heßelmann and G. Jansen, [[https://dx.doi.org/10.1016/S0009-2614(02)00538-9|Chem. Phys. Lett.]] **357**, 464 (2002).\\ $[4]$ A. Heßelmann and G. Jansen, [[https://dx.doi.org/10.1016/S0009-2614(02)01097-7|Chem. Phys. Lett.]] **362**, 319 (2002).\\ $[5]$ A. Heßelmann and G. Jansen, [[https://dx.doi.org/10.1016/S0009-2614(02)01796-7|Chem. Phys. Lett.]] **367**, 778 (2003).\\ $[6]$ A. Heßelmann and G. Jansen, [[https://dx.doi.org/10.1039/B310529F|Phys. Chem. Chem. Phys.]] **5**, 5010 (2003).\\ **Density fitting DFT-SAPT (DF-DFT-SAPT):**\\ $[7]$ A. Heßelmann, G. Jansen and M. Schütz, [[https://dx.doi.org/10.1063/1.1824898|J. Chem. Phys.]] **122**, 014103 (2005). **Density fitting DFT-SAPT for arbitrary monomer basis sets (DFSAPT):**\\ $[8]$ A. Heßelmann and T. Korona, [[https://dx.doi.org/10.1063/1.4893990|J. Chem. Phys.]] **141**, 094107 (2014). **DFT-SAPT employing exact-exchange response kernels (EXX-SAPT):**\\ $[9]$ A. Heßelmann, [[https://dx.doi.org/10.1021/acs.jctc.7b01233|J. Chem. Theory Comput.]] **14**, 1943 (2018) **Single-determinant-based SAPT without single-exchange approximation:**\\ $[10]$ R. Schäffer and G. Jansen, //Theor. Chem. Acc.// **131** (2012) 1235\\ $[11]$ R. Schäffer and G. Jansen, [[https://dx.doi.org/10.1080/00268976.2013.827253|Mol. Phys.]] **111**, 2570 (2013) **SAPT with regularised nuclear attraction integrals:**\\ $[12]$ A. J. Misquitta, [[https://dx.doi.org/10.1021/ct400704a|J. Chem. Theory Comput.]] **9**, 5313 (2013) (See also:\\ K. Szalewicz, K. Patkowski and B. Jeziorski, Struct. Bond **116**, 43 (2005)\\ K. Szalewicz, WIREs Comput. Mol. Sci. **2**, 254 (2011)\\ and references therein for a related approach to DFT-SAPT termed SAPT(DFT)) ===== First example ===== A typical input for SAPT has the following form: r=5.6 geometry={nosym; he1; he2,he1,r} basis=avqz !wf records ca=2101.2 cb=2102.2 !monomer A dummy,he2 {hf; save,$ca} sapt;monomerA !monomer B dummy,he1 {hf; start,atdens; save,$cb} sapt;monomerB !interaction contributions sapt;intermol,ca=$ca,cb=$cb Here the ''%%sapt;monomerA/B%%'' store some informations about the two monomers which are needed in the subsequent SAPT calculation invoked by ''%%sapt;intermol%%''. The individual interaction energy terms are stored (in millihartree) in distinct variables and may be collected in arrays for producing potential energy surfaces. For example the input geometry={nosym; he1; he2,he1,r} basis=avtz !wf records ca=2101.2 cb=2102.2 !distances dist=[4.5,5.0,5.5,5.6,6.0,6.5,7.0] do i=1,#dist r=dist(i) !monomer A dummy,he2 {hf; save,$ca} sapt;monomerA !monomer B dummy,he1 {hf; start,atdens; save,$cb} sapt;monomerB !interaction contributions sapt;intermol,ca=$ca,cb=$cb elst(i)=E1pol; exch(i)=E1ex ind(i)=E2ind; exind(i)=E2exind disp(i)=E2disp; exdisp(i)=E2exdisp etot(i)=E12tot data,truncate,$ca enddo {table,dist,elst,exch,ind,exind,disp,exdisp,etot ftyp,d,d,d,d,d,d,d,d,d plot} yields the plot {{he2.svg?500}} Currently SAPT only accepts single-determinant wave functions for the monomers, i.e. from Hartree-Fock or Kohn-Sham DFT (see next section) calculations. This means that if Hartree-Fock wave functions are used for monomer, the following quantity is obtained (zero in superscript denotes that no intramonomer correlation is accounted for) [1]. $$\begin{aligned} E_{\text{SAPT}}=E^{(10)}_{\text{pol}}+E^{(10)}_{\text{exch}}+ E^{(20)}_{\text{ind,resp.}}+ E^{(20)}_{\text{exch-ind,resp.}}+E^{(20)}_{\text{disp}}+ E^{(20)}_{\text{exch-disp}}\end{aligned}$$ No point group symmetry can be exploited in a SAPT calculation. ===== DFT-SAPT ===== It is of crucial importance to account for the //intra//molecular correlation effects of the individual SAPT terms since Hartree-Fock theory often yields poor first- and second-order electrostatic properties. While this can be done using many-body perturbation theory [1] (in a double perturbation theory ansatz) a more efficient way is to use static and time-dependent DFT theory. This variant of SAPT, termed as DFT-SAPT [2-6], has in contrast to Hartree-Fock-SAPT the appealing feature that the polarisation terms ($E_{\text{pol}}^{(1)}$, $E_{\text{ind}}^{(2)}$, $E_{\text{disp}}^{(2)}$) are potentially exact, i.e. they come out exactly if the exact exchange-correlation (xc) potential and the exact (frequency-dependent) xc response kernel of the monomers were known. On the other hand, this does not hold for the exchange terms since Kohn-Sham theory can at best give a good approximation to the exact density matrix of a many-body system. It has been shown [6] that this is indeed the the case and therefore DFT-SAPT has the potential to produce highly accurate interaction energies comparable to high-level supermolecular many-body perturbation or coupled cluster theory. However, in order to achieve this accuracy, it is of crucial importance to correct the wrong asymptotic behaviour of the xc potential in current DFT functionals [3-5]. This can be done by using e.g.: {ks,lda; asymp,} which activates the gradient-regulated asymptotic correction approach of Grüning //et al.// ([[https://dx.doi.org/10.1063/1.1327260|J. Chem. Phys.]] **114**, 652 (2001)) for the respective monomer calculation. The user has to supply a ''shift'' parameter ($\Delta_{\text{xc}}$) for the bulk potential which should approximate the difference between the HOMO energy ($\varepsilon_{\text{HOMO}}$) obtained from the respective standard Kohn-Sham calculation and the (negative) ionisation potential of the monomer ($\mathrm{IP}$): $$\Delta_{\text{xc}}=\varepsilon_{\text{HOMO}}-(-\mathrm{IP})$$ This method accounts for the derivative discontinuity of the exact xc-potential and that is missing in approximate ones. Note that this needs to be done only once for each system. (See also section [[symmetry-adapted intermolecular perturbation theory#DFT-SAPT calculation of the NeAr dimer using the delta-HF correction|DFT-SAPT calculation of the NeAr dimer using the delta-HF correction]] for an explicit example). It is also possible to use an alternative asymptotic correction method described in Ref. [9] by using the ''EXX'' program (see section [[the density functional program#optimised Effective Potential (OEP) method|optimised Effective Potential (OEP) method]]). Here, the derivative discontinuity shift is determined through an interpolation between the HOMO energy of the EXX method and the HOMO energy of a corresponding hybrid functional. Currently, this technique should only be used in conjuction with hybrid methods employing 25% exact exchange. Concerning the more technical parameters in the DFT monomer calculations it is recommended to use lower convergence thresholds and larger intergration grids compared to standard Kohn-Sham calculations. ===== High order terms ===== It has been found that third and higher-order terms become quite important if one or both monomers are polar. As no higher than second-order terms are currently implemented in SAPT, one may use a non-correlated estimation of those terms by using supermolecular Hartree-Fock (see e.g. [7]). This can be done by adapting the following template: !dimer hf edm=energy !monomer A dummy, {hf; save,$ca} ema=energy sapt;monomerA !monomer B dummy, {hf; start,atdens; save,$cb} emb=energy sapt;monomerB !interaction contributions sapt,sapt_level=2;intermol,ca=$ca,cb=$cb esup=(edm-ema-emb)*1000. mH dHF=esup-e1pol-e1ex-e2ind-e2exind which stores the resulting $\delta$(HF) term in ''dHF''. ===== Exchange interaction energies with infinite expansion in the intermolecular overlap ===== The exchange interaction energies are commonly expanded in powers of the intermolecular overlap $S$: $$\begin{aligned} E_{\text{exch}}^{(1)}(S^\infty)&=&E_{\text{exch}}^{(1)}(S^2)+\Delta E_{\text{exch}}^{(1)}(S^4) +\Delta E_{\text{exch}}^{(1)}(S^6)\dots \\ E_{\text{exch-ind}}^{(2)}(S^\infty)&=&E_{\text{exch-ind}}^{(2)}(S^2)+\Delta E_{\text{exch-ind}}^{(2)}(S^4) +\Delta E_{\text{exch-ind}}^{(2)}(S^6)+\dots \\ E_{\text{exch-disp}}^{(2)}(S^\infty)&=&E_{\text{exch-disp}}^{(2)}(S^2)+\Delta E_{\text{exch-disp}}^{(2)}(S^4) +\Delta E_{\text{exch-disp}}^{(2)}(S^6)+\dots \\\end{aligned}$$ where $E_{\text{exch}}^{(1)}(S^2)$ corresponds to the exchange interaction energy in which the intermolecular overlap occurs only in single and quadratic powers and $\Delta E_{\text{exch}}^{(1)}(S^4)$ collects all terms up to the fourth power in $S$ excluding the terms in $E_{\text{exch}}^{(1)}(S^2)$. The above expansions are commonly truncated after the first terms, leading to the so-called $S^2$ approximation. Accordingly, for very small distances of the monomers, the exchange interaction energies become inaccurate and eventually the $S^2$ approximation may even break down, see Refs. [10,11]. While this problem may be circumvented by adding the $\delta$(HF) term to the SAPT interaction energy (see section [[symmetry-adapted intermolecular perturbation theory#high order terms|high order terms]]), a derivation of the infinite expansions of the first-order exchange energy was given by Jeziorski, Bulski and Piela (//Int. J. Quantum Chem.// **10** (1976) 281). More recently, formulas for the second-order exchange induction and dispersion energies were derived by Schäffer and Jansen, see Refs. [10,11]. In Molpro, the calculation of the exchange interaction energies in the $S^2$ approximation or with an infinite expansion in $S$ ($S^\infty$) can be controlled through the switch sapt,sinf=<0,1,2> where ''sinf=0'' will calculate the exchange interaction energies within the $S^2$ approximation and ''sinf=1'' calls the routines to calculate all exchange interaction terms with infinite overlap expansion ($S^\infty$). Using ''sinf=2'' will calculate both the $S^2$ and the $S^\infty$ contributions. Note that the first-order exchange energy is always calculated in both variants. ===== Using exact exchange Kohn-Sham response kernels ===== The exact exchange Kohn-Sham response kernel can be used in the computations of the second order (exchange-)induction and (exchange-)dispersion contributions, see Ref. [9]. In order to do so, the $\langle i|v_x^{\text{NL}}-v_x|j\rangle$ matrix elements have to be precomputed for each monomer and supplied to the SAPT program. The following input template shows how to compute DFT-SAPT interaction energy contributions using the (A)TDEXX response kernels and the LPBE0AC xc potentials (cf section [[the density functional program#optimised Effective Potential (OEP) method|optimised Effective Potential (OEP) method]]): dummy, !monomer A {ks,lda; save,2100.2; start,atdens} {exxhyb,orb=2100.2,oep=7,homo=1,charge=1,maxit=30,scfthr=1d-9,add=1,xfac=0.25,shift=1 func,pbex,pbec; dftfac,0.75,1.0} {oep,orb=2101.2,vxdiff=5100.2,gamma=3d-4,dfit=1,homo=1} {sapt;monomerA} dummy, !monomer B {ks,lda; save,2200.2; start,atdens} {exxhyb,orb=2200.2,oep=7,homo=1,charge=1,maxit=30,scfthr=1d-9,add=1,xfac=0.25,shift=1 func,pbex,pbec; dftfac,0.75,1.0} {oep,orb=2201.2,vxdiff=5200.2,gamma=3d-4,dfit=1,homo=1} {sapt;monomerB} sapt;intermol,ca=2101.2,cb=2201.2,icpks=0,vxdiffA=5100.2,vxdiffB=5200.2,\ nlexfac=1.0,tdexx= The input parameter ''TDEXXMODE'' takes the values 1 or 2. With ''TDEXXMODE=1'' the full nonadiabatic TDEXX response kernel and with ''TDEXXMODE=2'' only the adiabatic TDEXX response kernel is used in the calculation of the response matrices of the monomers. Notice also from the above template that it is important to set ''nlexfac=1.0'' which sets the scaling factor of exact exchange in the SAPT program. If ''nlexfac'' is set to a value lower than one, a hybrid ALDA-TDEXX kernel will be used in the calculation. All reported calculations using this program should cite Ref. [9]. Note that unlike in the conventional hybrid-DFT case, DFT-SAPT calculations can be performed also using density fitting techniques for (A)TDEXX kernels with Molpro. A corresponding efficient DF implementation which works in conjunction with (A)TDEXX kernels is the DFSAPT program, see section [[symmetry-adapted intermolecular perturbation theory#DFSAPT: a density-fitting DFT-SAPT program for arbitrary monomer basis sets|DFSAPT: a density-fitting DFT-SAPT program for arbitrary monomer basis sets]]. ===== SAPT with regularised nuclear attraction integrals ===== The SAPT expansion of the intermolecular interaction energy does not yield the charge-transfer (CT) interaction energy as a separate term, rather it is included in the sum of the second-order induction and exchange-induction energies: $$E_{\text{int}}^{\text{CT}} \subset E_{\text{ind}}^{(2)}+E_{\text{exch-ind}}^{(2)}$$ Due to Misquitta, see Ref. [12], the CT interaction energy can, however, be ’estimated’ by using a regularised nuclear attraction potential for the monomers in order to calculate (exchange-)induction energies excluding the (short-ranged) CT interaction contribution. The CT interaction energy can then be obtained by calculating the difference between the full (exchange-)induction interaction energy and the (exchange-)induction interaction energy as obtained by the regularised SAPT method: $$E_{\text{int}}^{\text{CT}} \approx E_{\text{ind}}^{(2)}+E_{\text{exch-ind}}^{(2)} -E_{\text{ind}}^{(2)}[regul]-E_{\text{exch-ind}}^{(2)}[regul]$$ For obtaining regularised nuclear attraction integrals Misquitta proposes the form [12]: $$V_{\text{regul}}=\frac{1}{r}\Big(1-e^{-\alpha r^2}\Big)$$ with $\alpha=3.0$. In Molpro, regularised SAPT calculations can be performed by calling the ''vreg'' module ’before’ saving the monomer geometries with ''%%{sapt; monomerA,B}%%''. An input template for the calculation of the CT interaction energy is given by: dummy, {ks,; orbital,2100.2} {sapt;monomerA} dummy, {ks,; start,atdens; orbital,2200.2} {sapt;monomerB} sapt,sapt_level=2;intermol,ca=2100.2,cb=2200.2,icpks=0 eIND1=e2ind+e2exind dummy, vreg,alpha=3d0,nq=1 {sapt;monomerA} dummy, vreg,alpha=3d0,nq=1 {sapt;monomerB} sapt,sapt_level=2;intermol,ca=2100.2,cb=2200.2,icpks=0 eIND2=e2ind+e2exind eint_CT=eIND1-eIND2 The ''vreg'' module accepts the following options: * **''ALPHA''** Exponent of the exponential function of the damping factor. * **''NQ''** Global exponent of the damping factor (type: integer). * **''SAVE''** Record (''record.file'') in which to save the integrals. If not used, the standard record that holds the nuclear attraction integrals is overwritten. The regularised SAPT can also be used in order to resolve the breakdown of the conventional SAPT expansion for the interaction between heavy atoms, see, e.g., Liu //et al.//, [[https://dx.doi.org/10.1021/ct200243s|J. Chem. Theory Comput.]] **7**, 2399 (2011). A regularised SAPT method was, e.g., used in T. Clark, A. Heßelmann, [[https://dx.doi.org/10.1039/C8CP03079K|Phys. Chem. Chem. Phys.]] **20**, 22849 (2018) for calculating the bonding between carbon tetrahalides and halide anions. Note that in this case the addition of the $\delta$(HF) term (see section [[symmetry-adapted intermolecular perturbation theory#high order terms|high order terms]]) for taking into account short-ranged polarisation effects is crucial. ===== Density fitting ===== In order to be able to study interactions between extended monomers one can use density fitting to approximate the integrals in SAPT [7]. For this one may use the input: {sapt;intermol,ca=$ca,cb=$cb,fitlevel=3 dfit,basis_coul=jkfit,basis_exch=jkfit,basis_mp2=mp2fit,cfit_scf=3} with in the basis section defined ''jkfit'' and ''mp2fit'' fitting basis sets (see section [[density fitting]]). Currently only ALDA or (A)TDEXX xc-kernels are implemented for the case ''SAPT_LEVEL''=3 and ''SAPT_FITLEVEL''=3. This means that a corresponding SAPT calculation would be uncompatible with (conventional) hybrid-DFT monomer orbitals/orbital energies. Therefore it is recommended to use nonhybrid functionals in the case the dispersion/exchange-dispersion energy terms are requested in a DF-DFT-SAPT run. Another possibility is to localise the xc-potential via, e.g., the OEP method (see section [[the density functional program#exact exchange Kohn-Sham methods|exact exchange Kohn-Sham methods]] and the example in section [[symmetry-adapted intermolecular perturbation theory#DF-DFT-SAPT calculation of the NeAr dimer using the delta-HF correction|DF-DFT-SAPT calculation of the NeAr dimer using the delta-HF correction]]). New density fitting routines for the calculation of $E^{(2)}_{\text{IND}}(S^\infty)$ and $E^{(2)}_{\text{DISP}}(S^\infty)$ can be activated with the switch sapt,df2=.true.,... For performing density-fitting DFT-SAPT calculations with arbitrary monomer basis sets see section [[symmetry-adapted intermolecular perturbation theory#DFSAPT: a density-fitting DFT-SAPT program for arbitrary monomer basis sets|DFSAPT: a density-fitting DFT-SAPT program for arbitrary monomer basis sets]]. ===== SAPT with ECP’s ===== If effective core potentials (ECP’s) are used in the monomer calculations, it is important to add the $\delta$(HF) term to the SAPT interaction energy (see K. Patkowski, K. Szalewicz, [[https://dx.doi.org/10.1063/1.2784391|J. Chem. Phys.]] **127**, 164103 (2007)). For examples for the calculation of $\delta$(HF) see sections [[symmetry-adapted intermolecular perturbation theory#high order terms|high order terms]] and [[symmetry-adapted intermolecular perturbation theory#examples|examples]]. ===== Examples ===== ==== HF-SAPT calculation of the H2O dimer using the delta-HF correction ==== gthresh,energy=1.d-8,orbital=1.d-8,grid=1.d-8 symmetry,nosym orient,noorient GEOMETRY={ 1,O1,,0.00000000,0.00000000,0.00000000 2,H1,,0.00000000,0.00000000,1.83606000 3,H2,,1.77604000,0.00000000,-0.4656040 4,O2,,-0.6605540,0.00000000,5.54064000 5,H3,,-1.6582100,-1.4536300,6.05324000 6,H4,,-1.6582100,1.45363000,6.05324000 } basis=avdz !sapt files ca=2101.2 cb=2102.2 !dimer hf edm=energy !monomer A dummy,o2,h3,h4 {hf; save,$ca} ema=energy {sapt;monomerA} !monomer B dummy,o1,h1,h2 {hf; start,atdens; save,$cb} emb=energy {sapt;monomerB} !interaction contributions {sapt,SAPT_LEVEL=3;intermol,ca=$ca,cb=$cb,icpks=1} !HF supermolecular interaction energy and delta(HF) contribution eint_hf=(edm-ema-emb)*1000 mH delta_hf=eint_hf-e1pol-e1ex-e2ind-e2exind !add E2disp + E2exch-disp to HF interaction energy eint_sapt=eint_hf+e2disp+e2exdisp ==== DFT-SAPT calculation of the NeAr dimer using the delta-HF correction ==== 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=avtz !=========delta(HF) contribution for higher order interaction terms==== !sapt files ca=2101.2 cb=2102.2 !dimer hf edm=energy !monomer A dummy,ar {hf; save,$ca} ema=energy {sapt;monomerA} !monomer B dummy,ne {hf; start,atdens; save,$cb} emb=energy {sapt;monomerB} !interaction contributions {sapt,SAPT_LEVEL=2;intermol,ca=$ca,cb=$cb,icpks=1} !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==== !sapt files ca=2103.2 cb=2104.2 !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 dummy,ar {ks,pbe0; asymp,shift_ne; save,$ca} {sapt;monomerA} !monomer B dummy,ne {ks,pbe0; start,atdens; asymp,shift_ar; save,$cb} {sapt;monomerB} !interaction contributions {sapt;intermol,ca=$ca,cb=$cb,icpks=0} !add high-order approximation to obtain the total interaction energy eint_dftsapt=e12tot+delta_hf ==== DF-DFT-SAPT calculation of the NeAr dimer using the delta-HF correction ==== 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 ===== Options ===== * **''SAPT_LEVEL''** Set to 1 for first-order terms ($E_{\text{pol}}^{(1)}$ and $E_{\text{exch}}^{(1)}$), to 2 for additional second order (exchange-)induction terms ($E_{\text{ind}}^{(2)}$ and $E_{\text{exch-ind}}^{(2)}$) and 3 for all first- and second-order terms (including then also $E_{\text{disp}}^{(2)}$ and $E_{\text{exch-disp}}^{(2)}$) (default 3) * **''SINF''** Set to 0 to calculate exchange interaction energies within the $S^2$ approximation and to 1 for calculating all exchange terms with an infinite overlap expansion (default 1). * **''SAPT_FITLEVEL''** Level of density fitting approximations in SAPT which can have values 0 to 3 (default 0) * **''SAPT_ICPKS''** Switch between iterative (=1) and non-iterative (=0) solution of coupled-perturbed Kohn-Sham equations (default 0) * **''SAPT_CPKSTHR''** Threshold for density matrix convergency in the coupled-perturbed Kohn-Sham program (default 1.d-6). * **''SAPT_CPKMAXIT''** Maximum number of iterations in the coupled-perturbed Kohn-Sham program (default 50). * **''SAPT_FROZENA''** Number of frozen electrons in the response calculations for monomer A (default 0) * **''SAPT_FROZENB''** see abobe The following parameters are of importance if ''SAPT_FITLEVEL''$>$0: * **''SAPT_NFRQ_DISP''** Number of frequencies for the Casimir-Polder integration (default 10) * **''SAPT_NORM_DISP''** Norm for the density fitting which can be either ''COULOMB'' or ''NATURAL'' (default ''COULOMB'') * **''SAPT_DISP_N4''** Can speedup the calculation of the dispersion energy by $N^4$ scaling (default 1) * **''THR_XCKERN''** Density threshold for the xc kernel matrix elements (default 1.d-8) * **''FIT_XCKERN''** Fit both sides of the xc kernel (default 0) * **''SAPT_DISK''** If 0 write all dimer amplitudes to file, if 1 write 3-index response propagators to file and if 2 write 3-index response propagators compressed to file. The latter two variants save disk space but need more CPU time to compute $E^{(2)}_{\text{exch-disp}}$ (default 0) * **''COMPRESS_THR''** If ''SAPT_DISK''=2 this value determines the compression cutoff (default 1d-12) * **''UNCOUPLED''** If ''SAPT_DISK''$>$0 calculate also uncoupled (exchange-)dispersion energies (default false) * **''THRAO''** Threshold for AO 3-index integrals (default 1.d-12) * **''THRMO''** Threshold for MO 3-index integrals (default 1.d-8) * **''THROV''** Threshold for AO 2-index integrals (default 1.d-10) * **''THRPROD''** Product threshold for first half transformation (default 1.d-8) * **''THRSW''** Threshold for Schwarz screening (default 1.d-5) * **''C6''** Calculate dispersion coefficients for the two monomers (Note that the full dimer basis set is used in each case and that a closer distance of the monomers can perturb the result). * **''XCKERN_NBLOCK''** number of grid points treated together as a block for (aux$|$f$_{\text{xc}}|$occ$\times$virt) integrals (default 128) * **''CFAC''** factor for VWN correlation in ALDA xc-kernel (default 1.d0) * **''SAPT_QUAD''** quadrature type in frequency integration (0: Chebyshev, 1: Gauss-Legendre (default)) * **''W0''** parameter in Gauss-Legendre quadrature (default 2.d0) * **''GRIDTHR''** Threshold for grid accuracy in some routines which generate xc-kernel integrals. * **''THRDUM''** Threshold for density calculated at dummy function gridpoints (which can serve as a singularity correction for the calculation of the XC kernel matrix) (default: ''1d-12''). * **''LSING''** Singularity correction for XC kernel matrix (default: ''.false.'') * **''DF2''** Logical switch for new DF routines. Note that with ''DF2=.true.'' the exchange interaction contributions are only calculated without the single-exchange approximation. (default: ''.false.'') * **''DF2MOD''** Controls calculation mode of $E^{(2)}_{\text{DISP}}$: with ''DF2MOD=0'' all 3idx integrals are calculated and stored on disk in advance. In case of ''DF2MOD=1'' this is done in two subsequent steps at the expense, however, of the fact that then the amplitudes have to be constructed twice. (default: ''0'') * **''X0THR''** SVD threshold for inversion of $\chi_0$ if EXX kernel used (default: ''1d-8'') * **''NAUX_MAX''** Maximal length of auxiliary function batches in EXX kernel routine (default: ''50''). * **''PVVDISK''** Write $(P|ab)$ ints to disk in EXX kernel calculation. This can speed up the calculation on the expense of the disk space requirements (scaling as: $N_{\text{aux}}^A$ $\times N_{\text{vir}}^A$ $\times (N_{\text{vir}}^A+1)/2$ $+N_{\text{aux}}^B$ $\times N_{\text{vir}}^B$ $\times (N_{\text{vir}}^B+1)/2$ (default=''.false.''). The last threshold values for the 2- and 3-index integrals should not be set higher in density fitting calculations as this can cause lower accuracies in the interaction terms. In addition SAPT knows the following subcommands: * **''MONOMERA''** Stores informations (like number of electrons, etc.) about previous monomer A calculation * **''MONOMERB''** See above * **''INTERMOL''** Starts the SAPT calculation ''INTERMOL'' may have the following subkeywords: * **''CA''** Record number of wave function for monomer A (always needed) * **''CB''** Record number of wave function for monomer B (always needed) * **''SAPTLEVEL''** See above * **''FITLEVEL''** See above * **''ICPKS''** See above * **''FROZA''** See above * **''FROZB''** See above * **''NLEXFAC''** Amount of nonlocal exact exchange in hybrid DFT-SAPT calculations * **''CPKSTHR''** Threshold for density matrix convergency in the coupled-perturbed Kohn-Sham program. * **''CPKSMAXIT''** Maximum number of iterations in the coupled-perturbed Kohn-Sham program. ===== DFSAPT: a density-fitting DFT-SAPT program for arbitrary monomer basis sets ===== DFSAPT is an alternative density-fitting DFT-SAPT program which works with arbitrary monomer basis sets and which can be used for extended dimer systems of about 800 electrons, see Ref.[8]. For a standard DFSAPT calculation using a dimer centered basis set the Molpro input file looks very similar to a corresponding DF-DFT-SAPT claculation, see [[symmetry-adapted intermolecular perturbation theory#DFSAPT calculation of the water dimer using a DCBS basis set|DFSAPT calculation of the water dimer using a DCBS basis set]]. There are, however, two important points to be taken care of: * ****important #1**** Set ''oldnorm=1'' before any call to the DFSAPT program. * ****important #2**** Add numbers to atom labels in the geometry blocks. These have to have a subsequent ordering, see examples in [[symmetry-adapted intermolecular perturbation theory#DFSAPT calculation of the water dimer using a DCBS basis set|DFSAPT calculation of the water dimer using a DCBS basis set]] to [[symmetry-adapted intermolecular perturbation theory#DFSAPT calculation of the water dimer with additional d(HF) correction|DFSAPT calculation of the water dimer with additional d(HF) correction]]. A corresponding calculation, in which the water monomers are calculated in their own local basis set (monomer centered basis set calculation) is presented in [[symmetry-adapted intermolecular perturbation theory#DFSAPT calculation of the water dimer using an MCBS basis set|DFSAPT calculation of the water dimer using an MCBS basis set]]. Note, however, that very large basis sets are required to converge the individual interaction energy contributions. In order to improve the basis set convergence, additional midbond functions can be added as is exemplified in [[symmetry-adapted intermolecular perturbation theory#DFSAPT calculation of the water dimer using an MC+BS basis set|DFSAPT calculation of the water dimer using an MC+BS basis set]]. Here the def2-TZVP basis functions for the Neon atom were positioned in between the O$\cdots$H hydrogen bridge. A complete (DCBS) DFSAPT calculation with additional $\delta$(HF) correction and using the PBE0AC(BRJ) xc potential (with automatic asymptotic correction) is shown in section [[symmetry-adapted intermolecular perturbation theory#DFSAPT calculation of the water dimer with additional d(HF) correction|DFSAPT calculation of the water dimer with additional d(HF) correction]]. Note that in contrast to the DFT-SAPT/DF-DFT-SAPT programs the exchange-dispersion energy is not calculated on the coupled level. Instead, the coupled $E^{(2)}_{exch-disp}$ energies are estimated by scaling the uncoupled ones with a factor which has been obtained from a linear fit to exact results, see Ref.[8]. The following list shows all options which can be used with the DFSAPT program: * **''E1P''** Calculate $E^{(1)}_{pol}$ (default: ''1''). * **''E1X''** Calculate $E^{(1)}_{exch}$ (default: ''1''). * **''E2I''** Calculate $E^{(2)}_{ind}$ (default: ''0''). * **''E2XI''** Calculate $E^{(2)}_{exch-ind}$ (default: ''0''). * **''E2D''** Calculate $E^{(2)}_{disp}$ (default: ''0''). * **''E2XD''** Calculate $E^{(2)}_{exch-disp}$ (default: ''0''). * **''ALL''** Calculate all SAPT contributions (up to 2nd order) (default: ''0''). * **''AUXBAS''** Auxiliary basis set (default: ''MP2FIT''). * **''NFREQ''** Number of frequencies for computing the monomer response matrices (default: ''10''). * **''GRIDTHR''** Threshold for grid accuracy (default: ''1d-10''). * **''THRKERN''** Threshold for density in xc kernel calculation (default: ''1d-12''). * **''PRINT''** Print flag (default: ''0''). * **''W0''** Value effects Gauss-Legendre quadrature for computing the response function (see Amos,Handy,Knowles,Rice,Stone,//J. Phys.Chem.// **89** (1985) 2186) (default: ''2d0'') * **''FEXX''** Exact exchange factor in response calculations (default: ''0d0''). * **''FXC''** Exchange-correlation factor in response calculations (default: ''0d0''). * **''FCOUL''** Coulomb contribution factor in response calculations (default: ''1d0''). * **''XFAC''** Exchange factor in XC kernel (default: ''1d0''). * **''CFAC''** Correlation factor in XC kernel (default: ''1d0''). * **''AUXMAX''** Maximum block size of auxiliary function space for three-index integrals of the type (aux$|$orb.orb) (default: ''50''). * **''MOMAX''** Maximum block size of orbital function space for three-index integrals of the type (orb$|$orb.aux) (default: ''50''). * **''FXCSING''** Singularity correction for XC kernel matrix (default: ''0'') * **''THRDUM''** Threshold for density calculated at dummy function gridpoints (which can serve as a singularity correction for the calculation of the XC kernel matrix) (default: ''1d-4''). * **''COMP''** Write 3-index integrals compressed to disk (to save disk space) (default: ''0''). * **''THRC''** Threshold for 3-index integral calculations (integrals which are smaller than the given threshold are not written to disk) (default: ''1d-10''). * **''MINB''** Minimal bytelength for integral comression (default: ''1''). * **''MODE''** Switch for alternative DFSAPT driver which can use integral compression (to be evoked by ''MODE=2'') (default: ''1''). * **''FXCFIT''** Use density fitting for computation of XC contributions in CPKS calculation (default: ''0''). * **''CPKS''** Switch for different CPKS solvers (default: ''1''). * **''CDIRECT''** Use integral direct algorithm for computing the (exact) exchange contributions in CPKS calculations, otherwise intergrals of the type (aux$|$vir.vir) are written to disk (default: ''0''). * **''MAXIT''** Maximum number of iteractions in CPKS (default: ''60''). * **''THRCG''** Threshold for CG/GMRES solver in CPKS calculations (default: ''1d-8''). * **''NSYSOLVE''** Switch for BiCG (''NSYSOLVE=1'') or GMRES (''NSYSOLVE=2'') solver in CPKS calculations (default: ''2''). * **''NSUB''** Subspace size in GMRES solver (default: ''50''). * **''T678''** Switch for 3 different computation modes for $E^{(2)}_{exch-disp}$. ''T678=3'' is recommended for the calculation of large systems but requires more disk space to store intermediate quantities (default: ''1''). * **''TDEXX''** set to 1 for using the adiabatic TDEXX kernel and to 2 for using the nonadiabatic TDEXX kernel * **''TDEXXFAC''** prefactor for TDEXX response kernel * **''X0THR''** threshold for SVD inversion of the KS response matrix * **''KERNMOD''** if set to 1, DFSAPT computes (aux$|$virt$\times$virt) integrals for use in the TDEXX kernel calculation and stores them on disk. If set to 2, these integrals are directly contracted to reduce I/O operations. * **''RESTART''** when set to 1, integral records are not deleted and can be used for subsequent calculations. This is possible by using ''RESTART=2''. * **''PLUSOV''** take terms that contain $\langle occ|v_x^{\text{NL}}-v_x|virt\rangle$ matrix elements into account in the calculation of the kernel (normally very small) The DFSAPT program can also be used to compute the monomer response functions within the exact-exchange Kohn-Sham response (TDEXX) approach, see section [[symmetry-adapted intermolecular perturbation theory#using exact exchange Kohn-Sham response kernels|using exact exchange Kohn-Sham response kernels]]. As in the standard DFT-SAPT program, the $\langle i|v_x^{\text{NL}}-v_x|j\rangle$ have to be computed for each monomer using the OEP program (section [[the density functional program#optimised Effective Potential (OEP) method|optimised Effective Potential (OEP) method]]). An input template for this would be given by: dummy, !monomer A {df-ks,lda; save,2100.2; start,atdens} {exxhyb,orb=2100.2,oep=7,homo=1,charge=1,maxit=30,scfthr=1d-9,add=1,\ xfac=0.25,shift=1,dfit=1 func,pbex,pbec; dftfac,0.75,1.0} {oep,orb=2101.2,vxdiff=5100.2,gamma=3d-4,dfit=1,homo=1} {sapt;monomerA} dummy, !monomer B {ks,lda; save,2200.2; start,atdens} {exxhyb,orb=2200.2,oep=7,homo=1,charge=1,maxit=30,scfthr=1d-9,add=1,\ xfac=0.25,shift=1,dfit=1 func,pbex,pbec; dftfac,0.75,1.0} {oep,orb=2201.2,vxdiff=5200.2,gamma=3d-4,dfit=1,homo=1} {sapt;monomerB} dfsapt,all=1,fcoul=1.0,fxc=1.0,cpks=2,x0thr=1d-8,tdexx=2,nfreq=10, \ tdexxfac=1d0,kernmod=1,auxmax=100 which again presents an example where the orbitals of the monomers are calculated with the LPBE0AC xc potential. All publications resulting from DFSAPT calculations employing the (A)TDEXX response kernel should cite Ref. [9]. ==== DFSAPT calculation of the water dimer using a DCBS basis set ==== ***,h2o-dimer DCBS gthresh,energy=1.d-8,orbital=1.d-8,grid=1.d-8 oldnorm=1 ca=2101.2 cb=2102.2 basis={ set,orbital default,def2-tzvpp set,mp2fit default,def2-tzvpp/mp2fit} symmetry,nosym angstrom geometry={ 1,O1,, -1.551007, -0.114520, 0.000000 2,H2,, -1.934259, 0.762503, 0.000000 3,H3,, -0.599677, 0.040712, 0.000000 4,O4,, 1.350625, 0.111469, 0.000000 5,H5,, 1.680398, -0.373741, -0.758561 6,H6,, 1.680398, -0.373741, 0.758561} dummy,4,5,6 {ks,pbe; save,$ca} monomerA,orb=$ca,core=1 dummy,1,2,3 {ks,pbe; start,atdens; save,$cb} monomerB,orb=$cb,core=1 {grid,gridthr=1d-8} dfsapt,all=1,fcoul=1.0,fxc=1.0 ==== DFSAPT calculation of the water dimer using an MCBS basis set ==== ***,h2o-dimer MCBS gthresh,energy=1.d-8,orbital=1.d-8,grid=1.d-8 oldnorm=1 ca=2101.2 cb=2102.2 symmetry,nosym angstrom noorient geometry={ 1,O1,, -1.551007, -0.114520, 0.000000 2,H2,, -1.934259, 0.762503, 0.000000 3,H3,, -0.599677, 0.040712, 0.000000} basis={set,orbital; default,def2-tzvpp} {ks,pbe; save,$ca} monomerA,orb=$ca,core=1 symmetry,nosym angstrom noorient geometry={ 4,O4,, 1.350625, 0.111469, 0.000000 5,H5,, 1.680398, -0.373741, -0.758561 6,H6,, 1.680398, -0.373741, 0.758561} {ks,pbe; start,atdens; save,$cb} monomerB,orb=$cb,core=1 symmetry,nosym angstrom noorient geometry={ 1,O1,, -1.551007, -0.114520, 0.000000 2,H2,, -1.934259, 0.762503, 0.000000 3,H3,, -0.599677, 0.040712, 0.000000 4,O4,, 1.350625, 0.111469, 0.000000 5,H5,, 1.680398, -0.373741, -0.758561 6,H6,, 1.680398, -0.373741, 0.758561} basis={ set,orbital; default,def2-tzvpp set,jkfit; default,def2-tzvpp/jkfit set,mp2fit; default,def2-tzvpp/mp2fit} {grid,gridthr=1d-8} dfsapt,all=1,fcoul=1.0,fxc=1.0 ==== DFSAPT calculation of the water dimer using an MC+BS basis set ==== ***,h2o-dimer MC+BS gthresh,energy=1.d-8,orbital=1.d-8,grid=1.d-8 oldnorm=1 ca=2101.2 cb=2102.2 basis={set,orbital; default,def2-tzvpp} symmetry,nosym angstrom noorient geometry={ 1,O1,, -1.551007, -0.114520, 0.000000 2,H2,, -1.934259, 0.762503, 0.000000 3,H3,, -0.599677, 0.040712, 0.000000 4,Ne7,, 0.375474, 0.076091, 0.000000} dummy,ne7 {ks,pbe; save,$ca} monomerA,orb=$ca,core=1 symmetry,nosym angstrom noorient geometry={ 4,O4,, 1.350625, 0.111469, 0.000000 5,H5,, 1.680398, -0.373741, -0.758561 6,H6,, 1.680398, -0.373741, 0.758561 7,Ne7,, 0.375474, 0.076091, 0.000000} dummy,ne7 {ks,pbe; start,atdens; save,$cb} monomerB,orb=$cb,core=1 symmetry,nosym angstrom noorient geometry={ 1,O1,, -1.551007, -0.114520, 0.000000 2,H2,, -1.934259, 0.762503, 0.000000 3,H3,, -0.599677, 0.040712, 0.000000 4,O4,, 1.350625, 0.111469, 0.000000 5,H5,, 1.680398, -0.373741, -0.758561 6,H6,, 1.680398, -0.373741, 0.758561 7,Ne7,, 0.375474, 0.076091, 0.000000} basis={ set,orbital; default,def2-tzvpp set,jkfit; default,def2-tzvpp/jkfit set,mp2fit; default,def2-tzvpp/mp2fit} dummy,ne7 {grid,gridthr=1d-8} dfsapt,all=1,fcoul=1.0,fxc=1.0 ==== DFSAPT calculation of the water dimer with additional d(HF) correction ==== gdirect gthresh,energy=1d-8,orbital=1d-7,grid=1d-8 symmetry,nosym angstrom geometry={ 1,O1,, -1.551007, -0.114520, 0.000000 2,H2,, -1.934259, 0.762503, 0.000000 3,H3,, -0.599677, 0.040712, 0.000000 4,O4,, 1.350625, 0.111469, 0.000000 5,H5,, 1.680398, -0.373741, -0.758561 6,H6,, 1.680398, -0.373741, 0.758561} basis={ set,orbital; default,avdz set,jkfit; default,avdz/jkfit set,mp2fit; default,avdz/mp2fit} oldnorm=1 ca=2110.2; cb=2120.2 nkerna=1; nkernb=1 {df-hf,basis=jkfit; save,2100.2} edm_hf=energy dummy,4,5,6 {df-hf,basis=jkfit; save,$ca; start,atdens} ema_hf=energy monomerA,orb=$ca,core=$nkerna dummy,1,2,3 {df-hf,basis=jkfit; save,$cb; start,atdens} emb_hf=energy monomerB,orb=$cb,core=$nkernb {grid,gridthr=1d-6} {dfsapt,e1p=1,e1x=1,e2i=1,e2xi=1,auxbas='jkfit',fexx=1.0,fcoul=1.0,mode=1,cdirect=1,thrcg=1d-6} {grid,gridthr=1d-8} eint_hf=(edm_hf-ema_hf-emb_hf)*1000e0 mH delta_hf=eint_hf-e12tot dummy,4,5,6 lxgamma=1.15 {df-ks,lxbecke,pbex,pbec,basis=jkfit; dftfac,0.25,0.75,1.0; asymp,0d0; save,$ca; start,atdens} monomerA,orb=$ca,core=$nkerna dummy,1,2,3 lxgamma=1.15 {df-ks,lxbecke,pbex,pbec,basis=jkfit; dftfac,0.25,0.75,1.0; asymp,0d0; save,$cb; start,atdens} monomerB,orb=$cb,core=$nkernb {grid,gridthr=1d-6} {dfsapt,all=1,auxbas='mp2fit',fxc=1.0,fcoul=1.0,mode=1,cpks=2,fxcfit=1,fxcsing=1,t678=3} eint_sapt=e12tot+delta_hf ===== SAPT(CCSD) ===== SAPT with monomers described on the CCSD level is available for small complexes. In SAPT(CCSD) monomer density matrices and density-density matrix response functions from expectation-value CCSD theory are utilized. A high cost of SAPT(CCSD) results from the necessity of calculation of the CCSD response functions. Cumulant contributions from two-electron density matrices for $E_{\text{exch}}^{(1)}$ and $E_{\text{exch-ind}}^{(2)}$ are also available. The calculations can be performed for dimer-centered or monomer-centered-plus basis sets (the latter for all components but $E^{(2)}_{\text{exch-disp}}$). $E^{(2)}_{\text{disp}}$ and $E^{(2)}_{\text{exch-disp}}$ are usually calculated from density-fitted response functions in order to reduce the CPU time (from ${\cal O}({\cal N}^8)$ to ${\cal O}({\cal N}^7)$, where ${\cal N}$ is the molecular size). Note that $E^{(2)}_{\text{disp}}$ and $E^{(2)}_{\text{exch-disp}}$ are calculated separately from other SAPT components. It is important to note that SAPT(CCSD) and DFT-SAPT from Chapter [[symmetry-adapted intermolecular perturbation theory#introduction|introduction]] are completely different codes. References: **Review of SAPT(CCSD):**\\ $[1]$ T. Korona, in //Recent Progress in Coupled Cluster Methods//, Eds. P. Čársky, J. Paldus, J. Pittner, Springer-Verlag (2010), //Coupled cluster treatment of intramonomer correlation effects in intermolecular interactions//, p. 267\\ $[2]$ T. Korona, M. Przybytek, B. Jeziorski, [[https://dx.doi.org/10.1080/00268970600673975|Mol. Phys.]] **104**, 2303 (2006)\\ $[3]$ T. Korona, B. Jeziorski, [[https://dx.doi.org/10.1063/1.2364489|J. Chem. Phys.]] **125**, 184109 (2006)\\ $[4]$ T. Korona, B. Jeziorski, [[https://dx.doi.org/10.1063/1.2889006|J. Chem. Phys.]] **128**, 144107 (2008)\\ $[5]$ T. Korona, [[https://dx.doi.org/10.1063/1.2933312|J. Chem. Phys.]] **128**, 224104 (2008)\\ $[6]$ T. Korona, [[https://dx.doi.org/10.1039/B804513E|Phys. Chem. Chem. Phys.]] **10**, 5698 (2008)\\ $[7]$ T. Korona, [[https://dx.doi.org/10.1039/B807329E|Phys. Chem. Chem. Phys.]] **10**, 6509 (2008)\\ $[8]$ T. Korona, [[https://dx.doi.org/10.1021/ct900232j|J. Chem. Theory Comput.]] **5**, 2663 (2009) ===== MP2-coupled (MP2C) ===== The standard second-order Møller-Plesset perturbation theory (MP2) method fails to describe long-range correlation energies accurately. The reason for this stems from the fact that the MP2 method does not take into account intramolecular correlation effects of two interacting monomers, i.e., describes dispersion interactions on an uncoupled Hartree-Fock level [1]. To remedy this, the supermolecular MP2 interaction energy can be corrected in a hybrid supermolecular-perturbation theory approach by subtracting the uncoupled HF (UCHF) dispersion energy and adding the dispersion energy from a more accurate response theory method instead. In the MP2C (MP2 coupled) method, the latter is computed using the time-dependent density-functional theory (TDDFT) approach employing static response functions from the EXX method and coupled response functions computed by using the ALDA exchange kernel: $$\Delta E_{\text{int}}^{\text{MP2C}}=\Delta E_{\text{int}}^{\text{MP2}} -E_{\text{disp}}^{(20)}[\mbox{UCHF}] +E_{\text{disp}}^{(2)}[\mbox{TDDFT}]$$ In Molpro, MP2C calculations can be performed with the aid of the TDDFT module that computes frequency-dependent coupled or uncoupled response functions. Using the LHF (Local Hartree-Fock) method (section [[the density functional program#local Hartree-Fock (LHF) method|local Hartree-Fock (LHF) method]]) an input for computing the MP2C interaction energy could read (using density-fitting within the HF/MP2 calculations and using 10 frequency integration points) basis={ set,orbital; default, set,jkfit; default,/jkfit set,mp2fit; default,/mp2fit set,dflhf; default,/jkfit set,tddft; default,/mp2fit} {df-hf,basis=jkfit} !Dimer {df-mp2} edm_mp2=energy dummy, !Monomer A {df-hf,basis=jkfit; save,} {df-mp2} ema_mp2=energy {dftresp; respmat,orb=,save=,order=0,nfreq=10,ncore=,init=1} dummy, !Monomer B {df-hf,basis=jkfit; save,} {df-mp2} ema_mp2=energy {dftresp; respmat,orb=,save=,order=0,nfreq=10,ncore=} {dftresp; disp} edisp_uhf=e2disp !uncoupled dispersion energy dummy, !Monomer A {cfit,basis=jkfit,locorb=0} dflhf,orb=,lhffit=1 {dftresp; respmat,orb=,save=,nfreq=10,ncore=,cfac=0.0,init=1} dummy, !Monomer B {cfit,basis=jkfit,locorb=0} dflhf,orb=,lhffit=1 {dftresp; respmat,orb=,save=,nfreq=10,ncore=,cfac=0.0} {dftresp; disp} edisp_lhf=e2disp !coupled dispersion energy esup_mp2disp=esup_mp2+(-edisp_uhf+edisp_lhf)*1000.0 mH show,esup_mp2 !MP2 interaction energy show,edisp_uhf !uncoupled dispersion energy show,edisp_lhf !coupled dispersion energy show,esup_mp2disp !MP2C interaction energy Notice that in the basis input section a number of basis sets have to be specified for the different parts in the calculation. Alternatively, the MP2C method can also be used in combination with the OEP method described in section [[the density functional program#optimised Effective Potential (OEP) method|optimised Effective Potential (OEP) method]]. In this case, the latter part where the coupled dispersion energy is calculated would have to be replaced by something like: dummy, !Monomer A {cfit,basis=jkfit,locorb=0} oep,orb=,gamma=3d-4,dfit=1 {dftresp; respmat,orb=,save=,nfreq=10,core=,cfac=0.0,init=1} dummy, !Monomer B {cfit,basis=jkfit,locorb=0} oep,orb=,gamma=3d-4,dfit=1 {dftresp; respmat,orb=,save=,nfreq=10,core=,cfac=0.0} The ''TDDFT'' program used to compute the response functions can have the following further options: * **ORB** record containing orbital coefficients (mandatory) * **SAVE** record number used to save the response matrix (in the auxiliary basis) for each frequency * **ORDER** can be used to denote the order of the electron-electron interaction that is taken into account in the calculation of the response matrix, so ''ORDER=0'' corresponds to the uncoupled case. If set to a negative value (default) the electron-electron interactions are taken into account through infinite order (fully coupled case, the default setting) * **NORM** may be used to change the fitting norm (default ’J’, should better not be changed) * **3IDX** a logical flag (choices: 0/1) to specify if the xc-kernel matrix should be computed by using 3-indexed (aux$|f_{\text{xc}}|$occ$\times$virt) integrals * **GRIDTHR** threshold for quadrature integration * **GRIDFREE** set to $\ne$ 0 if the xc kernel matrix shall be computed with a gridfree approach (see Ref. [2]). This is more efficient but can be less accurate and is only recommended if used in conjunction with large auxiliary basis sets. * **THRKERN** a threshold that is used to remove the singularities of the $(P|f_{\text{xc}}|Q)$ matrix, see Ref. [3] * **XFAC** scaling factor of ALDA exchange * **CFAC** scaling factor of ALDA correlation * **INIT** set to $\ne$ 0 in order to initialse an MP2C calculation. In subsequent calculations (where ''INIT''=0) the response matrices are stored for each individual call to ''TDDFT'' * **NCORE** number of core orbitals * **NBLOCK** num,ber of batch points in the numerical integration * **THRED** can be used to compress the size of the xc kernel values on the grid in order to increase the efficiency (experimental) * **FCX3MOD** switch between different computation modes for the calculation of the (aux$|f_{\text{xc}}|$occ$\times$virt) integrals * **FXC2MOD** switch between different computation modes for the calculation of the (aux$|f_{\text{xc}}|$aux) integrals * **THRDUM** when integrating over a dummy centre and when the electron density is lower than this value, skip the integration point (used to eliminate the singularities of the xc kernel * **SING** set to $\ne$ 0 in order to employ a singularity correction to the xc kernel matrix (is used by default) * **POL** print polarisabilities * **ESHIFT** parameter $\alpha$ in Eq. (5) in Ref. [3] * **FSCAL** parameter $\beta$ in Eq. (5) in Ref. [3] * **MATINV** flag for choosing matrix inversion routine to take inverse of the metric matrix (1: LU decomposition and 2: SVD decomposition). * **THRINV** threshold for SVD matrix inversion All publications resulting from use of this program to compute MP2C interaction energies should acknowledge Refs. [1] and [2]. If the singularity correction to the xc kernel matrix is used, additionally also Ref. [3] can be cited. Bibliography:\\ $[1]$ A. Heßelmann, [[https://dx.doi.org/10.1063/1.2905808|J. Chem. Phys.]] **128**, 144112 (2008)\\ $[2]$ M. Pitonak, A. Heßelmann, [[https://dx.doi.org/10.1021/ct9005882|J. Chem. Theory Comput.]] **6**, 168 (2010)\\ $[3]$ A. Heßelmann and T. Korona, [[https://dx.doi.org/10.1063/1.4893990|J. Chem. Phys.]] **141**, 094107 (2014).