16 SAPT(DFT)—SAPT based on coupled Kohn-Sham treatment of monomers

16.1 Introduction

SAPT(DFT) is an extension of the SAPT theory. In SAPT(DFT), the monomers are described in terms of Kohn-Sham (KS) orbitals and orbital energies as well as of TD-DFT response functions. A complete description of the theory together with references to the historical developments of the method, as well as numerical examples, can be found in Ref. 14.

SAPT(DFT) consists of two steps. First is the so-called SAPT(KS), where electrostatics, first-order exchange, induction, exchange-induction, dispersion, and exchange-dispersion are obtained by using SAPT terms of the zeroth order in W with SCF orbitals and orbital energies replaced by their KS equivalents. For meaningful results, SAPT(KS) requires asymptotically corrected (AC) Kohn-Sham calculations. SAPT(KS) does not reproduce dispersion correctly. This problem was found to be due to the use of a formula asymptotically related to uncoupled dynamic polarizabilities. Instead, the dispersion (and induction) energies should be calculated from frequency-dependent density susceptibility (FDDS) functions, also referred to as propagators, obtained from TD-DFT, i.e., at the coupled Kohn-Sham (CKS) level. The SAPT(KS) and CKS steps form the complete SAPT(DFT). The total SAPT(DFT) interaction energy (up to second order in V ) can be defined as [14]

ESAP T(DFT) = E (1)(KS )+ E (1) (KS ) + E(2)(CKS )+ ˜E(2) (CKS ) int elst exch ind exch-ind +E (d2is)p(CKS )+ E˜(2ex)ch-disp(CKS ), (14)
where the terms with the CKS label result from the coupled Kohn–Sham approach. The exch-ind(2)(CKS) and exch-disp(2)(CKS) are approximated by scaling their KS counterparts as described in Ref. 14. The method can be shown to be potentially exact for all major components of the interaction energy (asymptotically for exchange interactions) in the sense that these components would be exact if the DFT description of the monomers were exact.

The nominal scaling of SAPT(KS) is O(N5) and that of SAPT(DFT) is O(N6), in both cases significantly better than the O(N7) scaling of the regular SAPT. The scaling of SAPT(DFT) can be lowered to O(N5) by using density fitting [1214155762]. Furthermore, although the density-fitting transformation still scales as O(N5), it has a greatly reduced prefactor.

16.2 Installation and usage

SAPT(DFT) requires dalton 2.0 [37] or orca [53] for monomer DFT calculations (in the latter case, only the density-fitting version has been interfaced, because orca does not provide two-electron integrals needed for the standard version. See Sec. 16.4.2 for the detailed description). The users of SAPT(DFT) need to obtain a separate license and download the package from
http://www.kjemi.uio.no/software/dalton/dalton.html. Consult the dalton manual for installation requirements. Some of the SAPT(DFT) functionality is available via a patch dalton.diff distributed together with sapt2012. First, compile and install the standard dalton 2.0. After testing the installation, update the path to the dalton directory in the patchdalton script located in the main sapt2012 directory and run the script. If the patching is successful, compile the patched dalton. The recompilation should last significantly shorter than the original compilation. If you have not compiled sapt2012 with dalton support, set the DALTON variable of the Compall script to the proper path to the dalton executable, make sure the variable SAPTDFT=YES, and recompile sapt2012. If the compilation is error-free, the SAPT(DFT) code is ready to use (both the regular and density-fitted versions are created).

Before performing calculations for systems of interest to you, check the examples from examples/saptdft directory. The SAPT(DFT) run consists of the monomers calculation, transformation, the SAPT(KS) step, and the CKS step. Notice that some of the errors printed by Dalton are harmless in the SAPT(DFT) calculations. In particular,

  --- SEVERE ERROR, PROGRAM WILL BE ABORTED ---  
   SIRIUS NORMAL STOP AFTER ORBITAL ORTHONORMALIZATION.

is a normal stop after calculating the necessary one-electron integrals. Warning. Patched dalton should not be used for any ather purpose except for SCF/DFT calculations. The effect of the patch on other parts of the dalton code has not been tested.

The keywords discussed below are for both the regular and density-fitted versions of SAPT(DFT). Some additional keywords for the latter case are described in Sec. 16.4 . SAPT(DFT) requires SAPTKS=T keyword in the INPUTCOR namelist of the nameP.data file and proper keywords in the TRN namelist prepared like for the regular runs (see also Sec. 9.5 ). For CKS calculations that are a part of SAPT(DFT), a separate SAPTDFT namelist with keywords CKSDISP=T and CKSIND=T should also be present at the minimum. The complete set of SAPTDFT namelist keywords is:

  1. CKSDISPEdisp(2)(CKS)
  2. CKSINDEind(2)(CKS)
  3. MAXMEM — maximum memory in words to be used in the CKS calculation. The program will allocate necessary memory not greater than MAXMEM. If the required memory is greater than MAXMEM, the program will exit. If not set, up to 100 megawords will be used.
  4. IQUADTYP : The type of quadrature scheme to be used in performing the ω-integral in the Casimir-Polder dispersion formula:
    • IQUADTYP=1 sets the Gauss-Legendre quadrature with the transformation of the integral variable ω = ω0(1+t) (1-t). This is the default.
    • IQUADTYP=2 sets the Gauss-Legendre quadrature with the transformation ω = ω0 tan(t).
    • IQUADTYP=3 sets the Gauss-Laguerre quadrature scheme.
  5. NQUAD : The variable NQUAD sets the number of quadrature points to be used for the integration. The default is 8.
  6. OMEGA0 and ALPHA : The transformation used in the Gauss-Legendre quadrature schemes involves a constant ω0. The variable OMEGA0 in the namelist SAPTDFT allows you to set this constant. It is typically between 0.3 and 0.5 and the default is 0.5. The Gauss-Laguerre quadrature scheme involves the constant α. For the integrals encountered here, one should set ALPHA=0.0.

The dalton patch enables some new options in dalton which are required for SAPT(DFT) calculations. All the new keywords belong to the *DFT INPUT submodule. The following new keywords are implemented:

  1. .DFTAC — the Fermi-Amaldi asymptotic correction with the Tozer-Handy splicing scheme [63]. The subsequent line should contain 4 numbers: DFTIPT, DFTBR1, DFTBR2, DELTAIPT. The first parameter, DFTIPT, is the ionization potential of the monomer in atomic units. It is recommended to use an experimental value. The parameters DFTBR1 and DFTBR2 are related to the Tozer-Handy switching function and are distances, in Bragg radii, where the switching takes place. Tozer et al. recommend values of 3 and 4 [64]. The last parameter, DELTAIPT, is reserved for the open-shell program and should be set to zero for all closed-shell systems. For open-shell calculation it is equal to IPβ-IPα, where IPβ are IPα are ionization potentials for β and α electrons, respectively.
  2. .CKS — Calculates TD-DFT integrals for the CKS program. Mandatory if the CKS code is used. By default, the LDA kernel is used in TD-DFT (see point 3 below).
  3. .GGAKER — Use generalized gradient approximation (GGA) instead of LDA in the TD-DFT kernel. It is significantly slower and is not recommended for large systems. More details and a discussion of accuracy of the LDA kernel is given in Ref. 14.
  4. .GRAC — Gradient-regulated asymptotic correction (GRAC) of Ref. 36. The subsequent line should contain 4 numbers. The first and last one are identical to the Fermi-Amaldi asymptotic correction (see above) and describe the ionization potential. The second and third numbers describe switching function parameters. Recomended values are 0.5 and 40, respectively.

For the DFT calculations, we recommend the standard Dalton grid or one of the denser grids. For the .THRESH parameter in *HF INPUT, we recommend values of about 10-7.

Since SAPT(DFT) currently works with dalton interface only (and with orca in the density-fitted version), a few scripts are located in misc/daltutil to help in creating dalton input files. Script atmol2dalton converts atmol name.intinp input files into name.mol required by dalton, altergeo.awk updates the name.mol geometry, and with createinputs one can construct name.mol files using basis sets included in the dalton distribution. The usage of the scripts is explained in the sources. The name.dal files can be taken from the examples.

16.3 Terms beyond second order in the interaction operator

The SAPT(DFT) formalism presented above has been implemented up to second order in the operator V . For systems with a significant induction component (e.g., the water dimer), some higher-order terms constitute a large fraction of the total interaction energy. Those terms can be estimated from the supermolecular HF-SCF approach, with δEint,respHF defined in Eq. (5 ). Since no such term can be computed from a pure DFT calculation, this term has to be extracted from a separate SAPT/HF-SCF run. For the calculation of δEint,respHF, no expensive SAPT terms are needed and the scaling is O(N5). Nevertheless, the time of this calculation would be similar to that of the SAPT(DFT) step since the costly transformation step must be done again. To calculate only δEint,respHF, one should prepare inputs as for the regular SAPT runs (any interface can be used as long as the geometry and basis sets are the same as in the SAPT(DFT) run), set DELTASCF=T in the INPUTCOR namelist (cf. Sec. 10.2.3 ), and use scfcp option for the SAPT script. The δEint,respHF value should then be added to the SAPT(DFT) interaction energy. An example of such a calculation is in the directory examples/saptdft/H2O-aTZ-MC.

16.4 Density-fitting version of SAPT(DFT)

Density-fitting (DF) method (also known as resolution of identity) has been implemented in SAPT(DFT) [1557]. The resulting approach has the scaling reduced from O(N6) to O(N5). Substantial savings are also achieved for O(N5) terms. The code with density fitting is also more memory and disk efficient.

Running DF-SAPT(DFT) requires some modifications of the input files and the creation of a file with an auxiliary basis set. The latter file should be named name.aux and contains coefficients of the auxiliary basis function for each atom. The format is the following (see also examples/dfsaptdft):

0 z  
charge 0.0 0.0 0.0  
coefficients  
...  
0 z  
charge 0.0 0.0 0.0  
coefficients

where charge is the charge of an atom and the coefficients are in the turbomole format. We recommend basis sets from Ref. 61. These basis sets are available in turbomole format at ftp://ftp.chemie.uni-karlsruhe.de/pub/cbasen/. The number and ordering of the atoms must be identical to the main basis set used. The quality of the auxiliary basis set is critical for the accuracy of the DF approach and using auxiliary basis set optimized for a different main basis would result in poor accuracy.

In the TRN namelist of the nameP.data, one should specify T2EL=F to suppress the standard two-electron transformation. There is also an additional DF namelist. The only parametr used in this memlist is MEMTRAN=x keyword, where x is equal to either the number of requested memory words or zero. In the latter case, the memory is allocated automatically to the lowest reasonable level. Finally, for running DF-SAPT(DFT), one uses SAPTdf and RunlotDALTONdf scripts instead of SAPT and RunlotDALTON, respectively. Example inputs and scripts for running DF-SAPT(DFT) (analogous to non-DF SAPT(DFT) counterparts in examples/saptdft) are located in examples/dfsaptdft, along with the pertinent outputs. Notice small DF errors when comparing the results of DF and non-DF calculations.

For larger auxiliary basis sets, where the linear dependencies start to emerge, the numerical errors start to amplify and the accuracy of the electrostatic energy is diminished. The problem has been described in Ref. 57. An effective solution is to perform the inversion of the auxiliary J matrix with quadruple precision (QP). For compilers supporting QP, one should use this option to improve the accuracy. The QP code has been tested with ifort and ibm64 platforms. To compile the QP code, one should include QUADINV=YES in the Compall script. If the compiler does not support quadruple precision, the compilation will fail. Since QP calculations are more time consuming than double precision ones, the QP inversion is not turned on by default. To use it, one needs to also specify quad option for SAPTdf or RunlotDALTONdf scripts. With this option, after the SAPT(DFT) calculation is complete, an extra calculation of the electrostatic energy is performed. This result replaces the standard double precision results in the summary table. By comparing the double and quadruple precision values, one can estimate the severity of the quasi-linear dependence. Additionally, slightly better accuracy of density-fitting of the electrostatic energy is obtained by using basis sets of Ref. 65 rather than those of Ref. [61]. One can include such basis set in name-elst.aux and it will be used instead of name.aux in the quadruple precision calculation of the electrostatic energy. If the file is missing, the standard name.aux will be used.

DF-SAPT(DFT) is also partially parallelized. All calculations performed by the (patched) dalton program can be run in parallel. Also the calculations of the CKS kernel integrals are parallelized. The rest of the programs is not parallelized but can be used with parallel dalton.

DF-SAPT(DFT) is limited to up to g functions in the auxiliary basis set. This restriction is due to the use of gamess integrals. In addition, when using the orca interface, the main basis set is currently limited to up to f functions (in practice, both these restrictions are equivalent, because the auxiliary basis set should contain functions with the highest angular momentum L at least by one larger than the main basis set).

The density-fitting code contains some atomic integrals code from gamess. The gamess authors requested that the users of density-fitting SAPT(DFT) should cite gamess [22] along with SAPT.

16.4.1 Using dalton for monomer DFT calculations

The dalton input files required in a DF-SAPT(DFT) run are almost identical with those in regular SAPT(DFT). However, in name.dal files, CKS keyword should be replaced by CKSAUX. This keyword works only with the LDA kernel, therefore, GGAKER keyword cannot be used. It it recommended to include .NOTWO keyword in **INTEGRALS of nameA.dal and nameB.dal when monomer basis sets are used.

16.4.2 Using orca for monomer DFT calculations

Orca is a quantum chemistry package developed by Neese and coworkers [53]. Its DFT module has very efficient density-fitting techniques and we recommend this front-end for DF-SAPT(DFT) calculations, especially in the case of monomers larger than about a dozen of atoms. The Linux, Windows, and Mac-OS X binaries are available at http://www.thch.uni-bonn.de/tc/orca, free of charge for academic users (a registration and acquiring a license is necessary). The orca source files are not needed by DF-SAPT(DFT). Version 2.9.0 or later is required, because earlier versions contain errors in the GRAC asymptotic correction code and lead to incorrect SAPT results. The installation of orca can be done before or after the installation of SAPT, as both are completely independent, but in the latter case one has to edit the file vars.cfg so that the variable ORCADIR is set to the directory with the orca executables.

To use DF-SAPT(DFT) with orca, ISITORCA=T must be set in the TRN namelist of the nameP.data file. The required monomer and/or dimer orca input files follow the same naming pattern (ending with .data) as in the case of gaussian (see Sec. 10.1 for the detailed description). Each such file, besides containing regular orca input in the first section, must end with a special second section starting with a line containing the string META as its first four characters. This section always defines the molecular geometry and basis set used in the DFT calculation (even if a basis set is already defined in the first section, it gets overridden here). Additionally, if one of the density-fitting DFT algorithms implemented in orca is invoked, the META section can be used to override the auxiliary (density fitting) basis set defined before or provided by orca as a default. The example shown below can be used to performed the DFT calculation of the first monomer within the water dimer, using the full dimer-centered basis set (file nameA.data).

! PBE def2-svp def2-svp/J RI TightSCF bohrs  
 
%method  
Grid 4 FinalGrid 5  
gracLB true  
ip 12.62  
end  
 
META  
  6     0 1  
8       0.1088092367800     0.0000000000000    -0.0607755435800   8.0  
1      -0.1616146814400     0.0000000000000     1.7553074350400   1.0  
1      -1.5793331071000     0.0000000000000    -0.7828987377400   1.0  
8      -0.0964834851100     0.0000000000000     6.5351482873200   0.0  
1       0.7718678809000     1.4536519623000     7.2451823417100   0.0  
1       0.7718678809000    -1.4536519623000     7.2451823417100   0.0  
 
BAS 1  
3  
 0 3  
         13.0107010000      0.0334854882  
          1.9622572000      0.2347218706  
          0.4445379600      0.8137702843  
 0 1  
          0.1219496200      1.0000000000  
 1 1  
          0.8000000000      1.0000000000  
 
BAS 8  
6  
 0 5  
        2266.1767785000     -0.0053893504  
         340.8701019100     -0.0402347214  
          77.3631351670     -0.1800818421  
          21.4796449400     -0.4682885766  
           6.6589433124     -0.4469261716  
 0 1  
           0.8097597567      1.0000000000  
 0 1  
           0.2553077223      1.0000000000  
 1 3  
          17.7215043170      0.0626302488  
           3.8635505440      0.3333113849  
           1.0480920883      0.7414863830  
 1 1  
           0.2764154441      1.0000000000  
 2 1  
           1.2000000000      1.0000000000

The first line declares the use of the PBE functional, def2-SVP as the main basis set and def2-SVP/J as the Coulomb-fitting auxiliary basis set, as well as tight convergence criteria for the SCF procedure and the use of bohr units in the geometry section (see the orca manual for details). The gracLB true setting is essential for any SAPT(DFT) calculation as it turns on the GRAC asymptotic correction of the exchange-correlation functional (the other type of asymptotic correction, Fermi-Amaldi, is curently not implemented in orca). The correction requires the ionization potential of the molecule, which must be given in electron-volt units (as opposed to the hartree units used in dalton). The META section starts with a line containing three integers: number of atoms (including ghosts), total charge and multiplicity. For each atom, there is a line containing: atomic number, x, y, and z Cartesian coordinates, and atomic charge (zero for ghost atoms). The main basis set for each element is defined in a block (separated by exactly one empty line from previous blocks) starting with a line containing the string BAS as the first three characters and the atomic number. The next line of the BAS block gives the total number of contracted orbitals. For each such orbital, there is a sub-block starting with two integers (angular momentum and number of primitive functions) and then containing, in separates lines, pairs of exponential Gaussian parameters and contraction coefficients. The order of elements (BAS blocks) is arbitrary but all elements contained in the geometry specification must have an explicitly defined basis set (in particular, it can be specified as empty by setting the number of contracted orbitals to zero), even if this definition is redundant with respect to the basis set keyword in the general orca input. This requirement stems from the fact that the META section information is used not only by orca but also by a separate one-electron integral module, based on routines extracted from the gamess(us) code. Currently, only Gaussian basis functions with angular momenta up to f can be used.

The def2-SVP/J auxiliary basis set requested in the orca input could be replaced by another basis set, by using in the META section additional blocks, with the keyword BAS replaced by AUX. The order of these blocks is also arbitrary but no empty AUX blocks are allowed (this is an orca restriction). Note also that the RI-J Coulomb fitting (requested by the RI keyword in our example) can only be used with pure density functionals. For hybrid functionals, orca has two other methods known as RIJCOSX and RI-JK, the latter requiring “/JK” (rather than “/J”) auxiliary basis sets. All these DFT fitting basis sets must not be confused with auxiliary basis sets required by DF-SAPT(DFT) and defined in the file name.aux, which are designed to fit the correlation effects and must be provided by the user, as in the case of dalton (such basis sets are actually used in orca at the MP2 level and denoted by “/C”).

The orca route in DF-SAPT(DFT) requires a preparation of one input file not occuring in the dalton route, called kernel.data. This file contains just two lines of the form

n_rad n_aug mem  
ksi

where n_rad and n_aug are the numbers of radial and angular grid points around each atom in the quadrature of the CKS kernel integral. In general, much looser grids can be used here compared to the main DFT grids without any significant loss of accuracy. Typically, the values of n_rad = 50 and n_ang = 50 lead to completely negligible errors, while values of n_rad = 30 and n_ang = 26 introduce errors in the dispersion energy of just a few tenths of one percent. Parameter mem specifies the maximum amount of memory (in 8-byte words) to be used in the kernel generation. For larger systems, limiting memory use requires slower multipass calculations. The parameter ksi is the fraction of exact exchange of the employed functional (for instance, 0.25 for PBE0). As the last difference compared to the dalton route, the script dfSAPT.orca rather than dfSAPT is used.