New scaling relations to compute atom-in-material polarizabilities and dispersion coefficients: part 1. Theory and accuracy

Polarizabilities and London dispersion forces are important to many chemical processes. Force fields for classical atomistic simulations can be constructed using atom-in-material polarizabilities and Cn (n = 6, 8, 9, 10…) dispersion coefficients. This article addresses the key question of how to efficiently assign these parameters to constituent atoms in a material so that properties of the whole material are better reproduced. We develop a new set of scaling laws and computational algorithms (called MCLF) to do this in an accurate and computationally efficient manner across diverse material types. We introduce a conduction limit upper bound and m-scaling to describe the different behaviors of surface and buried atoms. We validate MCLF by comparing results to high-level benchmarks for isolated neutral and charged atoms, diverse diatomic molecules, various polyatomic molecules (e.g., polyacenes, fullerenes, and small organic and inorganic molecules), and dense solids (including metallic, covalent, and ionic). We also present results for the HIV reverse transcriptase enzyme complexed with an inhibitor molecule. MCLF provides the non-directionally screened polarizabilities required to construct force fields, the directionally-screened static polarizability tensor components and eigenvalues, and environmentally screened C6 coefficients. Overall, MCLF has improved accuracy compared to the TS-SCS method. For TS-SCS, we compared charge partitioning methods and show DDEC6 partitioning yields more accurate results than Hirshfeld partitioning. MCLF also gives approximations for C8, C9, and C10 dispersion coefficients and quantum Drude oscillator parameters. This method should find widespread applications to parameterize classical force fields and density functional theory (DFT) + dispersion methods.


Introduction
When combined with large-scale density functional theory (DFT) calculations, the DDEC method has been shown to be suitable for assigning atom-centered point charges for exible molecular mechanics force-eld design. [1][2][3][4][5][6][7][8] The assignment of C 6 coefficients and atomic polarizabilities is another active area of research in force eld design. 3,[9][10][11][12] Polarization effects are especially important for simulating materials containing ions. [13][14][15][16][17][18][19][20] When considered alongside the importance of accurate theoretical methods to study van der Waals interactions at the nanoscale, 21 it is clear that a crucial feature of new methods to compute these important quantities is the ability to scale to large system sizes in reasonable computational time.
In this article, we develop new scaling laws and an associated method to compute polarizabilities and dispersion coefficients for atoms-in-materials (AIMs). These new scaling laws and computational method give good results for isolated atoms, diatomic molecules, polyatomic molecules, nanostructured materials, solids, and other materials. This new method is abbreviated MCLF according to the authors' last name initials (where the C is both for Chen and Cole). We performed tests on isolated neutral and charged atoms, small molecules, fullerenes, polyatomic molecules, solids, and a large biomolecule with MCLF. The results were compared with experimental data, high level CCSD calculations, time-dependent DFT (TD-DFT) calculations, or published force-eld parameters.
As discussed in several recent reviews and perspectives, the dispersion interaction is a long-range, non-local interaction caused by uctuating multipoles between atoms in materials. [22][23][24][25][26] It is especially important in (i) condensed phases including liquids, supercritical uids, solids, and colloids, (ii) nanostructure binding such as the graphene layers forming graphite, and (iii) the formation of noble gas dimers. The dispersion interaction is closely related to AIM polarizabilities. The dispersion energy can be described by an expansion series. The leading term is inversely proportional to R 6 , where R is the distance between two atoms. The coefficient of this term is called the C 6 dispersion coefficient, and this term quanties uctuating dipole-dipole interactions between two atoms. The intermolecular C 6 coefficient is given by the sum of all interatomic contributions where M 1 and M 2 refer to the rst and the second molecules. Higher-order terms represent different interactions. 27 For example, the eighth-order (C 8 ) term describes the uctuating dipole-quadrupole interaction between two atoms. The ninthorder (C 9 ) term describes the uctuating dipole-dipole-dipole interaction between three atoms. The tenth-order (C 10 ) term describes the uctuating quadrupole-quadrupole and dipoleoctupole interaction between two atoms. Methods for computing polarizabilities and dispersion coefficients can be divided into two broad classes: (A) quantum chemistry methods that explicitly compute the system response to an electric eld (e.g., TD-DFT, CCSD perturbation response theory, etc.) and (B) AIM models. Class A methods can be highly accurate for computing polarizabilities and dispersion coefficients of a whole molecule, but they do not provide AIM properties. Therefore, class A methods cannot be regarded as more accurate versions of class B methods. Parameterizing a molecular mechanics force-eld from quantum mechanics requires an AIM (i.e., class B) model. Our goal here is not merely to develop a computationally cheaper method than TD-DFT or CCSD perturbation response theory to compute accurate system polarizabilities and dispersion coefficients, but rather to exceed the capabilities of both of those methods by providing accurate AIM properties for force-eld parameterization.
One must distinguish the polarizability due to electron cloud distortion from the polarization due to molecular orientation (and other geometry changes) such as occurs in ferroelectric materials. Herein, we are solely interested in the polarizability that occurs due to electron cloud distortions. Electric polarization due to molecular orientation and other geometric changes can be described by a geometric quantum phase. [28][29][30] There are several existing frameworks for calculating AIM polarizabilities and/or dispersion coefficients. Applequist et al. 31 introduced a formalism that represents the molecular polarizability tensor in terms of AIM polarizabilities via the dipole interaction tensor. Thole 32 rened this formalism by replacing atomic point dipoles with shape functions to avoid innite interaction energies between adjacent atoms. Applequist's and Thole's methods use empirical atomic polarizability tting to reproduce observed polarizability tensors of small molecules. 31, 32 Mayer et al. added charge-charge and dipole-charge interaction terms to calculate more accurate polarizabilities of conducting materials. 33, 34 Grimme et al. presented the D3 geometry-based method to calculate C 6 , C 8 , and C 9 dispersion coefficients and dispersion energies, but this method does not yield polarizability tensors. 11 The recently formulated D4 models are geometry-based methods that extend the D3 formalism by including atomic-charge dependence. [35][36][37] The exchange-dipole model (XDM) is an orbital-dependent approach that yields AIM dipole, quadrupole, and octupole polarizabilities and C 6 , C 8 , and C 10 dispersion coefficients. [38][39][40][41][42] Several density-dependent XDM variants have been formulated. 43,44 In 2009, Tkatchenko and Scheffler introduced the TS method for isotropic AIM polarizabilities and C 6 coefficients. 12 Both the XDM and TS methods yielded isotropic AIM polarizabilities rather than molecular polarizability tensors. 12,[38][39][40][41][42][43][44] In both the XDM and TS methods, the AIM unscreened polarizability is given by (2) where "ref" refers to the reference value for an isolated neutral atom of the same chemical element, "AIM" refers to the partitioned atom-in-material value, and hr 3 i refers to the r-cubed radial moment. Both the XDM and TS methods originally used Hirshfeld 45 partitioning to compute the hr 3 i AIM values. 12,42 In 2012, Tkatchenko et al. introduced self-consistent screening (TS-SCS) via the dipole interaction tensor to yield the molecular polarizability tensor and screened C 6 coefficients. 46 The TS-SCS dipole interaction tensor uses a quantum harmonic oscillator (QHO) model similar to that used by Mayer but extended over imaginary frequencies and omitting chargedipole and charge-charge terms. 33,34,46 That same article also used a multibody dispersion (MBD 46,47 ) energy model based on a coupled uctuating dipole model (CFDM 47,48 ). The TS-SCS screened static polarizability and characteristic frequency for each atom are fed into the CFDM model to obtain the MBD energy. 46 The TS-SCS approach has advantages of yielding a molecular polarizability tensor and AIM screened polarizabilities and AIM screened C 6 coefficients using only the system's electron density distribution as input. 46 The TS-SCS method has several key limitations. Two key assumptions of the TS-SCS method are: (i) for a specic chemical element, the unscreened atomic polarizability is proportional to the atom's hr 3 i moment, and (ii) for a specic chemical element, the unscreened C 6 coefficient is proportional to the atom's polarizability squared. 12,46 However, in their work these hypotheses were not directly tested. 12,46 Later, Gould tested these two assumptions and found them inaccurate for describing isolated neutral atoms placed in a connement potential. 49 Hirshfeld partitioning was used in the TS-SCS method 12,46 to compute the hr 3 i AIM . Because Hirshfeld partitioning uses isolated neutral atoms as references, 45 the Hirshfeld method typically severely underestimates net atomic charge magnitudes. [50][51][52][53] The TS-SCS method assumes a constant unscreened polarizability-to-hr 3 i ratio and constant characteristic frequency (wp) for all charge states of a chemical element, 46 but these assumptions are not realistic. Due to these assumptions, the TS and TS-SCS methods are inaccurate for systems with charged atoms. [54][55][56][57] Bucko et al. showed the TS and TS-SCS methods severely overestimate polarizabilities for dense solids. 57 The TS-SCS method has also not been optimized to work with conducting materials, and we show in Section 5.2 below the TS and TS-SCS methods sometimes predict erroneous polarizabilities even greater than for a perfect conductor. As discussed in Sections 4 and 5 below, the TS-SCS method overestimates directional alignment of uctuating dipoles at large interatomic distances. We also show the TS-SCS method sometimes gives asymmetric AIM polarizability tensors and unphysical negative AIM polarizabilities.
Several research groups developed improvements to the TS-SCS method. Ambrosetti et al. introduced range separation to avoid double counting the long-range interactions in TS-SCS followed by MBD (aka MBD@rsSCS). 58 MBD@rsSCS improves the accuracy of describing directional alignments of uctuating dipoles at large interatomic distances. 58 Bucko et al. used Iterative Hirshfeld (IH 51 ) partitioning in place of Hirshfeld partitioning to compute hr 3 i AIM . 55,57 While this was an improvement, it did not address several problems mentioned above. For example, the TS-SCS/IH method still unrealistically assumed the unscreened polarizability-to-hr 3 i ratio is the same for various charge states of a chemical element. 55,57 This assumption was removed in the subsequent Fractionally Ionic (FI) method. 54 However, the FI approach requires separate reference polarizabilities and C 6 dispersion coefficients for all charge states of a chemical element. 54 This is extremely problematic, because some anions that exist in condensed materials (e.g., O À2 ) have unbound electrons in isolation. 59 Although methods to compute charge-compensated reference ion densities have been developed, 53,55,59-64 those methods do not presently extend to computing charge-compensated polarizabilities and C 6 coefficients of ions. Thus, several problems with the TS-SCS approach have not been satisfactorily resolved in the prior literature.
In Gould et al.'s FI method, reference free atom polarizabilities were computed for various whole numbers of electrons and interpolated to nd fractionally charged free atom reference polarizabilities. 54 This yields different polarizability-to-hr 3 i ratios for different charge states of the same chemical element. 54 Due to the instability of highly charged anions, the À1 states of halogens were the only anions Gould et al. computed self-consistently. 54,65 For other anions, Gould and Bucko resorted to using DFT orbitals from the neutral atoms to build non-self-consistent anions for polarizability and C 6 calculations, 65 but this severely underestimates the diffuseness of anions (i.e., underestimates their polarizabilities and C 6 coefficients). Unfortunately, this problem cannot be easily resolved by performing self-consistent calculations for all charged states, because some highly charged anions (e.g., O 2À ) 59 contain unbound electrons. This makes the FI method problematic for materials containing highly charged anions. Because FI was not included in VASP version 6b that we currently have access to, we did not perform FI calculations for comparison in this work.
In this article, we develop a new approach that resolves these issues. Our method uses DDEC6 (ref. 61 and 66-68) partitioning to provide accurate net atomic charges (NACs), atomic volumes, and radial moments as inputs. Our method has new scaling laws for the unscreened atomic polarizabilities, characteristic frequency (wp), and C 6 dispersion coefficients. The different scaling behaviors of surface and buried atoms are included via m-scaling. Our approach accurately handles the variability in polarizability-to-hr 3 i moment ratio for charged surface atoms while only requiring reference atomic polarizabilities, reference C 6 coefficients, and reference radial moments for isolated neutral atoms. It uses a new self-consistent screening procedure to compute screened polarizability tensors and C 6 coefficients. Our approach separates non-directional screening from directional screening of the dipole interaction tensor. This allows a conduction limit upper bound to be applied between nondirectional and directional screening to ensure buried atoms do not have a screened polarizability above the conduction limit. Our method yields three different types of dipole polarizabilities: (a) induced static polarizabilities corresponding to a uniform applied external electric eld, (b) isotropic screened polarizabilities suitable as input into polarizable force-elds, and (c) uctuating polarizabilities that are used to compute C 6 dispersion coefficients via the Casimir-Polder integral. When computing C 6 dispersion coefficients, we use multi-body screening to taper off the dipole directional alignment at large interatomic distances. The self-consistent screening is applied incrementally. Richardson extrapolation provides high numeric precision. Through quantum Drude oscillator (QDO) parameterization, our method also yields higher-order polarizabilities (e.g., quadrupolar, octupolar, etc.) and higher-order dispersion coefficients (e.g., C 8 , C 9 , C 10 , etc.) for AIMs. Other important improvements include: improved damping radii, proportional partitioning of shared polarizability components to avoid negative AIM polarizabilities, iterative update of the spherical Gaussian dipole width, and AIM polarizabilities are symmetric tensors.
Our method was designed to satisfy the following criteria: (1) The method should require only the system's electron and spin density distributions as input; (2) The method should work for materials with 0, 1, 2, or 3 periodic boundary conditions; (3) The method should give accurate results for charged atoms in materials while only requiring reference polarizabilities and reference C 6 coefficients for neutral free atoms; (In this context, reference polarizabilities and reference C 6 coefficients refer to those polarizabilities and C 6 coefficients stored within the soware program that are used during application of the method.) (4) The method should give accurate results for diverse materials types: isolated atoms; small and large molecules; nanostructured materials; ionic, covalent, and metallic solids, etc.; (5) The method should give accurate results for both surface and buried atoms; (6) The method should yield both static polarizability tensors and polarizabilities suitable for constructing molecular mechanics force-elds; (7) The method should accurately describe both short-and long-range ordering of dipole polarizabilities and C 6 coefficients; (8) The method should have less than approximately 10% average unsigned error on C 6 coefficients and dipole polarizabilities for the benchmark sets studied here; (9) The method should include estimates for higher-order AIM polarizabilities (e.g., quadrupolar, octupolar, etc.) and dispersion coefficients (e.g., C 8 , C 9 , C 10 , etc.); (10) The method should have low computational cost for both small and large systems.
There are two main applications for this MCLF method. First, the polarizabilities and dispersion coefficients can be used to partially parameterize molecular mechanics force-elds. In addition to polarizabilities and C 6 dispersion coefficients, those force-elds would also need to include net atomic charges (NACs), exibility parameters (e.g., bond, angle, and torsion terms), exchange-repulsion parameters, (optionally) charge penetration parameters, and optionally other parameters. Second, the dispersion coefficients can be used to partially parameterize DFT + dispersion methods. 23 In addition to the C 6 dispersion coefficients, an accurate DFT + dispersion method should also include higherorder dispersion (e.g., C 8 , C 9 , and/or C 10 terms) or multi-body dispersion (MBD) combined with an accurate damping function. 22,23 (Partially analogous to the MBD@rsSCS method, 58 range separation would be required to avoid double counting dispersion interactions when combining a MCLF variant with a MBD Hamiltonian.) Because DFT and molecular mechanics are widely used in computational chemistry, our new method can have widespread applications.
The remainder of this article is organized as follows. Section 2 contains the background information. Section 3 contains the new isolated atom scaling laws developed for MCLF. Section 4 describes the theory of the MCLF method. Section 5 contains calculation results of C 6 coefficients and polarizabilities and comparisons to benchmark data. Section 6 is the conclusions.

Benchmarking methods
Experimental data and high-level quantum chemistry calculations were used as references in this work. Section 3.1 below describes reference polarizabilities and dispersion coefficients (C 6 , C 8 , and C 10 ) for the isolated atoms. For diatomic molecules in Section 5.1, we computed static polarizability tensors using CCSD calculations combined with the "polar" keyword in Gaussian09 (ref. 69). As explained in Section 5.2 below, we set the reference static polarizability for dense solids to the lesser of the Clausius-Mossotti relation value and the conduction limit upper bound based on the experimental crystal structure geometry and dielectric constant.
For the small molecules in Section 5.3, reference static polarizabilities were obtained from published experiments. Experimental isotropic polarizabilities were extracted from dielectric constant or refractive index measurements having approximately 0.5% or less error. 70 Refractive index can be measured by passing a light ray through a gas-phase sample. 71 The polarizability a(n) of the sample at frequency n can then be calculated using where h is the refractive index, m is the dipole moment magnitude, T is absolute temperature, L is the volume per molecule, and k B is Boltzmann's constant. 70 (The last term in eqn (3) accounts for orientational polarizability.) Static or low-frequency dielectric constants k were obtained by measuring the ratio of the capacitance of a set of electrodes with the sample material inbetween to the capacitance of the same electrodes with vacuum in-between. 72 The polarizability of a gas-phase sample can then be calculated using the Clausius-Mossotti relation: 73 Reference C 6 coefficients for the atom/molecule pairs (Section 5.4) were taken from the experimentally extracted dipole oscillator strength distribution (DOSD) data of Meath and co-workers 74,75 as tabulated by Bucko et al. 55 Time-dependent DFT (TD-DFT) and time-dependent Hartree-Fock (TD-HF) can be used to compute benchmark polarizabilities and dispersion coefficients. The Casimir-Polder integral is used to calculate C 6 coefficients from polarizabilities at imaginary frequencies (imfreqs). 76 For polyacenes (Section 5.5), reference C 6 coefficients and isotropic static polarizabilities are from the TD-DFT calculations of Marques et al. 77 For selected polyacenes, static polarizability tensor components were available as reference from Jiemchooroj et al.'s TD-DFT calculations. 78 Jiemchooroj et al. found their TD-DFT results were similar to TD-HF, experimental (where available), and CCSD (where available) results. For fullerenes (Section 5.5), the reference C 6 coefficients and isotropic static polarizabilities are from the TD-HF calculations of Kauczor et al. 116 Kauczor et al. also obtained similar results using TD-DFT. 116

Notation
A system may have either 0, 1, 2, or 3 periodic boundary conditions. In periodic materials, the term 'image' refers to a translated image of the reference unit cell. Each image is designated by translation integers (L 1 , L 2 , L 3 ) that quantify the unit cell translation along the lattice vectors. The reference unit cell is the image designated by (L 1 , L 2 , L 3 ) ¼ (0, 0, 0). ÀN # L i # N along a periodic direction. L i ¼ 0 along a non-periodic direction. Similar to the notation previously used in the bond order article, 67 a capital letter (A, B,.) designates an atom in the reference unit cell and a lowercase letter (a, b,.) designates an image atom. For example, b ¼ (B, L 1 , L 2 , L 3 ) denotes a translated image of atom B.
LetR B represent the nuclear position of atom B in the reference unit cell. Then, the nuclear position of a translated image is whereṽ (1) ,ṽ (2) , andṽ (3) are the lattice vectors. The distance between the nuclear position of atom A and the translated image of atom B is Cartesian components (s ¼ x, y, z) of the vector from atom A's nuclear position to image b's nuclear position are represented by r A is the vector from the image of atom A's nuclear position to the spatial positionr: The length of this vector is represented by A stockholder partitioning method assigns a set of atomic electron densities {r A (r A )} in proportion to atomic weighting factors {w A (r A )} so that all sum to the total electron density where summation over A, L means summation over all atoms in the material. The number of electrons N A and net atomic charge (q A ) assigned to atom A are where Q A is the atomic number for atom A. As discussed in Section 2.4 below, different ways of dening {w A (r A )} produce different stockholder methods. The AIM radial moment of order f is hr 2 i, hr 3 i, and hr 4 i are shorthand for h(r A ) 2 i, h(r A ) 3 i, and h(r A ) 4 i, respectively. Since a particular dispersion-polarization model can be combined with different charge partitioning methods, we indicate the combination by the dispersion-polarization model followed by '/' followed by the charge partitioning method. For example, TS-SCS/IH indicates the TS-SCS dispersionpolarization model using IH charge partitioning. Where our computed data are simply labeled 'TS-SCS', the DDEC6 charge partitioning method was used. Since all of the MCLF results reported in this paper used DDEC6 charge partitioning, we used the less precise but shorter term 'MCLF' in place of the full 'MCLF/DDEC6' designation. More generally, the MCLF dispersion-polarization model could potentially be combined with other charge partitioning methods (e.g., MCLF/IH), but that is beyond the scope of the present study.
Calculating dispersion coefficients involves integrating polarizabilities over imfreqs. This is inconvenient in two respects. First, it is easier to deal with real-valued variables rather than imaginary-valued ones. Second, numeric integration from zero to innite imaginary frequency is inconvenient, because innity cannot be readily divided into nite intervals for numeric integration. Letting u represent an imaginary frequency magnitude, we used the following variable transformation to solve these two problems: This conveniently transforms integration limits u ¼ [0, N) into u ¼ [Nimfreqs, 0), which upon changing the integrand's sign gives integration limits u ¼ (0, Nimfreqs]. As shown in the companion article, this allows convenient Rhomberg integration using integration points u ¼ 1, 2,. Nimfreqs. 79 In this article, a(u) denotes the polarizability at the imaginary frequency whose magnitude equals u(u).  (15). For atoms A and B, they dene a many-body polarizability matrix P, with its inverse Q having the form (16) where s, t designate Cartesian components. Square matrices P and Q have x, y, and z spatial indices for every atom to give a total of 3Natoms rows. The last term on the right-hand side is the dipole interaction tensor which has the form

Details of the TS-SCS methodology
where summation over L means summation over all periodic translation images (if any). s AB (u) is the attenuation length for the pair of atoms A and B The spherical Gaussian dipole width is obtained from and the isotropic dynamical atomic polarizability is In the TS and TS-SCS methods, where a TS is calculated by eqn (2). AIM polarizability tensors are computed using the partial contraction with the static polarizability tensor corresponding to u ¼ Nimfreqs. The screened frequency-dependent isotropic polarizability is computed as one third of the trace of the three-by-three polarizability tensor obtained by partial contraction of P These are fed into the Casimir-Polder integral expressed in terms of u (see companion article for derivation 79 ) to compute C 6,A .

Electron density partitioning methods
In Hirshfeld partitioning introduced in 1977, atoms are partitioned to resemble the neutral reference atom. 45 This makes the atoms tend to have lower charge than they should have. 50 The iterative Hirshfeld partition (IH) keeps updating the reference with the charge state of the atom. 51 However, this approach leads to the runaway charge problem in some cases. 61 As shown in ref. 57 and Section 5.2 below, using TS or TS-SCS with Hirshfeld or IH partitioning overestimates polarizabilities for dense solids. Manz and Sholl presented DDEC1 and 2 atomic population analysis methods in 2010. 60 By simultaneously optimizing the AIM density distributions to be close to spherically symmetric and to resemble charge-compensated reference ion densities, this method can give chemically meaningful NACs and accurate electrostatic potential for some materials, but was later found to give runaway charges in other materials. In 2012, Manz and Sholl presented the DDEC3 method that partially xes the runaway charge problem by increasing the optimization landscape curvature via conditioned reference densities and imposing an exponential decay constraint on each atom's electron density tail. 53 Manz and Limas presented DDEC6 partitioning in 2016. 61,66 This method xes the runaway charge problem. Also, new constraints are added to the decay rate of the buried atom tails. The weighted spherical average improves the effect of spherical averaging during charge partitioning. Along with guaranteed convergence in seven steps, this method is very accurate, cost efficient, and produces chemically meaningful NACs. 61,66,68 In 2017, Manz published a new method for computing bond orders, which is based on DDEC6. 67

Reference data
The reference polarizabilities (a CCSD ) used in this work are our calculated polarizabilities using the CCSD method with def2QZVPPDD basis set (the def2QZVPPDD basis set is dened in the ESI †). We tested two different methods: (a) using Gaussian09 (ref. 69) keyword 'polar' to compute the molecular static polarizability tensor using perturbation response theory and (b) using Gaussian09 keyword 'eld' to manually apply a small constant external electric eldẼ in order to calculate the molecular static polarizability tensor as a the molecular dipole moment. However, many of the manual (i.e., keyword ¼ 'eld') calculations did not converge and the converged results were not as consistent with Gould and Bucko's data 65 as the perturbation response calculations. So we decided to use the perturbation response calculations (i.e., keyword ¼ 'polar') for all elements except Y. For Y, the keyword ¼ 'eld' polarizability was used, because the perturbation response calculation gave an unreasonably low polarizability of 88.98 compared to a Gould ¼ 163 while the manually applied eld polarizability of 158.81 was close to Gould and Bucko's value and followed the trend of neighboring elements: a Sr,CCSD ¼ 204.51 a Zr,CCSD ¼ 143.47. Fig. 2 shows that our calculated polarizabilities are mostly consistent with Gould and Bucko's values. We used a CCSD rather than a Gould as the reference free atom polarizability, because our radial moments come from the same CCSD calculations as used to compute a CCSD . For elements using a relativistic effective core potential (RECP) in the def2-QZVPPDD basis set, the radial moments of core electrons replaced by the RECP are added back in using a reference core density library; thus yielding effective all-electron radial moments. Since Gould and Bucko used the auau principle for electron congurations of transition metal atoms, their calculations do not necessarily correspond to the ground state spin multiplicity for transition metal atoms. 65 Recently, Schwerdtfeger and Nagle published a set of recommended polarizabilities for chemical elements 1 to 120 (except livermorium), but those were not used in our study because they do not have a corresponding set of C 6 values. 80 CCSD in Gaussian09 does not have the capability of calculating C 6 coefficients. Therefore, to maximize consistency between the free atom reference radial moments, polarizabilities, and C 6 coef-cients, our reference C 6 coefficients were calculated as where a Gould and C Gould 6 are Gould and Bucko's values using TD-DFT. 65 This C 6 rescaling makes C ref 6 correspond to a CCSD , which corresponds to the computed radial moments. The 3/2 power occurs in eqn (25), because a and C 6 for a free atom are approximately proportional to the free atom's effective radius to the fourth and sixth powers, respectively (see Table S1 of ESI †).
The reference a, C 6 , wp, r moments, and damping radii are listed in the ESI. † This dataset contains neutral elements 1-86 except the f-block elements (58)(59)(60)(61)(62)(63)(64)(65)(66)(67)(68)(69)(70)(71). The reason for excluding the f-block and heavier elements is the def2QZVPPDD basis set we are using does not include these elements. The dataset also includes +1 cations of elements 3-7, 11-17, 19-57 and 72-86 and À1 anions of F, Cl, Br, I, and At. These ions were also selfconsistently calculated by Gould and Bucko. 65 Gould and Bucko included additional anions which were not computed selfconsistently, and we omit these because self-consistent polarizabilities are unavailable. 65 The reference wp were calculated from the CCSD polarizability and the C ref 6 using 81 The reference C 8 and C 10 are from several sources compiled by Porsev and Derevianko 82 and Tao et al. 83 and are listed in the ESI. † This dataset contains H, Li, Na, K, Rb, Cs, He, Ne, Ar, Kr, Xe, Be, Mg and Ca. Most of these reference values are based on correlated quantum chemistry calculations.

Deriving the new scaling laws
Johnson and Becke assumed that for a given chemical element, the polarizability of atoms in a material should be proportional to the hr 3 i moment of the atom-in-material. 42 This assumption was subsequently adopted by Tkatchenko and Scheffler when formulating the TS method. 12 Of course, this is not the same as assuming polarizabilities of isolated atoms across different chemical elements should be proportional to their hr 3 i moments. As pointed out by Gould, the isolated atom polarizabilities are not proportional to their hr 3 i moments. 49 Fig. 3 is a plot of ln(polarizability) versus ln(hr 3 i) for the isolated atoms. This plot shows a weak correlation (R 2 ¼ 0.7706). This motivated us to develop a new polarizability scaling law that applies both to isolated atoms and to atoms-in-materials across different chemical elements.
We tested seven models containing electron number and different combinations of the r moments as the independent variables and a, C 6 , and wp as the dependent functions. Table 1 lists the R 2 values of the 7 models. The coefficients were obtained by least squares tting of a linear combination of the log values of r moments and electron numbers to a, C 6 , or wp using  a Matlab program we wrote. For example, the entry in row 2 and column 2 in Table 1 is the R 2 value of 0.6626 obtained by tting log(hri) and log of electron number to log(a CCSD ). The results show that using only one r moment does not yield high R 2 value. Combinations of two or more r moments give higher R 2 values, with the hr 3 i & hr 4 i model giving the best average performance. Table 2  To test the robustness and transferability of the different models, the following tests were performed as shown in Table 3. The "PW91 retted" column are the R 2 values obtained by retting the model parameters with r moments from PW91. Manz and Limas performed these PW91 calculations in Gaussian09 using an all-electron fourth-order Douglas-Kroll-Hess relativistic Hamiltonian with spin-orbit coupling (DKHSO) and the MUGBS basis set near the complete basis set limit employing a nite-size Gaussian nuclear model. 60,61 The "PW91 predicted" column are the R 2 values calculated using CCSD model parameters from Table 1 but with PW91 r moments instead of CCSD r moments. Table 3 shows that the hr 3 i & hr 4 i model has the highest R 2 in both tests; therefore, this model is the most robust and transferable.
Hence, the new polarizability scaling law for an isolated atom is where N is the number of electrons, the superscript "ref" means the value of the neutral atom reference, and "AIM" means the value for atom-in-material aer partitioning. The new wp scaling law for an isolated atom is   for an isolated surface atom is then computed via eqn (26). These scaling laws allow different charge states of an atom to be accurately described while using only reference polarizability and wp values for a neutral free atom of the same element. For a, the effective power of the effective radius is 4 Â 3.3372 À 3 Â 3.1657 ¼ 3.8517, which is approximately 4. For wp, the effective power of the effective radius is 3 Â 3.7003 À 4 Â 3.2228 ¼ À1.7903, which is approximately À2. Scaling laws for non-isolated atoms will be addressed in Section 4 below.

Higher-order dispersion coefficients and quantum Drude oscillator parameters
In this section, we consider higher-order dispersion coefficients C 8 , C 9 , and C 10 and their mixing rules. The contribution of the three-body C 9 term to the dispersion energy is typically less than 10%. 11 Nevertheless, McDaniel and Schmidt 84 showed that in order to obtain accurate results from force-eld simulations for condensed phases, the three-body term (E ABC ) should be included. Tang and Toennies showed that the attractive potential at well depth for two free atoms is mainly composed of C 6 , C 8 , and C 10 with contributions of roughly 65%, 25%, and 7% respectively. 85 The rest comes from higher-order terms. Because the C 8 , C 9 , and C 10 terms have modest contributions, we decided to include them in our model.
The C 8,A coefficient describes the uctuating-dipole-uctuating-quadrupole two-body dispersion interaction between atoms of the same type. We dened two groups (with all quantities expressed in atomic units) for least-squares tting to obtain a model for C 8,A : The reason for using hr 3 i and hr 4 i is that these are the same r moments used in models discussed above. Since C 8,A describes the uctuating dipole-quadrupole coupling while C 6,A describes the uctuating dipole-dipole coupling, there is no reason to believe directional effects on C 8,A follow those on C 6,A . Therefore, our correlations for higher-order dispersion coefficients (i.e., C 8 , C 9 , and C 10 ) do not include directional coupling. C non-dir 6,A is obtained using the imfreq-dependent non-directionally screened atomic polarizability a non-dir A (u) in the Casimir-Polder integral. Linear tting was performed to obtain the slope and intercept for group 2 as a function of group 1. The results were 0.8305 for the slope and 1.7327 for the intercept yielding: The top le panel of Fig. 7 shows strong correlation between the model predicted C 8,A and the reference data 82,83 for selected isolated atoms with MARE of 14.6%. The Quantum Drude Oscillator (QDO) model provides a natural framework for describing multibody polarizability and dispersion interactions beyond the dipole approximation, including quadrupolar, octupolar, and high-order interactions. 9,86,87 A QDO consists of a negative pseudoparticle coupled via a harmonic potential to a pseudonucleus. 9,86,87 This harmonic coupling produces a Gaussian charge distribution. 9,86 In our model, one QDO is centered on each atom in the material. Each QDO is completely described by three parameters: (a) an effective mass (m QDO ), (b) an effective charge (q QDO ), and (c) an effective frequency (wp QDO ). 9,86 Using literature relations 9 applied to our non-directionally screened quantities yields Fig. 8 Illustration of MCLF atom-in-material polarizability tensors for the carbonyl sulfide molecule. The sulfur atom had much higher polarizability than the carbon and oxygen atoms. All three atoms showed enhanced polarizability along the bond direction. Only the relative sizes and shapes of the ellipsoids were drawn to scale. Table 4 Comparison of the % error in the polarizability of 28 solids as a function of C value. "NU" means no upper bound is applied and "U" means the upper bound is applied a xx /a a yy /a a zz /a a xx /a a yy /a a zz /a a xx /a a yy /a a zz /a The C 9,ABC QDO mixing rule 9 is similar to that described much earlier by Tang using a Padé approximation. 81 The three-body uctuating dipole-dipole-dipole interaction energy term is E 9 ¼ When constructing a force-eld using the MCLF dispersion coefficients, care should be taken not to double-count the threebody dipole-dipole-dipole interactions. Specically, the MCLF directional screening (i.e., C screened 6 ) already includes the three-body dipole-dipole-dipole interactions for the calculated system. (These were included via the directional dipole interaction tensor which was used in turn to compute C screened 6 .) For example, if C screened 6 is computed using the MCLF method for a single benzene molecule, then the intramolecular dipole-dipole-dipole interactions are already included via the C screened 6 coefficient term, but the intermolecular dipole-dipole-dipole interactions are not already included and should be added in the force-eld using the C 9 coef-cient term. In this case, the force-eld's C 9 three-body term should be constructed to include atom triplets from two or three molecules, but not from a single molecule. Fig. 7 compares model predicted to reference C 8 and C 10 dispersion coefficients. Reference values are from Porsev and Derevianko 82 and Tao et al. 83 and sources cited therein. Chemical elements included in Fig. 7 are: (a) the alkali metals Li, Na, K, Rb, Cs, (b) the noble gases He, Ne, Ar, Kr, Xe, (c) hydrogen H, and (d) the alkaline earths Be, Mg, Ca. Plotted C 8,AB and C 10,AB included 73 different heteroatomic pairs (see ESI † for details) formed from a subset of these elements.

Extension to element 109
The CCSD reference data are missing for elements 58-71 and 87-109 because Def2QZVPPDD basis set is not available for those elements. In order for the method to be more applicable, we used eqn (26)- (28) to estimate the reference polarizability, wp, and C 6 for these elements. This is a reasonable approach since Table 3 already showed that CCSD model with PW91 r moments yields reasonable results. For these elements, the PW91 hr 3 i and hr 4 i moments required as inputs were taken from all-electron fourth-order DKHSO calculations with MUGBS basis set near the complete basis set limit employing a nite-size Gaussian nuclear model. 60,61 Furthermore, data for La (element 57, an f block element) is available for both CCSD and PW91. Comparing its PW91 a and wp with CCSD ones, the change is 5.4% and 2.7% respectively. The reference data for elements 58-71 and 87-109 are listed in ESI. † Fig. 9 Flow diagram for MCLF method. (a) the C 6 coefficient from the Casimir-Polder integral will be unphysical, (b) corresponding vdW radii in the TS-SCS method will be complex (since they depend on the cubed root of the polarizability), 12 and (c) MBD, MBD@rsSCS, and related methods require non-negative AIM polarizabilities as input. 46,54,58 In TS-SCS, the partial contraction of P follows the assumed form of Applequist et al.: 31 This sometimes yields negative AIM polarizabilities, because the mixed contribution P AB st (which might be negative) between an atom A with small polarizability and an atom B with large polarizability can surpass the magnitude of P AA st . In our new method, atoms with larger pre-screened polarizability get a proportionally larger piece of the screening mixed polarizability contribution: Eqn (45)- (47) correspond to the non-directional, uctuating, and static polarizabilities, respectively. Aer this new partition is applied, the oxygen in ZrO has the MCLF polarizability of 6.754. Also, eqn (46) ensures the AIM polarizability P A st is a symmetric tensor just like the total polarizability tensors 88,89 of all molecules. In contrast, the TS-SCS polarizability tensor for an atom-inmaterial is sometimes asymmetric with respect to the spatial coordinates (e.g., the xy and yx components are different). As an example, the ESI † contains the TS-SCS and MCLF results les for dibromomethane, where the TS-SCS yz and zy polarizability components for the last Br atom were 4.05 and 7.31, respectively; the MCLF method gave 5.74 for both components.
The symmetric AIM polarizability tensor can be visualized by plotting the ellipsoid 88,89 r where l 1 , l 2 , and l 3 are its three eigenvalues andê I ,ê II , andê III are the corresponding mutually orthogonal eigenvectors. Fig. 8 plots AIM polarizability tensors for the carbonyl sulde molecule. All three atoms showed enhanced polarizability along the bond direction. The atom-in-material polarizability along a unit directionk is quantied by the projectionk$ a ! ! A $k. Choosingk parallel to the inter-nuclear direction gives the bond polarizability of two bonded atoms:k$ð a ! ! A þ a ! ! B Þ$k. 88,89 Of interest, Raman spectrum peak intensities are proportional to the change in projected polarizability as vibration occurs. [90][91][92]

Polarizability upper bound
Consider a perfectly conducting plate with thickness d in an external electric eld E ext applied perpendicular to its surface. From Gauss' Law, the two faces perpendicular to E ext develop surface charge densities s 0 ¼ AE3 0 E ext and form a dipole moment The local polarizability must be computed based on the electric eld (E loc ) felt by each local volume element within the material. In this geometry, E loc = E ext /2 is the average of the eld before (E ext ) and aer (E nal = 0 inside the conductor) polarization. Its local polarizability is p ¼ m/E loc ¼ 23 0 d Â area, which in atomic units (4p3 0 ¼ 1 in atomic units) gives a ¼ p/(4p3 0 ) ¼ d Â area/(2p), also referred to as the polarizability volume since it carries volume units. Thus, the local polarizability-to-volume ratio of the perfectly conducting plate is 1/(2p) in atomic units. For this geometry, the overall polarizability-to-volume ratio (p ¼ m/E ext ) is 1/(4p).
As a second example, consider a sphere of radius R and dielectric constant k placed in a constant externally applied electric eld E ext along the z-direction. This sphere will develop a dipole moment of 93 Therefore, its overall polarizability in atomic units is R 3 (k À 1)/(k + 2), and its overall polarizability-to-volume ratio is (3/(4p))((k À 1)/(k + 2)) in atomic units. The solution for a perfect conductor can be obtained by taking the limit k / N which yields 3/(4p) as the overall polarizability-to-volume ratio for the conducting sphere.
In theory, the polarizability caused by distortion of the electron cloud for xed nuclear positions should be less than or equal to that of a perfect conductor. Comparing the above results for the conducting plate and conducting sphere shows the overall polarizability-to-volume ratio of a perfect conductor is shapedependent; this is due to directional interactions within the material. To address this issue, we apply a conduction limit upper bound during non-directional screening, before any directional screening occurs. Applying the conduction limit upper bound before directional screening allows it to be based on the local polarizablity-to-volume ratio instead of the overall polarizability-tovolume ratio. The local polarizability-to-volume ratio should be approximately independent of the material's shape. Therefore, the local polarizability-to-volume ratio of the conducting plate is operational for the non-directional screening. So we dened the polarizability upper bound of atom A as where V A is the volume of atom A in the material. During MCLF calculation, if the calculated non-directionally screened polarizability is higher than this upper bound, the result will be replaced by the upper bound polarizability. This procedure is carried out by a smooth minimum function where Y ¼ 25 because this value provides a smooth curve with less than 0.58% deviation from max(min (a, b), 0). If either a or b is #0, the function is set to return a value of 0. Using a smooth curve, instead of simply min(a, b), is required to ensure the forces (i.e., energy derivatives with respect to atom displacements) are continuous functions.

M scaling to describe both surface and buried atoms
Section 3.2 above shows the polarizability of an isolated atom scales approximately proportional to its effective radius to the 4th power. Brinck et al. showed the polarizability of a polyatomic molecule is approximately proportional to its volume divided by an effective electron removal energy. 94 As shown in Section 4.2, the polarizability of a conducting plate or sphere is proportional to its volume. For non-conducting uids, the Clausius-Mossotti relationship describes a polarizability-to-volume ratio that is a weak function of the dielectric constant. 73 This implies the polarizability-to-volume ratio of a buried atom is approximately proportional to its volume. This means the atom-in-material polarizability transitions from approximately 4th power to 3rd power dependence on the atom's effective radius as the atom goes from isolated to buried. We dene m to quantify how exposed an atom is: where r A (r A ) is the electron density of atom A, w A (r A )/W(r) is the fraction of the total electron density r(r) at positionr that is assigned to atom A, and h i r A means the spherical average. m equals 1 for an isolated atom and goes towards zero when an atom gets more and more buried. The modied unscreened scaling law for a now has the form where C is a constant to be determined. When m ¼ 1, it becomes the isolated atom scaling law in Section 3. When m ¼ 0, it becomes There are two primary justications for eqn (57). First, the Clausius-Mossotti relation and conduction limit upper bound show the polarizability of a condensed phase (e.g., liquid or solid) is approximately proportional to its volume. Here, hr 3 i is a proxy for the volume of an atom, so the sum of hr 3 i moments for atoms in a material is a proxy for the material's volume. The Clausius-Mossotti relation and conduction limit upper bound show the polarizability-to-volume ratio of a material is not strongly elementdependent and has modest dependence on the material's dielectric constant (see eqn (75) and (76)). Second, in condensed phases electrons undergo chemical potential equilibration that transfers some electron density from the least electronegative elements to the most electronegative elements. Accordingly, the chemical potential of the equilibrated condensed phase is not as extreme as either the most electropositive elements nor the most electronegative elements. For atoms in condensed materials, the polarizability to hr 3 i moment ratios will typically be lower than an isolated neutral alkali metal atom (extremely electropositive) on the one hand and higher than an isolated neutral uorine atom (extremely electronegative) on the other hand. The polarizability to hr 3 i moment ratios of atoms in condensed materials thus exhibit a narrower range of values than for isolated atoms. Group 14 elements have approximately equal tendency to gain or lose electrons, thereby readily forming both positive and negative oxidation states. The isolated atom polarizability to hr 3 i moment ratios of Group 14 elements are approximately independent of the periodic row: C (0.34), Si (0.37), Ge (0.34), Sn (0.32), and Pb (0.31). Accordingly, we expect the same constant C in eqn (57) will work for both light and heavy elements.
Tests on the polarizabilities of 28 solids were performed to optimize C. The criteria are that without any upper bound imposed, the MRE should be 0-10% and aer the upper bound is imposed, the MARE should be #$10%. These criteria make sure the scaling law is as accurate as possible without the upper bound and the upper bound will not decrease the accuracy. Table 4 shows that C ¼ 0.4 is optimal. This value is slightly higher than the polarizability to hr 3 i moment ratios of isolated Group 14 atoms.
The characteristic frequency wp also has different expressions for isolated and buried atoms. For an isolated atom (i.e., m ¼ 1), wp should be proportional to a À1/2 which has been captured by the isolated atom scaling law. For a buried atom (i.e., m ( 1), the polarizability becomes nearly proportional to hr 3 i AIM and C 6 remains proportional to the atom's effective radius to the sixth power; therefore, wp is less sensitive to changes in a when the atom is buried. So the scaling for wp has been modied to: z ¼ U when m ¼ 1, eqn (60) reduces to the isolated atom scaling law in Section 3. A complete explanation and derivation of eqn (60) is given in Section S3 of ESI. †

Iterative polarizability screening
The non-directional, uctuating, and induced dipole interaction tensors are dened in eqn (61) the subscripts s and t refer to spatial indices x, y, and z. This enables the non-directional and directional components to be separately screened. As described in Section 4.2, the conduction limit upper bound applies to the non-directionally screened (not the directionally screened) components. Computing static polarizabilities (i.e., system response to constant externally applied electric eld) requires directionally screened polarizabilities. As discussed in Section 4.7, parameterizing polarizable force elds requires non-directionally screened polarizabilities. f cutoff (d Ab ) is a smooth cutoff function that smoothly turns off dipole-dipole interactions between atoms as the distance between them increases: where d Ab ¼ r AB,L is the distance between atom A and the image b in bohr. d cutoff is the dipole interaction cutoff length. H is the Heaviside step function. As shown in Fig. S1 of ESI, † this function smoothly decreases from $1 at d Ab ¼ 0 to zero when d Ab $ d cutoff . The power of three in the exponential ensures that both the rst and second derivatives are continuous at d Ab ¼ d cutoff , which is required for frequency calculations. The factor of 20 in the exponent provides a balanced cutoff function steepness. Specically, f cutoff (d Ab # 0.5 d cutoff ) $ 0.9179 ensuring that all positions within half the cutoff distance are counted at nearly full strength. Imagine a parameter 0 # l # 1 that continuously turns on a particular screening type (e.g., non-directional, uctuating, or static). When l ¼ 0, the corresponding screening type is fully off. When l ¼ 1, the corresponding screening type is fully on. Thus, the corresponding screening process can be envisioned as transitioning continuously from l ¼ 0 to l ¼ 1. We expect the partially screened system (i.e. 0 < l < 1) to have a polarizability intermediate between that for l ¼ 0 and l ¼ 1. Since the Gaussian width s AB (u) depends on the polarizability (eqn (19)), the MCLF method continuously updates s AB (u) during the screening process. In contrast, the TS-SCS method uses s AB (u) corresponding to l ¼ 0 for the entire screening process, which does not allow the Gaussian width to update during the screening process.
The screening of each tensor is now an iterative process. In each iteration, we compute where D j is the screening increment. Then P j+1 is computed as the inverse of Q j+1 . The screening process is divided into nondirectional and directional screening. Directional screening has two separate parts: screening for the static polarizability and screening for the uctuating polarizabilities. Dividing the screening process into these three parts allows us to compute three different types of polarizabilities having distinct uses: (a) force-eld polarizabilities that are suitable inputs for polarizable force-elds, (b) uctuating polarizabilities for computing C 6 coefficients via the Casimir-Polder integral, and (c) static induced polarizabilities corresponding to a constant externally applied electric eld. As described in Section 4.2, this division also allows a polarizability upper bound to be applied. These three parts of the screening process are summarized below. Non-directional screening. For non-directional screening, both D and s are Natoms by Natoms matrices. In the rst iteration, j ¼ 1 and D is constructed using (a unscreened (u)) À1 for the respective atoms along the diagonal and zeros elsewhere, and s AB (u) is calculated from a unscreened (u). s j is dened in eqn (61) where s AB (u) is computed using P A,j and P B,j in each iteration j, D j is the matrix which has (P non-dir A,j (u)) À1 of each atom on the diagonal and the rest is zero. The new D j and s j are fed into eqn (65) for the next iteration. The calculation continues until X On the last iteration, a raw_non-dir A (u) is calculated via eqn (45). Because the polarizability upper bound is calculated for uniform applied electric eld and does not include any directional interactions between dipoles, it is applied at the end of nondirectional screening and before directional screening. For each u value, the polarizability upper bound was applied to the raw nondirectionally screened polarizability (i.e., a raw_non-dir A (u)) via the following equation to generate note that Directional screening for static polarizability. This procedure is only done at zero frequency. For directional screening, D and s (eqn (63)) are 3Natoms by 3Natoms matrices. D is a block diagonal matrix which has the inverse polarizability tensor of each atom on the block diagonal and the rest is zero. In the rst iteration, D and s AB are obtained using a force-eld A . P j undergoes a partitioned partial contraction (eqn (47)) to obtain partially screened polarizability tensors for each atom. The isotropic polarizability of each atom is 1/3 of the trace of their tensors. From this partially screened polarizability, the associated Gaussian screening width (eqn (19)) is obtained and fed into s j . The polarizability tensor of each atom is inverted to construct a new D matrix. The new D j and s j are fed into eqn (65) for the next iteration. The calculation continues until D j sums to 1. The anisotropic polarizability correction is then applied as described in the next section.
Directional screening for uctuating polarizabilities. D and s (eqn (62)) are 3Natoms by 3Natoms matrices. D is a block diagonal matrix which has the inverse polarizability tensor of each atom on the block diagonal and the rest is zero. In the rst iteration, D and s AB (u) are calculated using a non-dir A (u). P j undergoes a partitioned partial contraction (eqn (46)) to obtain partially screened polarizability tensors for each atom. The isotropic polarizability of each atom is 1/3 of the trace of their tensors. From this partially screened polarizability, the associated Gaussian screening width (eqn (19)) is obtained and fed into s j . The polarizability tensor of each atom is inverted to construct a new D matrix. The new D j and s j are fed into eqn (65) for the next iteration. The calculation continues until D j sums to 1. On the last iteration, a screened A (u) is calculated via eqn (46). This procedure is done at multiple frequencies and the resulting {a screened A (u)} are fed into the Casimir-Polder integral to calculate the C 6 dispersion coefficient of each atom. Note that

Anisotropic polarizability correction
We found that both the TS-SCS and MCLF screening processes oen overestimate anisotropy of the static polarizability tensor. Physically, this occurs because the TS-SCS and MCLF methods use spherical Gaussian dipoles in their screening models. In reality, the effective screening width should be larger along the direction of increased polarizability and smaller along the direction of decreased polarizability. Consequently, using a spherical Gaussian dipole model under-screens along the direction of increased polarizability and over-screens along the direction of decreased polarizability. This inates the polarizability tensor anisotropy.
where the correction factor (C.F.) is between 0 and 1. Because the MCLF AIM static polarizability tensors have all non-negative eigenvalues, it follows from eqn (71) that the projection of a ! ! static A =a static A onto any direction exceeds C.F. To nd a reasonable value for C.F., polarizability anisotropy ratios for four polyacenes were compared in Table 5. The reference ratios were calculated from Jiemchooroj et al.'s TD-DFT results. 78 As shown in Table 5, C.F. ¼ 0.2 gave the lowest mean relative error (MRE) and mean absolute relative error (MARE) for these materials. This value performed better than C.F. ¼ 0 and better than the TS-SCS method. Therefore, we used C.F. ¼ 0.2 as the regular value for the MCLF method. Except where otherwise specied, all reported MCLF results used this C.F. value.
As a second example, we chose a linear C 6 H 2 molecule. This molecule's structure was optimized using B3LYP functional and def2QZVPPDD basis set in Gaussian09. Optimized bond lengths from one end to another were 1.06 (H-C), 1.21 (C-C), 1.35 (C-C), 1.21 (C-C), 1.35 (C-C), 1.21 (C-C), 1.06 (C-H)Å. The calculated molecular polarizability along the bonds is 181. 15  Additional examples were provided by two test sets described in the following sections. For a test set containing 57 diatomic molecules (see Section 5.1), the MRE and MARE for eigenvalues are À6.90% and 11.77% with C.F. ¼ 0.2, and À7.96% and 13.67% with C.F. ¼ 0. For a test set containing small molecules (see Section 5.3), the MRE and MARE for eigenvalues are 5.45% and 8.10% with C.F. ¼ 0.2, and 4.85% and 9.40% with C.F. ¼ 0.

Multibody screening parameter (MBSP)
When a uniform external electric eld is applied, the atomic dipoles induced by the eld will align and atoms will interact with not only their neighbors but also atoms far away. The dispersion force, however, is caused by uctuating dipoles. The uctuating dipoles of the atoms will align with their neighbors but out of sync with atoms far away. The MBSP controls the length scale over which directional alignment persists: First-order exponential decay (eqn (73)) is the natural choice for the directional alignment function, because if 50% of the directional alignment persists over a distance d half , then over a distance 2 Â d half the expected persistence of directional alignment will be (50%) Â (50%) ¼ 25%. As dened in the ESI, † r unscreened damp,A and r unscreened damp,B are unscreened damping radii of the atoms-in-material. We optimized MBSP with the C 6 of the 12 polyacenes and 6 fullerenes, the same set studied in Section 5.5. Table 6 shows the MRE and MARE of different MBSP values of which 2.5 is optimal. Fig. 9 is the ow diagram for MCLF method. This method yields three kinds of polarizabilities: a force-eld , a static and a low_freq . a force-eld is the polarizability with no directional ordering of the electric eld and intended for use in force-eld simulations. a static is the static polarizability in a constant applied electric eld. a low_freq reects the short-range order and long-range disorder of the uctuating dipoles present in the dispersion interaction.

Flow diagram and explanation for three polarizability types
Because polarizability is a multibody effect, the polarizability of a molecule is not the sum of polarizabilities of isolated atoms. (In contrast to an atom in a material, the term isolated atom means an atom sitting alone in vacuum.) Directional interactions between atoms create components to the molecular polarizability that do not exist for the isolated atoms. 31,32 Directional dipole-dipole interactions 31,32 are considered during classical atomistic (e.g., molecular dynamics or Monte Carlo) simulations that utilize polarizable force elds. 16,95 To avoid double counting these directional dipole-dipole interactions, the force eld's atomic polarizability parameters must be the non-directionally screened values. [96][97][98][99] Many authors proposed direct additive partitioning of the quantum-mechanically computed molecular polarizability tensor into constituent atoms. 10,[100][101][102][103][104][105] Because those values already incorporate directional dipole-dipole interactions, their use in force elds would result in double counting the directional dipoledipole interactions; therefore, we do not recommend their use as force eld parameters. Our method provides both the nondirectionally screened and directionally screened AIM polarizabilities, and the non-directionally screened values should be used as polarizable force-eld parameters. The directional screening then arises during the course of the classical atomistic simulation.
Although the mixed C 6,AB dispersion coefficients could in principle be computed directly from the Casimir-Polder integral using the AIM polarizabilities at imfreqs, this would involve many integrations for materials containing thousands of atoms in the unit cell. Therefore, we used the following mixing formula which is consistent with both Padé approximation 81 and QDO 9 models: For MCLF, the polarizabilities appearing in eqn (74) must be a low_freq , because these are the polarizabilities associated with the dispersion interaction. Of note, the TS and TS-SCS methods use a similar mixing formula, except the polarizabilities appearing in the mixing formula are a TS and a TS-SCS , 12,46 because those methods do not yield a low_freq .

Results and discussion
Unless otherwise labeled, the following results are calculated using DDEC6 partitioning. The DDEC6 results were calculated using the Chargemol program. 61,[66][67][68] The MCLF results computed in this article avoided large matrix inversions. Details of the MCLF computational algorithm and computational parameters (i.e., dipole interaction cutoff length, Rhomberg integration order, and Richardson extrapolation of the screening increments) are presented in the companion article. 79 TS-SCS results in the present article were computed using matrix inversion via Gaussian elimination with partial pivoting (GEPP). The dipole interaction cutoff length for both MCLF and TS-SCS calculations was 50 bohr. Our MCLF calculations used a smooth cutoff function (eqn (64)). Our TS-SCS calculations used a hard cutoff (as consistent with prior literature 106 ). As explained in the companion article, Rhomberg integration for both MCLF and TS-SCS used 16 imfreq points. 79 For MCLF, Richardson extrapolation of the screening increments used the number of steps recommended in the companion article. 79

Diatomic molecules
We rst tested our model's sensitivity to the choice of exchangecorrelation functional and basis set used to compute the system's electron and spin density distributions. A set of 57 diatomic molecules with elements across the periodic table were chosen as the test set. Three sets of electron density distributions of the test set were generated using Gaussian09 with CCSD/def2QZVPPDD, Gaussian09 with B3LYP/ def2QZVPPDD, and VASP with PBE/planewave. The Chargemol program was then used to DDEC6 partition the electron density followed by MCLF analysis to obtain the static polarizabilities and C 6 coefficients. Fig. 10 shows polarizability and C 6 computed using CCSD/def2QZVPPDD densities versus polarizability and C 6 computed using the PBE/planewave and B3LYP/ def2QZVPPDD densities. This gure shows our model did not show strong dependence on the choice of exchange-correlation functional or basis set. Hence, MCLF gives similar results using electron density distributions from different procient quantum chemistry levels of theory.
For the same 57 diatomic molecules, the isotropic static polarizability and the three eigenvalues of the static polarizability tensor from TS-SCS and MCLF were compared to the reference data. The reference is CCSD calculations with def2QZVPPDD basis set. System-specic polarizabilities for CCSD reference, TS-SCS, and MCLF are listed in the ESI. † Table  7 summarizes the error statistics. The TS-SCS method gave large errors independent of the charge partitioning method used (Hirshfeld or DDEC6). On average, MCLF was four times more accurate than TS-SCS for these materials. Fig. 11 is the absolute percentage error of isotropic polarizability from TS-SCS versus DDEC6 NAC magnitude for these diatomic molecules. From the plot, we can see TS-SCS gives large errors when the NAC magnitude is high. This conrms the TS-SCS model is less accurate for charged atoms and MCLF xed this problem.

Dense periodic solids
Static polarizabilities were computed for 28 dense periodic solids including electric insulators, semi-conductors, and conductors. The geometries are from Inorganic Crystal Structure Database (ICSD). 107,110 We generated electron densities in VASP 108 using the PBE 109 functional. See Gabaldon-Limas and Manz 68 for a description of the VASP computational settings used.
Since the polarizability should not exceed the conduction limit, the reference static polarizability was set equal to the smaller value of the Clausius-Mosotti relation (eqn (75)) and conduction limit (eqn (76)): where volume is obtained from the ICSD 107,110 crystal structure and k is the experimental static dielectric constant. Results from different partitioning methods and screening methods were compared to this reference data. ICSD codes and computed results for individual materials are listed in the ESI. † As listed in the ESI, † the experimental dielectric constants were taken from Young and Frederikse 111 and other sources. Table 8 summarizes the error statistics. Comparing the MRE and MARE of the screened and unscreened polarizabilities with the same partitioning method (e.g., TS/H vs. TS-SCS/H), screening increases the accuracy for each partitioning method. Because all of the unscreened methods gave >100% error for some materials, screening is an essential step in polarizability calculations. Comparing TS-SCS with different partitioning methods: IH is more accurate than H, and DDEC6 is more accurate than IH. The tendency of TS-SCS/H and TS-SCS/IH to overestimate polarizabilities for solids was previously reported by Bucko et al. 57 with improved results reported using the fractionally ionic approach of Gould et al. 54 (Those studies included some of the same materials as here. 54,57 ) Using DDEC6 partitioning, MCLF gives  more accurate results than TS-SCS with MARE of 12% and 24%, respectively.

Polarizabilities of small molecules
A set of 22 molecules were selected from Thole 32 and Applequist et al. 31 that have experimentally measured isotropic static polarizability. Six of these also have experimentally measured polarizability tensor eigenvalues. 31,32 The geometries are from geometry optimization we performed in Gaussian09 with B3LYP functional and def2QZVPPDD basis set. Tables 9 and 10 show that both TS-SCS/DDEC6 and MCLF performed well for this test set.

C 6 coefficients of atom/molecule pairs
This test involves C 6 coefficients for pairs of atoms and molecules studied by Tkatchenko and Scheffler. 12 Our geometries are from geometry optimization we performed in Gaussian09 with B3LYP functional and def2QZVPPDD basis set. Fig. 12 compares the TS-SCS/DDEC6 and MCLF C 6 coefficients to the reference values derived from the dipole oscillator strength distribution (DOSD) data of Meath and co-workers 74,75 as tabulated by Bucko et al. 55 As shown in Fig. 12, TS-SCS predicts too large C 6 values for the larger molecules, while MCLF is consistently accurate for different sized molecules. For these 49 atoms/molecules, MCLF gave 1.15% MRE and 5.79% MARE while TS-SCS/DDEC6 gave 8.09% MRE and 10.29% MARE.
The same reference data source was also used for the 1225 pairs formed from these 49 atoms/molecules. Fig. 13 plots MCLF versus reference C 6 values for these pairs. These MCLF C 6,AB values were computed from the C 6,A values using eqn (74). Then, the molecular C 6 were computed from these C 6,AB using eqn (1 113 Other studies also investigated this test set.

Polyacenes and fullerenes
Polycyclic aromatic compounds, such as polyacenes, and fullerenes, have strong directional alignment of induced and   uctuating dipoles in the ring planes. A set of 12 polyacenes and a set of 6 fullerenes were selected as test sets. We generated electron densities in VASP 108 using the PBE 109 functional. See Gabaldon-Limas and Manz 68 for a description of the VASP computational settings used.
The polyacene geometries are from our VASP geometry optimization using the PBE functional. Polyacene reference static polarizabilities and C 6 coefficients are from Marques et al. 77 The polyacene structures are shown in Fig. 14. Computed results are summarized in Table 11. The MRE and MARE show MCLF was signicantly more accurate than TS-SCS for describing both the static polarizabilities and the C 6 coefficients of these materials. Table 12 summarizes calculation results for fullerenes. The fullerene geometries are from Saidi and Norman. 114 Tao et al. studied this set of fullerenes in 2016 and obtained excellent results using a hollow sphere model with modied single frequency approximation. 115 The reference static polarizabilities and C 6 coefficients are from Kauczor et al.'s TD-DFT calculations. 116 For this test set, TS-SCS systematically underestimates the polarizabilities by 18.7% and C 6 coefficients by 10.4%.
In these two test sets, MCLF has better overall performance than TS-SCS with all four MCLF MAREs under 10% compared to all four TS-SCS MAREs over 10%. In contrast to TS-SCS, MCLF uses (i) an iterative update of the Gaussian dipole width and (ii) a multi-body screening function (eqn (73)) to describe decay of the uctuating dipole directional order. These allow MCLF to describe induced and uctuating dipole directional alignment effects more accurately than TS-SCS.

Large biomolecule
Non-reactive molecular mechanics force elds are model potentials containing several different parameter types: dispersion-repulsion parameters (e.g., Lennard-Jones, Buckingham, or other forms), point charges or model atomic charge distributions (e.g., Gaussian, Slater), exibility parameters, and (optionally) polarizabilities. 117,118 These molecular mechanics force elds enable classical atomistic simulations, such as molecular dynamics or Monte Carlo, to be performed over larger distance and time scales than would be practical using quantum chemistry methods such as DFT. [118][119][120] These atomistic simulations are useful to estimate thermodynamic ensemble properties (e.g., density, vapor pressure, adsorption isotherms, etc.), transport properties (e.g., diffusion coefficients, viscosity, etc.), and structures (e.g., protein folding and other conformational changes). 118,119,[121][122][123][124] An approach oen used is to classify atoms into types, where similar atoms share the same atom type and same force-eld  parameters. 125,126 We will refer to these as Typed Force Fields (TFF). TFFs have been successively improved over the past few decades by rening their atom type denitions and parameter values to make them more accurate and robust. [127][128][129] Today, several TFFs perform reasonably well on various organic molecules and some small inorganic molecules. 122,123,[130][131][132][133] There are still areas for further improving force elds, especially for systems containing high chemical bonding diversity and charged ions. Metal-containing systems are especially prone to high chemical bonding diversity. For example, metal-organic frameworks (MOFs) can contain dozens of different metal elements in a plethora of different bonding motifs. [134][135][136] While some efforts have been made to dene new atom types for MOFs, 137 the high chemical diversity makes it difficult to completely parameterize force elds for all MOFs using atom types. Quantum-mechanically derived force-elds (QMDFFs) are ideal for modeling these systems, because QMDFFs do not require pre-dened atom types. 138,139 Machine learning is a recent approach in which force-eld parameterization is trained to a machine learning model using QM-derived parameters. 4,140,141 Advantages of the machine learning approach include high automation, fast computation, and the ability to handle large chemical diversity without having to manually dene atom types. 4,140,141 Non-polarizable force elds oen model charged ions with reduced effective charges, but this places articial constraints on the simulation. 142 To make the parameters more transferable between different chemical systems and compositions, the reduced effective charge model should be replaced with a polarizable force-eld. 16,19,[143][144][145] Kiss and Baranyai concluded for water that "It is impossible to describe the vapor-liquid coexistence properties consistently with nonpolarizable models, even if their critical temperature is correct." 15 Hence, an automated method like MCLF to assign atom-in-material polarizabilities is extremely important for modeling materials containing ions.
Useful insights can be gained by comparing TFF parameters for an atom type to QM-derived ones. Where the TFF parameters and the QM-derived ones are in good agreement, this validates the parameterization. Conversely, where the TFF parameters and the QM-derived ones differ substantially, this indicates areas for further study to potentially rene the atom type denitions and/or parameter values. A multimodal distribution or wide range of QM-derived parameter values suggests when to divide atoms into multiple atom types. A narrow distribution of QM-derived parameter values that is substantially offset from the TFF parameter value can indicate a need to update the TFF parameter value. In such a way, these comparisons can produce force-eld improvements.
In this section, we compare MCLF C 6 atom-in-material dispersion coefficients and polarizabilities to OPLS and AMOEBA force-eld parameters, respectively, for the Human Immunodeciency Virus reverse-transcriptase (HIV-RT) enzyme complexed with an inhibitor molecule. HIV is a retrovirus with RNA genome. 146 Retroviruses replicate in a host cell by using a reverse-transcriptase enzyme to transcribe the virus's RNA genome into the host cell's DNA that is subsequently replicated by the host. 146 Therefore, inhibition of the reverse-transcriptase enzyme is a potential way to slow virus replication, which is extremely important for controlling disease caused by the virus. 146 Bollini et al. and Cole et al. previously studied several oxazole derivatives as HIV-RT inhibitors. [147][148][149] As an example of the MCLF method applied to macromolecules, we study one of these HIV-RT inhibitors, C 20 H 16 ClF 2 N 3 O (CAS # 1422256-80-1), shown in Fig. 15 together with a signicant portion (2768 atoms) of its complex with wild-type HIV-RT. We computed the electron density distribution for this complex in ONETEP 150,151   using the PBE 109 exchange-correlation functional. Preparation of the input structure is described elsewhere. 147,148 It was constructed from the 1S9E PDB le 152 using the MCPRO 153 and BOMB 154 soware. The 178 amino acids closest to the ligand were retained. The complex was solvated in a 25Å water cap and equilibrated at room temperature for 40 million Monte Carlo steps using the MCPRO soware. Water molecules were stripped from the nal conguration, and the resulting structure ( Fig. 15) was used as input for the ONETEP calculations.
Interactions between electrons and nuclei were described by Opium norm-conserving pseudopotentials. 155 NGWFs were initialized as orbitals obtained from solving the Kohn-Sham equation for isolated atoms. 156 The NGWFs were expanded as a psinc basis set 157 with an equivalent plane-wave cutoff energy of approximately 1000 eV and the electron density was stored on a Cartesian grid of spacing 0.23 Bohr. The localization radii of the NGWFs were 10.0 Bohr. Calculations were performed using an implicit solvent model with a dielectric constant of 78.54 to mimic the water environment. 158 Electron density partitioning was performed using the DDEC6 method implemented in the Chargemol program. 61,66-68 Table 13 compares the computed C 6 coefficients in the HIV-RT complex for a number of frequently-occurring OPLS atom types, 122 including both backbone and side chain atoms. The C 6 coefficients computed using MCLF generally show a much smaller range than the corresponding TS-SCS/DDEC6 data. Also, the OPLS force eld C 6 coefficients are usually closer to the MCLF results than to the TS-SCS results. While one should not draw too many conclusions from this, the OPLS parameters have been carefully t over a number of decades to Table 15 The MCLF method is a set of guiding principles and their empirical implementations. The first eight guiding principles are physically motivated. The last three guiding principles are motivated by computational efficiency

Guiding principle Empirical implementation
Method should require reference polarizabilities and C 6 coefficients for neutral free atoms only (not charged free atoms) and still properly describe scaling laws for charged atoms a Isolated atom scaling laws (eqn (30) and (31)) Conduction limit upper bound for a material's non-directional polarizability at xed geometry and xed orientation Conduction limit upper bound function (eqn (52)) applied to nondirectional polarizabilities as a function of frequency (eqn (68)) using AIM volume (eqn (50)) Polarizability interaction not necessarily distributed equally between a pair of atoms  accurately reproduce experimental observables, such as organic liquid properties 122 and protein NMR measurements. 159 It is also noteworthy that the OPLS C 6 coefficient is slightly higher than the MCLF result. This is expected since the C 6 term in the OPLS force-eld must effectively compensate for higher-order dispersion (C 8 , C 10 ,.) terms that are not explicitly included in the OPLS force-eld. 160 The high value of C 6 on the backbone carbonyl carbon (C (bb)) has been noted previously, 3 and should be revisited in future force elds. Table 14 compares the three kinds of MCLF polarizabilities (i.e., force-eld, static, and low_freq) to parameters used in the AMOEBA 16,144 polarizable force eld. For the H atoms, the AMOEBA polarizabilities are larger than the MCLF polarizabilities, though both are small compared to the polarizabilities of C, N, and O atoms. For the N (bb) atom, the AMOEBA polarizability is similar to the MCLF force-eld polarizability. For the O (bb) atom, the AMOEBA polarizability (5.6) is smaller than the MCLF polarizability (6.9). Since the polarizability of an isolated neutral oxygen atom is $5.2, 65,80 these polarizabilities are in line with O (bb) being slightly negatively charged and more diffuse than an isolated neutral oxygen atom. For the C atoms, the AMOEBA polarizabilities are substantially larger than the MCLF force-eld polarizabilities and close to the MCLF static polarizabilities; this suggests the AMOEBA C atom polarizabilities include some directional screening effects. The AMOEBA polarizabilities use Thole 32 model atomic charge distributions 16 c 1 exp(Àc 2 (r A ) 3 ) while the MCLF polarizabilities use Gaussian model charge distributions, so the optimized atomic polarizabilities can be different in these two methods and still approximately reproduce the molecular polarizability. Importantly, computed results presented earlier in this article show the MCLF force-eld polarizabilities input into the dipole interaction tensor and solved yield good accuracy molecular static polarizabilities. Therefore, the MCLF force-eld polarizabilities are appropriate for use in polarizable force-elds.
OPLS atom types 122 were assigned using MCPRO soware 153 and have the following descriptions. Sidechain atom types included: CT(CH 3 ): alkane carbon bonded to three hydrogen atoms; CT(CH): alkane carbon atom bonded to one hydrogen atom; HC: alkane hydrogen atom; CA: carbon atom in an aromatic ring; HA: hydrogen bonded to an aromatic ring. Backbone atoms were classied according to: N: amide nitrogen; H: hydrogen bonded to amide nitrogen; C: amide carbon; O: amide oxygen; CT: alpha carbon atom. Only the most common atom types were analyzed to ensure adequate statistics. The substantial ranges of QM-derived parameter values in Tables 13 and 14 suggest there is room to further improve the atom type denitions. An alternative approach denes atom types as the set of unique local environment subgraphs that extend up to one or two bonds, and this approach has been successfully used with DDEC NACs. 4,68,161 We recommend the atom type denitions be explored further in future work.
AIM polarizabilities and higher-order AIM dispersion coefficients (e.g., C 8 , C 9 , C 10 ) and associated mixing rules. (8) The MCLF atom-in-material polarizability tensors are always symmetric, while the TS-SCS atom-in-material polarizability tensors are sometimes asymmetric. Symmetric polarizability tensors are more convenient, because they can be displayed as ellipsoids. (9) As explained in the companion article, the computational cost of MCLF scales linearly with increasing number of atoms in the unit cell for large systems. 79 This is achieved using a dipole interaction cutoff function combined with computational routines that avoid both large matrix inversions and large dense matrix multiplications. 79 Tests were performed on diverse material types: isolated atoms, diatomic molecules, periodic solids, small organic and inorganic molecules, fullerenes, polyacenes, and an HIV reverse transcriptase biomolecule. For each test set in this study, MCLF gave #12% MARE on the static polarizabilities, C 6 coefficients, and static polarizability eigenvalues. This substantially improves over the TS-SCS method. For the static polarizabilities of solids: (a) TS-SCS with H, IH, and DDEC6 partitioning gave MARE of 47%, 30%, and 24%, respectively, (b) MCLF gave MARE of 12%, and (c) all of the unscreened methods gave much larger errors than the screened methods. For the static polarizabilities of diatomic molecules: (a) TS-SCS/H gave 41% MARE with a largest error of 437%, (b) TS-SCS/DDEC6 gave 37% MARE with a largest error of 440%, and (c) MCLF gave 10% MARE with a largest error of 34%. We anticipate MCLF should be useful for parameterizing polarizable force elds and DFT + dispersion methods.
There are several key differences between the MCLF and D3 and D4 approaches. The MCLF and D4 (ref. 35 and 37) approaches incorporate atomic charge information, while the D3 (ref. 11) method does not. The D3 and D4 approaches are based on the molecular geometry without requiring a quantum-mechanically computed electron density distribution, 11,35,37 while MCLF uses a quantum-mechanically computed electron density distribution. (The atomic charges used in D4 method variants can be computed from a quantum-mechanically computed electron density distribution, but this is not required. 35,37 ) The MCLF method explicitly considers dipole-dipole tensor interactions, while the D4 method only does so when a damped MBD Hamiltonian is included. (Without the MBD Hamiltonian, some dipole-dipole tensor interactions may be implicitly included in the D3 and D4 methods via coordination number effects, but only to the extent those dipole-dipole interactions correlate to the atom's coordination number and resemble dipole-dipole interactions in the reference compounds.) When using classical electronegativity equilibration NACs, the D4 method facilitates computing analytic forces during DFT + dispersion calculations 37 MCLF yields AIM and system polarizability tensors, while D4 yields isotropic AIM and system polarizabilities. MCLF yields three types of polarizabilities: (a) nondirectionally screened polarizabilities suitable for use in polarizable force-elds, (b) uctuating polarizabilities that describe London dispersion interactions, and (c) static induced polarizabilities that include directional dipole-dipole interactions under a constant externally applied electric eld. Table 15, we believe it is most useful to regard the MCLF method as a set of guiding principles and their empirical implementations. Some of the guiding principles are non-empirical (i.e., derived from rst-principles) while their implementation involves some empirical aspects (i.e., tting functional forms to observational data). Both the non-empirical and the empirical aspects of the MCLF method are highly innovative and important. Therefore, we do not believe it is reasonable to categorize the entire MCLF method as being exclusively empirical or exclusively non-empirical. If one had to choose a single word to categorize it, "semiempirical" is the appropriate choice. The Merriam-Webster dictionary denes semiempirical as "partly empirical". 162 Notably, one could generate alternative empirical implementations of the same guiding principles. For example, the isolated atom scaling laws could be reconstructed using the hr 2 i and hr 4 i moments in place of the hr 3 i and hr 4 i moments, as shown in Table 1. This implementation change would usually produce minor differences in results, because both models are trained to the same underlying dataset of isolated atom polarizabilities and C 6 coefficients. As another example, using a different but similar smooth minimum function to impose the conduction limit upper bound could usually produce similar results, because both functions are ultimately controlled by the same physical limit of the non-directional polarizability of an ideal conductor at xed geometry and xed orientation. Furthermore, the MCLF method must be applied using some chosen partitioning method (e.g., DDEC6). Clearly, some partitioning methods are poor choices (e.g., Hirshfeld as shown in Table 8), while others (e.g., DDEC6) are good choices. The DDEC6 method is also a set of guiding principles (which contain some non-empirical aspects) and their empirical implementations that comprise a semiempirical method. 61,68

Conflicts of interest
There are no conicts to declare.