Describing Molecular Polarizability by a Bond Capacity Model

Subscriber access provided by OCCIDENTAL COLL

Molecular Mechanics

Describing Molecular Polarizability by a Bond Capacity Model. Pier Paolo Poier, and Frank Jensen J. Chem. Theory Comput., Just Accepted Manuscript • Publication Date (Web): 28 Mar 2019 Downloaded from http://pubs.acs.org on March 28, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Describing Molecular Polarizability by a Bond Capacity Model. Pier Paolo Poier, Frank Jensen* Department of Chemistry, Aarhus University, Langelandsgade 140, DK-8000 Aarhus, Denmark

ABSTRACT

We propose a bond capacity model for describing molecular polarization in force field energy functions at the charge-only level. Atomic charges are calculated by allowing charge to flow between atom pairs according to a bond capacity and a difference in electrostatic potential. The bond capacity is closely related to the bond order and decays to zero as the bond distance is increased. The electrostatic potential is composed of an intrinsic potential, identified as the electronegativity, and a screened Coulomb potential from all other charges. The bond capacity model leads to integer fragment charges upon bond dissociation and displays linear scaling of the polarizability with system size. Bond capacity parameters can be derived from reference molecular polarizabilities, while electronegativity parameters can be derived from reference atomic charges or a reference molecular electrostatic potential. Out-of-plane polarization for planar systems is modelled by off-nuclei charge sites. The model is shown to be able to reproduce anisotropic reference molecular polarizabilities with an accuracy of ~10% using a limited set of bond capacity parameters, and can describe both inter- and intra-molecular polarization.

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Introduction. Computational efficient functions for describing the geometry dependence of molecular energies are key ingredients for molecular simulations.1-3 Ab initio energy functions, typically of the density functional theory (DFT) type, can be used for simulating small systems for pico-seconds, but simulations of large biomolecular systems for nano- or micro-seconds require parameterized energy functions, usually denoted force fields. Common examples are CHARMM,4 AMBER5 and GROMOS6 which have been developed and refined over the last three decades and are widely used for modelling biomolecular systems. The bonded energy terms in these force fields are taken as harmonic expressions or low-order Fourier series, while the non-bonded energy is partitioned into a non-polar van der Waals term, usually taken as a Lennard-Jones potential, and a polar electrostatic term, described as fractional atomic charges interacting by a Coulomb expression. In some of the more recent versions, the electrostatic energy has been augmented with polarization terms, by means of atomic dipole polarizability in AMBER7 and the Drude particle model in CHARMM,8 with an earlier version employing the fluctuating charge model.9 The AMOEBA is an example of a modern force field where atomic electrostatic moments up to quadrupoles and dipole polarizability have been included from the onset, and parameterized against electronic structure results.10 The electrostatic energy is the major component in a force field energy for biomolecular systems, and is the only part considered in the present work. Electronic structure methods are widely used for parameterizing force fields, as there rarely are sufficient experimental data available for assigning the often very large number of parameters. Atomic charges can be assigned based on decomposition schemes for the electronic wave function, for example by the distributed multipole analysis (DMA)11 or LoProp12 methods, or from fitting to the molecular electrostatic potential (ESP),13-14 and these methods can be extended to atomic


2

Page 2 of 59

Page 3 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dipole and quadrupole moments. Atomic charges can alternatively be assigned based on partitioning the electron density into atomic contributions, and variations of the Hirshfeld procedure15 has been suggested for deriving force field electrostatic parameters.16-19 It is well known that atomic charges, and atomic electrical moments in general, depend on the molecular geometry and this can be considered as arising from two different underlying effects, changes in orbital overlaps and an electric field effect. The former is the dominating factor for the dependence of atomic charges on changes in bond lengths and valence angles,20 where the nature of the wave function changes due to changes in atomic orbital overlaps. The electric field effect, usually denoted polarization, occurs between moieties that are sufficiently far removed that wave function overlap can be ignored. The geometry dependence is in this case mediated by changing the electric potential/field, which in turn changes the wave function. At intermediate distances, where there is a small orbital overlap, the polarization effect can include a charge transfer component, and this is important for describing for example hydrogen bonding. The atomic charge dependence on bond distances and valence angles is a key component for modelling intensities in infrared spectroscopy,21-22 while intermolecular polarization is essential for accurately modelling intermolecular interactions. The conformational dependence of atomic charges due to changes in torsional angles is mainly an intramolecular polarization effect. Modelling the force field electrostatic energy requires a choice for assigning atomic charges and possibly higher order moments. The electron density, and thus the molecular ESP, is a physical observable that can be calculated with a good accuracy using medium level electronic structure methods, and reproducing the ESP is an objective quality measure. The DMA and LoProp methods require inclusion of up to atomic quadrupole or octupoles moments for an accurate reproduction.23 A lower error can be obtained by explicit fitting parameters to the ESP but we have shown that the


3


parameter space rapidly becomes redundant,24 and there is a large flexibility for choosing a nonredundant set of parameters, i.e. a specific combination of charge, dipole and quadrupole moments.25 The parameter redundancy is likely a component in the observation that atomic charges derived by ESP fitting often display significant dependence on the molecular geometry. The electron density can be partitioned into atomic contributions based on an integration with an atomic weight factor at each point in physical space. If the weight factor is restricted to values of 0 and 1, the electron density at a given point belongs to only one atom, and the spatial partitioning can be denoted hard-boundary. If the weight factor can take any value between 0 and 1, the electron density at a given point is distributed between all (in principle) atoms, and the spatial partitioning can be denoted soft-boundary. A number of closely related methods based on soft-boundary partition of the electron density into atomic contributions has been proposed for assigning less geometry dependent charges that provide a good representation of the ESP.26-27 The original Hirshfeld method15 produces atomic charges that are too small for modelling the ESP, but charges derived by the iterative Hirshfeld,28 the iterative Stockholder29 and their variations26, 30-32 are often found to provide ESP close to that obtained from ESP fitted charges.27 The Hirshfeld/Stockholder approaches, however, are not without problems,33 and for example rely on using artificial (basis set) bound electron densities for anions34 and restricted function choices for fitting atomic densities.26, 30 The DDEC (density derived electrostatic and chemical) family of methods can be considered as combinations of the iterative Hirshfeld and Stockholder approaches, where advantages and disadvantages of the underlying methods are inherited in the combination.35-37 The CMx family of atomic charge models is a combination of a wave function based atomic charges and an empirical function containing atomic and diatomic parameters for a better reproduction of the molecular dipole moment.38 Atomic charges derived by hard-boundary partitioning of the


4

Page 4 of 59

Page 5 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


electron density, such as quantum theory of atoms in molecules, are in general not suitable for force fields, unless a number of higher order multipole moments are also included.39 The ESP derived charges are for planar systems identical to the out-of-plane component of the atomic polar tensor, which determines the intensities in infrared spectroscopy within the double harmonic approximation.22, 40-41 The generalized atomic polar tensor charges defined as 1/3 of the trace of the atomic polar tensor42 are known to be poor at reproducing ESP, since they include average charge-flux contributions, which usually are significantly larger than the static contribution.22, 40, 43

Given that the parameter space for the static electrical moments rapidly becomes redundant,24 it is very likely that the corresponding polarization parameter space will also contain redundancies even at low orders. Non-polarizable force fields employing (only) atomic charges have been the workhorse for performing simulations for several decades but improved force fields need to consider both higher order multipole moments and polarization simultaneously, as these effects are of similar importance.24 The question is thus how to choose in a systematic and efficient fashion a consistent set of parameters that provide a predetermined accuracy as measured relative to a set of reference data. From a simple order expansion argument, it is expected that static atomic charge parameters should be more important that atomic dipole and quadrupole moments, and that charge polarization should be more important than dipole polarization, and the latter has been shown explicitly.44 Unke et al. have recently shown that an optimized distributed charge model with only a few off-nuclei charges is better at reproducing the ESP than a full atomic charge-dipolequadrupole model,45 and this emphasizes the need for a model that can include polarization at the charge-only level. Charge polarization has in previous work been modelled by the fluctuating charge method or the corresponding split-charge formulation discussed in the next section, but


5


Page 6 of 59

these models include non-physical effects such as charge transfer over infinite distance and nonlinear scaling of the polarizability with system size. The Drude particle and point dipole polarization methods aim at describing atomic dipole polarization, and do not account for charge polarization, while advanced force fields such as NEMO include also quadrupole polarizabilities.46 Force fields based on continuous representations of orbital electron densities, such as SIBFA,47 GEM,48 EFP49 and RexPoN50-51 are computationally significantly more demanding. In the present work, we propose a bond capacity model that is capable of describing charge polarization within a force field model, and that, in contrast to the fluctuating charge model, leads to integer fragment charges upon bond dissociation and displays linear scaling of the polarizability with system size. The bond capacity model is capable of reproducing molecular polarizabilities from electronic structure calculations with an accuracy of ~10% with a limited set of parameters. It provides a general framework for modelling electrostatic interactions and can be parameterized against polarizabilities and a selected set of atomic charges, of which some are discussed above, or directly against a molecular ESP. We will in the present case focus on the general framework, with illustrations for selected systems, while the use of the model for designing new force fields will require a careful selection of atom types and corresponding parameters, and reference data. We will in the next section review the fluctuating charge and split charge models for placing the present bond capacity model in context. Previous work, fluctuating charge and split-charge models. Several closely related methods have been proposed for modeling charge flux and charge polarization effects: fluctuating charge (FQ) or charge equilibration,9,

52-53

extended charge

equilibrium and its bond-order corrected version (EQeq and EQeq+C),54-55 electronegativity equalization method,56 atom-atom current equilibration (ACTT),57 split charge equilibration


6

Page 7 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(SQE),58-60 charge transfer with polarization current equilibration (QTPIE),61-62 chemical potential equilibration63-64 and atom-condensed Kohn-Sham DFT model approximated to second order (ACKS2).65 Some of these models are formally identical but formulated in different parameter spaces (atomic vs. bond parameters), and different acronyms are used for models differing only in the parameterization. The central idea in all these methods is a second order expansion of the electrostatic energy in terms of atomic charges Q. Eel   i Qi  12 ij Qi Q j i

(1)

ij

t

Eel  Q χ 

1 Qt ηQ 2

Since the atomic charge is directly related to the number of electrons assigned to the atom, the linear and quadratic terms for a given atom can be associated with the DFT concepts of electronegativity and hardness, respectively, being first and second derivatives of the energy with respect to the number of electrons.66-67 The energy derivatives within a DFT framework are discontinuous at integer number of electrons, which is not displayed by eq. (1), and eq. (1) is usually considered simply as a Taylor expansion with the electronegativity i and diagonal hardness ii elements being free parameters. The off-diagonal elements in the hardness matrix describe the interaction between charges and must approach 1/Rij at long distance to recover the Coulomb interaction, and can be combined with suitable screening/interpolation functions for intermediate distances.52, 59, 68-73 The equilibrium charges are obtained by making the energy stationary with inclusion of a Lagrange multiplier to ensure charge conservation.74



L  Q,    Qt χ  12 Qt ηQ   Qt 1  Qtot



Assuming that the hardness matrix is non-singular, eq. (2) can be solved explicitly.74


7

(2)


Q   η1  χ   1



1t η1χ  Qtot 1t η11

Page 8 of 59

(3)

(4)

The FQ model with constant electronegativity and hardness parameters has two non-physical implications: a non-zero charge transfer (CT) at infinite distance for atoms of different electronegativities61-62 and a cubic scaling of the polarizability with system size.57, 75 The Lagrange multiplier in eq. (3) can be interpreted as a shift in the electronegativity scale to ensure the correct total charge. For a system with fragments separated by a large distance, the -1 part of eq. (3) will be block diagonal and atomic charges in a given fragment will only depend on  and  elements within the fragment. The Lagrange multiplier in eq. (4), however, contains contributions from all atoms, regardless of interatomic distances, and it thus introduces a global coupling of all electronegativity parameters in eq. (3). A non-zero value for the Lagrange multiplier therefore leads to CT over infinite distance, and a key point is thus to design models that ensure the correct total charge without the need for additional constraints. The problem with CT at infinite distance in the FQ model has lead a number of groups to reformulate the model in terms of pairwise parameters associated with transferring charge between atom pairs, which makes it (formally) possible to restrict CT between atoms. The SQE formulation parameterizes the atomic charges in terms of pair-wise increments, often called split-charge parameters, where pji = -pij by symmetry.

Qi   p ji

(5)

j

The antisymmetry of pji means that the sum of atomic charges necessarily is zero, and this can therefore only be used for neutral systems. Extensions to charged systems have been done by eq.


8

Page 9 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(6), which spreads the total charge uniformly over all atoms,76 or by the SQE+Q0 model in eq. (7) that allows for assigning (usually integer) fixed net charges to atoms.17

Qi 

Qtot   p ji N atom j

Qi  Qi0   p ji

(6)

(7)

j

The SQE+Q0 model is formally equivalent to the ACKS2 model,65, 77 and as shown in ref. 77, is identical to the FQ model, which is also implied by the general equivalence of quadratic models discussed by Chen et al. (see eq. (9) below).78 The FQ parameterization has the advantage that only single-index (atom type) parameters are required (linear in the number of atom types), while the SQE parameterization is in terms of double-index (pairwise atom types) parameters (quadratic in the number of atom types). The splitcharge parameters can be non-zero for all atom pairs (thereby allowing CT), but are usually restricted to be non-zero only for bonded atom pairs, or made distance dependent, which effectively removes the parameters between distant atoms. Depending on the system, there may be more or fewer split-charge parameters than atomic charges (non-cyclic systems have N-1 bonds, monocyclic systems have N bonds, bicyclic systems have N+1 bonds, etc. up to a maximum of ½N(N-1) unique pair-wise terms for a fully connected system). The transformation between the two sets of parameters can be written as in eq. (8), where p is a vector containing the unique splitcharge parameters and T contains the transfer parameters (summation in eq. (5)) with elements 0, 1, i.e. T includes the anti-symmetry and p contains only the non-redundant parameters.77-78 Q  Tp

 

T1Q  Tt T

1

Tt Q  p


9

(8)


Page 10 of 59

Chen et al. have shown that for a fully connected system, any quadratic model in p-space can be transformed into Q-space (and vice versa).78 EFQ  Qt χ  12 Qt ηQ ; ESQE  pt L  12 pt Mp L  Tt χ ; M  Tt ηT

 

(9)

   

t

t

χ  T1 L ; η  T1 M T1

Since the number of Q and p parameters is not necessarily the same, the inverse transformation involves a generalized inverse matrix, eq. (8). For a non-cyclic bond-only connected system, the TtT matrix is non-singular, but since T is rectangular, T-1 is a generalized inverse. When the number of parameters is the same, T is a square matrix, but one of the eigenvalues of TtT is zero (corresponding to net charge of zero), and thus also involves a generalized inverse. When there are more split-charge parameters than atoms, there will be additional zero eigenvalues of the TtT matrix. The generalized inverse matrices in eq. (9) make the analysis of the relationship between SQE and FQ models non-trivial, since neither T-1T nor TT-1 are unit matrices. The transformation between the Q and p parameters, however, is linear, and one set can thus always be written as a linear combination of the other. The structure of the T matrix is easy, since it is given by the connectivity, but the structure of the T-1 matrix connecting a p element with Q is non-intuitive. In a fullyconnected system where all off-diagonal T-elements are 1, the connection can be written as pij = (Qi – Qj)/Natom, but in the general case, a given pij element is a linear combination of (all) Qi, i.e. a non-local connection. Discussions of split-charge models are with few exceptions limited to non-cyclic bond-connectedonly systems where there is a one-to-one mapping of the p and Q parameters. For all other systems,


10

Page 11 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


however, there are more split-charge parameters than atom parameters, and the actual number depends on the connectivity. A direct substitution of eq. (5) into eq. (1) gives eq. (10), where the anti-symmetry of pij has been utilized and the factor of ½ allows an unrestricted summation over both variables. Eel 

1 2

  i   j  pij  12 ij  pki plj ij

ij

(10)

kl

Various papers employing models within the SQE formalism usually recognize that this is just a reformulation of the FQ model, and thus should give the same result and have the same problems, but the SQE parameterization is convenient for introducing additional restrictions and terms to alleviate the FQ problems. The split-charge formulation (formally) allows restrictions on which atom pairs that are allowed to exchange charge by only allowing certain pij elements to be nonzero, and thereby controlling CT and the magnitude of the polarization. A common adjustment that is usually implied by the SQE acronym, is shown in eq. (11) where the energy expression is augmented with terms dependent on bond electronegativity ij and hardness ij parameters.58-59, 77 The split-charge parameters are analogous to the FQ model determined by minimizing the energy and subsequently using eq. (5) for determining the atomic charges.



Eel   i  pki  12 ij  pki plj   ij pij  12 ij pij2 i

k

ij

kl

ij



(11)

The last term in eq. (11) ensures that the polarizability for linear chains display the correct linear dependence, while the third term only serves to improve the fitting by making the electronegativity dependent on the neighbors.76 The bond electronegativity ij is taken as an electronegativity difference parameter, and required to be anti-symmetric in the indices. An alternative formulation shown in eq. (12) is a mixed parameterization in both charges and split-charges,59 where the FQ model is augmented with a quadratic split-charge term.


11


Eel   i Qi  12 ij Qi Q j  12 ij pij2 i

ij

Page 12 of 59

(12)

ij

Mathieu has suggested to make the bond hardness parameters ij geometry dependent, such that the last term effectively is an energy penalty for transferring charge across distance.79 Using eq. (9), the bond electronegativity parameters ij in eq. (11) can be decomposed into atomic contributions according to the linear coefficients in the T-1 matrix. The bond hardness parameters

ij in eqs. (11) and (12) similarly modify the diagonal and off-diagonal hardness elements in the FQ model upon back-transformation to Q-space, and this alters the Coulomb interaction at long distance. Including SQE penalty terms as in eqs. (11) and (12) thus modifies the Coulomb interaction between atoms in the FQ parameterization, which is incorrect in the long-distance limit. It should be further noted that adding penalty terms or setting specific pij elements = 0 only has an effect when there are more penalty terms than redundant p-parameters (assuming that each atom has at least one connection, and the whole system is connected). Most SQE papers only discuss the non-charged linear-connected case, where there are only N-1 (non-redundant) pij parameters. In the general case where the split-charge parameterization is redundant, however, penalty terms up to the number of redundancies have no effect, as charge can be transferred between any pairs of atoms by an indirect transport pathway, and this will effectively allow distant CT. This is perhaps easiest seen for a cyclic connected system, where an attempt of inhibiting CT between atoms 1 and N by setting p1N = 0 is circumvented by a sequence of p12, p23, …, pN-1N CT that arises when solving the linear equations. It should also be noted that this only depends on the number of non-zero entries in the T matrix, not on their magnitude. Attempts of reducing CT between two atoms by scaling a Tij element by a (small) constant will simply lead to the corresponding pij element increasing by the reciprocal of the constant, and thus result in the same set of atomic charges.


12

Page 13 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Mikulski et al. have generalized the split-charge formulation to limit the amount of transferable charge to a bond-dependent maximum value, and redefining the variable as fji (-1 < fji < +1), where fji includes the asymmetry factor.80 Qi   f ji p max ji

(13)

j

This solves the overpolarization when the overall energy expression is augmented with a tanh-1(fji) factor, but this factor lacks a physical interpretation. They in addition parameterize the maximum transferable bond charge in terms of the bond order, i.e. this becomes a geometry dependent parameter which goes to zero as the bond distance is increased. The transformation between charge and split-charge parameters in eq. (9) led Chen et al. to propose that atomic electronegativitiy parameters should be replaced by effective electronegativities including atomic overlap factors Sij, and this defines the QTPIE model.78 EQTPIE   Qi ieff  12 ij Qi Q j i

(14)

ij

ieff    i   j  j

 

Sij Rij

 

Sij Rij0

(15)

The QTPIE model eliminates the CT problem at dissociation for diatomic systems since the electronegativities go to zero as the distance increases. Mathieu has pointed out, however, that it does not solve the CT problem for multi-atom systems, since the effective electronegativity for a dissociating atom will approach zero, while the remaining atoms will have non-zero electronegativities, and thus leading to non-zero CT.79 The dipole polarization depends only on the hardness matrix, which with fixed hardness parameters predicts an unphysical cubic scaling with system size. The fundamental problem is that charge can migrate with no penalty (metallic behavior) from one end of a system to the other. This


13


overpolarization has in some models been corrected by making the hardness and electrostatic parameters charge dependent,81 which is equivalent to including higher order terms in the Taylor expansion eq. (1), or by an augmented energy expression as in eq. (12).59 In summary, the FQ model displays metallic conductance leading to overpolarization and to longrange CT. Split-charge models are in their native forms just reparameterizations of the FQ model containing additional redundant parameters, except for the special case of non-charged, non-cyclic systems, where there is a one-to-one correspondence between the two sets of parameters. Removing specific split-charge parameters has no effects on the calculated atomic charges until the number of non-zero parameters is less than Natom-1 for a non-fragmented system. Reducing the number of split-charge parameters to less than Natom-1 is equivalent to including additional Lagrange multipliers in eq. (2) for enforcing zero charge on selected fragments. Addition of penalty terms depending on split-charge parameters corresponds to modifying the interaction between atomic charges, which in the long-distance limit should be Coulombic, and addition of penalty terms may thus alter the underlying physics. Bond Capacity Model, Atomic Charges. The FQ model determines the atomic charges by minimizing an energy function (eq. (1)) which provides the electrostatic energy in the force field. A key feature in the present work is to separate the calculation of the atomic charges from the calculation of the electrostatic energy. We will describe the details of how we propose to calculate the atomic charges below, but once they are available, the electrostatic energy is calculated by a screened Coulomb expression, eq. (16). The screening function f can be based on bond connectivity, as in most present force fields, or taken as a distance dependent continuous function based on orbital overlap. The energy expression can be


14

Page 14 of 59

Page 15 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


modified to include charge penetration effects by partitioning the atomic charge Q into core and valence contributions.82-83 Eel 

1 2

 J ij QiQ j  12 

i j

1  fij Rij

i j

Qi Q j

(16)

In the spirit of classical electrostatics and the split-charge models, we propose a bond capacity (BC) model where the charge transferred between an atom pair is proportional to a bond capacity

ij and a difference in electrostatic potential.



pij  ij Vi  V j



(17)

The pair capacityij is parametrized as shown in eq. (18).

ij  Rij   ij0 g  Rij 

(18)

Here 0𝑖𝑗 is the bond capacity at the equilibrium geometry for a given pair of atom types, and it will be closely related to the bond order, as shown in the Results and Discussion section. The attenuation function g(Rij) takes a value of one at equilibrium distance and decays to zero for large distances which ensures zero CT between non-interacting atom pairs. The attenuation function can be considered as a normalized orbital overlap and will in the present work be parameterized in terms of overlaps between Slater type orbitals. The electrostatic potential Vi is written as a sum of three components as shown in eq. (19). Vi  i  i   i

(19)

Here i represents the intrinsic local potential, i is the potential due to the remaining (N-1) atomic charges while i is the potential due to an external field. The i is the main parameter of the BC


15


Page 16 of 59

model and is identified as the atomic electronegativity, while i provides intra- and inter-molecular charge polarization and i is included to derive an expression for the total molecular polarizability. The net atomic charge is parameterized by eq. (20), where the sign is chosen such that the more electronegative atom will carry a negative charge. N



Qi    ij Vi  V j j 1



(20)

Since the pair capacities are elements of the symmetric matrix , the antisymmetry given by the potential difference automatically guarantees charge neutrality of the system (see ref. 84 for similar ideas). Eq. (20) expresses the atomic charge as a sum of contributions from atom pairs within atomic wave function overlap (eq. (18)) and proportional to the difference in their electrostatic potentials. The double indexed potential differences in eq. (20) can be avoided by rewriting it in a matrix-vector notation as shown in eq. (21). Q  CV

(21)

The symmetric capacitance matrix C is constructed from bond capacities according to eq. (22).  N   ik Cij   k  j   ij

i j i j

(22)

An example for a three-atom system is given in eq. (23). 12 13   12  13   C    21  21   23  23    32 31  32  31 

(23)

Eq. (22) is a special form of the Maxwell capacitance matrix where the self-capacitance elements

ii are set to zero.85 The form of eq. (22) guarantees that the sum of all elements in a given row/column of C is zero. We note that the relation between charges and potentials in eq. (21) was


16

Page 17 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


also derived by Chen and Martinez in the circuit representation of QTPIE model for the special case of a diatomic system where the charge conservation is given implicitly by the Q1 = -Q2 condition.61-62 The molecular polarization potential i is given by a Coulomb interaction matrix J which includes a screening factor for atoms within bonding distance while the diagonal elements Jii describe selfinteraction potentials and are equivalent to the atomic hardness ii in eq. (1). The screening factor f (eq. (16)) will in the present work be taken to be the same as the attenuation function g in eq. (18) but could be chosen to have a different functional form. A corresponding term is used in the ReaxFF, where the screening factor is based on a distance cubed to the 1/3 power.86 N

N

1  fik Qk  J ii Qi k i Rik

i   J ik Qk   k 1

(24)

Combining eqs. (19), (21) and (24) leads to eq. (25), where I denotes the identity matrix.

 I  CJ  Q  C  χ  ψ 

(25)

Assuming that the (I + CJ) matrix is non-singular, the set of atomic charges can be obtained by matrix inversion, eq. (26). Q    I  CJ 

1

Cχ  ψ

(26)

Eq. (26) automatically enforces charge neutrality, which follows from eq. (20) and it is ensured by the properties of C, as shown explicitly in Appendix A. Eq. (26) represents a direct solution by matrix inversion, but eq. (25) can alternatively be solved by iterative methods for large systems. Charged systems can be handled by augmenting eq. (21) by an additive Q0 term as shown in eq. (27). Q  CV  Q0


17

(27)


Page 18 of 59

The N-dimensional vector Q0 contains the reference charges assigned to the atoms and the sum of the Q0 elements is the total molecular charge. We will in the present case take the Q0 elements as the formal charges assigned by a Lewis approach (e.g. -0.5 for the two oxygen atoms in a deprotonated carboxylic acid) but other choices could be made. The V potential in eq. (27) now contains a molecular polarization potential  arising from both Q and Q0 as shown in eq. (28).



  J Q  Q0



(28)

The combination of eqs. (19), (27) and (28) can after some rearrangements be written as in eq. (29), which is the analog of eq. (26) for non-neutral systems. Q    I  CJ 

1

C  χ  ψ    I  CJ 

1

 I  CJ  Q0

(29)

Eq. (29) describes a set of polarizable charges where the net molecular charge is given by eq. (30), as proven in Appendix A. Qtot  1t Q0

(30)

Appendix A shows that eq. (29) leads to a molecular dipole moment with the correct behavior for a change of origin for charged systems. Bond Capacity Model, Molecular Polarizability. The atomic charges in eq. (29) lead to the molecular dipole moment shown in eq. (31), where R is a (N,3) array containing the atomic positions. μ  R t Q   R t  I  CJ 

1

C  χ  ψ   R t  I  CJ 

1

 I  CJ  Q0

(31)

The external potential  can be used to identify the induced dipole component as shown in eq. (32). μind   R t  I  CJ 

1

Cψ


18

(32)

Page 19 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Eq. (32) shows that the induced dipole moment does not depend on the net charge of the molecular system. We now take the external potential  to be a linear electric field as shown in eq. (33).

i  



  x, y , z

Ri F

(33)

This allows the induced dipole moment to be written as in eq. (34) which by differentiation with respect the field F leads to the dipole polarizability tensor in eq. (35). μind  R t  I  CJ  α  R t  I  CJ 

1

1

CRF

CR

(34) (35)

Eq. (35) shows that the polarizability tensor depends on the capacitance and Coulomb matrices but it is independent on the net charge of the molecular system. The polarizability of neutral and cationic/anionic systems will in reality be significantly different, and this implies that bond capacity parameters must depend on the charge state of the molecule. The polarizability in eq. (35) has the correct translational and rotational symmetries and scales linearly with the system size, as shown in Appendix B. For a diatomic system at equilibrium distance and assuming no self-interaction potentials, the J, and thus CJ, matrices are zero due to the fact that the charges are completely screened (eq. (24)) and the C matrix is given by a single bond capacity , as shown in eq. (36).  C  

   

(36)

The diagonal and off-diagonal polarizability components are thus given by eqs. (37) and (38), where  and  denote Cartesian components.

    A   B 

2

    A   B  A   B 


19

(37) (38)


This agrees with the formulae derived by Holt et al. and provides insight into the bond contribution to the total molecular polarizability.87 Eq. (37) shows that the  polarizability component depends linearly on the bond capacity  and quadratically on the th component of the bond distance. This accounts for the increase in polarizabilities for a set of molecules such as HF, HCl, HBr, where the bond order and thus bond capacity is expected to be similar, but the equilibrium bond distance increases and thereby leads to an increase in the polarizability. Computational Details. All electronic structures calculations were performed with the Gaussian 09 program package88 using the B97X or B97X-D functionals89 and the aug-pcseg-390 and cc-pVDZ91 basis sets for polarizabilities and electron densities used to obtain atomic partial charges, respectively. The calculations of Iterative Hirshfeld charges were carried out with the Multiwfn software.92 The bond capacity calculations and parametrizations were performed with home-developed codes. Results and Discussion. The main parameters in the BC model are the bond capacities  and atomic electronegativities  and these allow calculation of the atomic charges according to eq. (29). It is clear, however, that there will be a strong coupling between these parameters, as the combination of large electronegativity differences and small bond capacities may lead to similar atomic charges as the combination of small electronegativity differences and large bond capacities (eq. (20)). It is therefore of interest to probe how these parameters can be assigned independently, and we propose that bond capacities parameters can be assigned to reproduce molecular polarizabilities, while electronegativity parameters are assigned subsequently to fit a set of reference atomic charges that models the molecular ESP. We note that the set of bond capacity and atomic electronegativity


20

Page 20 of 59

Page 21 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


parameters will depend on the choice of the Coulomb self-interaction (atomic hardness) and screening function f in eq. (24), but this dependence will be relatively small. We will in the present case keep the Coulomb screening function constant and show results corresponding to Coulomb self-interaction values of 0 and 1. The latter is chosen only to probe the behavior of the results towards a non-zero value, and the Coulomb self-interaction could in a large scale parameterization be chosen as atom-type specific parameters to allow further tuning on the BC model. Molecular polarizabilities using molecular optimized parameters Eq. (35) shows that molecular polarizabilities can be calculated from the molecular geometry and a set of bond capacities parameters. We will in this section probe the inherent limitations of the model by optimizing non-symmetry related bond capacity parameters to reproduce reference polarizabilities for each system, while the next section will probe how much accuracy is lost by assuming a common set of bond capacities based on atom types. We will take the self-interaction potential parameters Jii in eq. (24) to be zero, other choices will lead to the same final polarizabilities but influence the values of the bond capacity parameters. Table 1 shows the results for saturated hydrocarbons with reference values taken from B97XD/aug-pcseg-3 calculations, and where the column labelled “Bond types” shows the number of fitted parameters. The BC model clearly is capable of accurately reproducing the reference results when the bond capacity parameters are optimized for each molecular system. Butane shows the largest deviation of the considered hydrocarbons, and this is a consequence of using only one bond capacity for all C-H bonds. A better agreement can be obtained by differentiating between C-H bonds in methyl and methylene groups.


21


Modelling the ESP using only nuclear centered atomic charges is unable to describe out-of-plane polarization for planar molecules and this has been one of the motivations for developing the atomic dipole polarizability model. The static ESP can be improved by either including atomic dipole and quadrupole moments, or by including off-nuclei charges, as for example in the TIP4P and TIP5P water models.93 Including off-nuclei sites in the BC model allows modelling out-ofplane polarization for planar molecules and significantly improves polarization components along molecular direction involving lone pairs. Table 2 shows the results for water and ammonia with only nuclear centered charges and with a model including two or one off-nuclei dummy atoms as illustrated in Figure 1. The zz-polarizability components are poor without off-nuclei sites, but including dummy atoms in the classical lone-pair directions allows reproduction of the reference values. The results in Table 2 for methanol, methylamine, dimethylether and oxirane show that the requirement of lone-pair dummy atoms may be specific for small systems. When the alcohol/amine/ether functional group is coupled to methyl/methylene groups, there is sufficient flexibility in the parameter space to fit the molecular polarizability without including dummy atoms. This reflects that the number of reference data (three) is smaller than the number of free parameters, and the methyl/methylene groups are capable of accommodating the lack of polarizability in the sp3-like oxygen/nitrogen moiety. The use of off-nuclei dummy atoms is necessary for modelling of out-of-plane polarizabilities for planar system, with off-nuclei sites defined as shown in Figure 1. Table 3 shows the results for a selection of planar and unsaturated molecular systems, and it is significant that the BC model can accurately reproduce the anisotropic polarizabilities for these systems. Molecular polarizabilities using common parameters


22

Page 22 of 59

Page 23 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The results in Tables 1-3 show that molecular polarizabilities can be reproduced within a few percent by assigning molecular specific bond capacities for each non-symmetry related bond. As illustrated in Table 2, however, this raises the question of whether the obtained values are just a result of data fitting, rather than the model capturing the underlying physics. We have thus examined the optimized values of the bond capacity parameters for the systems in Tables 1-3 and extracted a set of common (averaged) values shown in Table 4, where parameters involving dummy atoms (X) have only been defined relative to sp2-atoms, as sp3 lone-pair dummy atoms appear unnecessary, except for water and ammonia (Table 2). Only bond types present in at least three molecules have been taken into account for averaging, and a complete list of optimized bond capacity parameters is given in Table S1 in the supporting information. As discussed in the Bond Capacity Model, Atomic Charges section, the BC model assumes that the bond capacities are related to bond orders, and this is nicely reflected by the values in Table 4. The bond capacities involving dummy atoms are significantly larger than those involving real atoms, but this is largely due to the choice of a C-X distance of 1 Bohr. According to Eq. (37), the principal polarizability components scale quadratically with the bond length while they are linear in the bond capacity parameter. This implies that in order to get the same polarizability components for a pair of bonds with lengths d and d/2, the bond capacity of the latter must be four times the former. The bond capacity parameters extracted from molecular polarizabilities are in qualitative agreement with ab-initio calculated bond capacities (Tables S1 and S2 in supporting information),94 and this supports that the bond capacity parameters carry a physical insight. Table 5 shows the calculated polarizabilities using the bond capacity parameters in Table 4, and displays, as expected, larger deviations than with molecular optimized parameters. The mean absolute error for the set of molecules in Table 5 relative to the reference values is 10.8% while


23


the corresponding value for the molecule specific bond capacities polarizability in Tables 1-3 is 0.4%. The common set of parameters in Table 4 reflects a typical classification into atom types within a force field model. We note, however, that the number of different atom types is a choice, which is tightly coupled to the compromise between accuracy and generality, where a high accuracy normally will require a finer graduation between atom types. Table 5 thus shows that the BC model with a set of physically motivated parameters is capable of reproducing anisotropic molecular polarizabilities with an accuracy of ~10%, but there is clearly ample room for refinements. Given a predefined target accuracy and a large set of reference molecules, it should be possible to deduce the required graduation in terms of atom types and derive suitable average bond capacity parameters for reproducing molecular polarizabilities. Intermolecular polarization. The previous section shows that bond capacity parameters can be chosen to reproduce reference molecular polarizabilities, and this section will show how atomic charges can be modelled by proper choices of atomic electronegativities, when the bond capacities are fixed at the values in Table 4. The molecular polarizability and ESP are in principle experimentally observable quantities, but atomic charges are not, and the choice of atomic charges for fitting the electronegativity parameters is therefore somewhat arbitrary. Atomic charges obtained by fitting to the molecular ESP are employed by many force fields, but the fitting procedure is often illconditioned and may display unphysical large changes in atomic charges with small geometry perturbations,50, 95 and this has led to restrained fitting procedures.96 Atomic charges obtained by the iterative Hirshfeld (HI) procedure have been shown to reproduce ab-initio computed ESP almost as well as ESP-fitted atomic charges.27 We will in the present case use HI charges as the


24

Page 24 of 59

Page 25 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


reference, but emphasize that the BC model can be parameterized to reproduce other choices as well, or directly parameterized to reproduce the reference ESP. Intermolecular polarization is in a charge-only force field modelled by allowing the atomic charges to depend on the intermolecular geometry. The difference relative to the situation in the previous section is that the static linear field perturbation is replaced by an anisotropic perturbation arising from the combined Coulomb potential of all atomic charges, and that the mutual polarization between molecules is also taken into account. We will as illustration take the system corresponding to the fluoromethane dimer shown in Figure 2, where the atomic electronegativity parameters of the monomer have been chosen to reproduce the HI charges calculated at the B97X/cc-pVDZ level. Figure 3 shows the change in the dipole moment of the A-monomer as a function of the interfragment distance between 6 and 13 Å, where the latter value has been chosen for aligning the result corresponding to choices of 0 and 1 for the Jii self-interaction parameters in eq. (24). Both Jii values lead to a qualitatively correct behavior and nicely bracket the reference results obtained using the HI charges. The polarization due to changes in the atomic charges obtained by setting all Jii = 0 is larger than with all Jii = 1, and this follows straight-forward from eq. (26). The atomic self-interaction parameters thus allow a tuning of the magnitude of the polarization response. As mentioned in the Results and Discussion section, we will not attempt to select optimum values for the self-interaction parameters but restrict ourselves to show the behavior for values of 0 and 1. For the CH3F molecule, there would be three atomic self-interaction parameters and clearly a common value of 0.5 would be a sensible choice for this system. We also refrain from probing polarization effects at distances shorter than 6 Å, as wave function overlap becomes important at short distances, with results depending on the exact choice of the Coulomb screening


25


function f in eq. (24) as well as allowing intermolecular CT effects from the g attenuation function in eq. (18). Figure 5 shows the change in the monomer A x- and z-dipole components when the position of monomer B is generated by a translation in the xz-plane with a constant radius vector r of 8 Å, as illustrated in Figure 4. The change in the z-component in Figure 5 mirrors the behavior in Figure 3 with the reference HI values between the Jii = 0 and Jii = 1 models. The change in the xcomponent in Figure 5 is qualitatively correct, but the results with the Jii = 0 model is here closest to the reference HI results. The Jii = 1 model underestimates the polarizability in the x-direction and this reflects the underestimation of the xx component shown in Table 5. A different choice of bond capacity parameters to better model the xx polarizability would increase the response of the x-dipole component, and this could also be obtained by including dummy atoms for modelling the lone-pairs on fluorine. We will in the present case refrain from such parameter tuning, but simply note that the BC model appears to capture the underlying physics and that model refinements could increase the agreement with reference data. The polarization of CH3F by CH3NH3+ (supporting information) is very similar to Figures 3 and 5, and shows that the BC model can handle neutral and charged systems within the same general framework. Intramolecular polarization. The intermolecular polarization discussed in the previous section should transfer to intramolecular polarization when the two moieties are sufficiently far removed in terms of bonding. The evidence from spectroscopic techniques such as NMR and IR suggests that through-bond effects become small beyond ~three bonds, with slightly longer decay distance for conjugated systems. In order to probe intramolecular polarization we have chosen 1-fluoro-decane where the fluoro-methylene group is rotated 180° to induce a charge polarization of the alkyl chain, as illustrated in Figure 6.


26

Page 26 of 59

Page 27 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Electronegativity parameters have again been optimized to reproduce the HI charges for the 0° geometry using the bond capacities in Table 4. Atomic charges are, as already mentioned, not observable quantities and only serve as a vehicle for modelling the molecular ESP. Since different sets of atomic charges may lead to very similar ESP, we will for the intramolecular case focus on the latter. We will use (R0 , Q0) and (R180, Q180) as notation for the geometry and charges before and after the rotation of the fluoro-methylene group and Q as the difference between the Q180 and Q0 charges. The through-space polarization effect is quantified by plotting the ESP in the xy-plane arising from Q for the C5-C10 atoms with the corresponding hydrogens, i.e. atoms that are at least five bonds removed from the perturbation. Figure 7 shows the reference Q potential obtained using HI charges and the Q potentials from the BC model using self-interaction parameters Jii of 0 and 1. The BC model with Jii = 0 leads to an overpolarization by roughly an order of magnitude and fails to qualitative reproduce the reference result. The Jii = 1 BC model, however, correctly reproduces both the qualitative behavior and the magnitude of the polarization. We again note that fine-tuning of atomic self-interaction parameters could further improve the agreement, but a much larger set of reference data will be required for probing the limitations. The corresponding results by rotating the methyleneammonium group in 1-ammonium-nonane are very similar and shown in Supplementary Material, and illustrate that neutral and charged systems can be handled in the BC model within the same framework. The results in this section show that through-space intra-molecular polarization can be treated analogous to intermolecular polarization by the BC model, in contrast to the FQ model that display metallic conductance and leads to overpolarization. The conclusion is that a non-zero value for the Coulomb self-interaction is required to avoid overpolarization in the BC model, but the details of


27


the numerical value, and whether a single or a number of atom specific parameters is required, will have to await testing on a larger selection of systems. We have not discussed how the BC model can incorporate through-bond intramolecular polarization, defined as changes in atomic charges within ~three bonds of the geometry perturbation. Modelling these effects will be sensitive to the exact nature of the attenuation function in eq. (18) as well as the Coulomb screening function in eq. (16)/(24), and this will be the subject of a separate study. Coulomb interactions between atoms separated by one or two bonds are ignored in most common force fields, and downscaled for atom pairs separated by three bonds in some force fields, and this essentially eliminates the need for modelling through-bond intramolecular polarization, but the BC model has potential for including such effects. We note that it recently has been shown that short-range polarization effects can be reproduced by a machine-learning approach.97-98 Conclusion. We propose a bond capacity model that is capable of modelling both inter- and intra-molecular polarization using only atomic charges. The model contains two main sets of parameters, bond capacities and atomic electronegativities. The bond capacity parameters control the flow of charge between atoms, with the atomic electronegativities providing the inherent potential for transferring charge. The bond capacity is shown to be closely associated with the bond order, and can be parameterized against molecular polarizabilities, with the derived parameters showing a good degree of transferability. Atomic electronegativities can subsequently be parameterized to reproduce molecular electrostatic potentials, either directly or indirectly by fitting to atomic charges derived from a reference wave function. The two sets of parameters can thus be derived from experimentally observable quantities, molecular polarizabilities and electrostatic potentials, rather than from results using non-unique decompositions into atomic quantities. The bond


28

Page 28 of 59

Page 29 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


capacity model is free from artifacts displayed by the fluctuating charge model, such as charge transfer over infinite distance and non-linear scaling of the polarizability with system size, and treats neutral and charged systems on the same footing. The bond capacity model describes charge polarization and is a natural lower order expansion than the polarizable dipole model, and it is of interest to establish its limitations before possibly including higher order terms. While a model with charges only on nuclear positions is unable to describe out-of-plane polarization for planar systems, we show that inclusion of off-nuclei charge sites allows an accurate representation of all the major polarization effects. We believe that the bond capacity model could be an interesting component of designing new force field models.


29


FIGURE CAPTIONS Figure 1. Illustrating the position of dummy atoms for water, ammonia and planar molecules. The distance to the dummy atom is in all cases taken as 1 Bohr. The molecular plane for water is xz, with dummy atoms placed in the zy-plane with an XOX angle of 97.8 while for ammonia the dummy atom is placed along the C3 axis. Structures A, B, C show the placement of dummy atoms for planar molecules/groups as benzene and ethane, carbonyl groups and boranes, respectively. Figure 2. The fluoromethane dimer for probing the intramolecular polarization shown in Figure 3. Figure 3. Change in the dipole moment for monomer A (Figure 2) as a function of intermolecular distance with results aligned at 13 Å Figure 4. Illustrative CH3F dimer geometries used to probe intermolecular polarization shown in Figure 5 as a function of angle. Figure 5. Change in the x- and z-dipole components with zero taken at  = 0 for monomer A with the geometry variation shown in Figure 4. Figure 6. 1-Fluoro-decane conformations used for probing intramolecular polarization shown in Figure 7 with the carbon backbone laying in the xy-plane. Figure 7. Electrostatic potential in atomic units in the xy-plane arising from the atomic charges difference for the C5-C10 atoms and bonded hydrogens for 1-fluoro-decane (Figure 6) for BC charges with Jii = 0, Jii = 1 and HI. Figure A1. Principal polarizability component calculated with the BC model for a system of up to 30 CO molecules aligned along the z-axis, with an interfragment distance of 3 Å and a bond capacity of 3.273 a.u.


30

Page 30 of 59

Page 31 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


TABLES Table 1. Polarizability components (atomic units) calculated by the BC model compared to reference values, where bond types indicate the number of fitted parameters. The mean absolute deviation for all the components is 0.59. Molecule BC polarizabilities B97X-D/aug-pcseg-3 Bond types Methane Ethane Propane Cyclopropane Butane Cyclobutane Fluoro-methane Fluoro-ethane

XX

YY

ZZ

XX

YY

ZZ

17.2 27.8 42.3 38.0 56.2 50.4 17.0 32.2

17.2 27.8 37.8 38.0 45.9 50.4 17.0 31.2

17.2 30.8 41.7 34.6 53.4 44.2 17.8 34.0

17.2 27.8 41.4 37.9 51.8 50.4 17.0 32.3

17.2 27.8 38.1 37.9 48.1 50.4 17.0 31.2

17.2 31.0 42.8 34.5 58.1 44.2 17.8 34.0

1 2 2 2 2 2 2 3

Table 2. Polarizability components (atomic units) calculated by the BC model compared to reference values, where bond types indicate the number of fitted parameters. With and without offnuclei refers to the dummy atoms shown in Figure 1. The mean absolute deviation for all the components is 1.93 without off-nuclei and 0.01 with off-nuclei sites. Molecule

Without off-nuclei

Water Ammonia Methanol Methylamine Dimethylether Oxirane

B97X-D/augpcseg-3

BC polarizabilities With off-nuclei

XX

YY

ZZ

XX

YY

ZZ

XX

YY

ZZ

11.8 16.6 20.8 26.5 40.3 29.0

0.0 16.6 21.8 25.8 34.7 29.6

3.0 1.4 22.0 30.0 44.2 26.7

10.9 16.3 20.8 26.6 40.3 29.0

9.8 16.3 21.8 25.8 34.7 29.4

10.3 17.0 21.9 30.0 44.4 26.6

10.9 16.3 20.8 26.6 40.3 29.0

9.8 16.3 21.8 25.8 34.7 29.5

10.3 17.0 21.9 30.0 44.4 26.7


31

Bond types With/without off-nuclei 1/2 1/2 3/4 3/4 3/4 3/4


Page 32 of 59

Table 3. Polarizability components (atomic units) calculated by the BC model compared to reference values, where bond types indicate the number of fitted parameters. The BC results are with dummy atoms as shown in Figure 1. The mean absolute deviation for all the components is 0.03. Molecule BC polarizabilities B97X-D/aug-pcseg-3 Bond types Ethene Benzene 1,1-Difluoroethene E-1,2-Difluoroethene Fluoro-benzene 1,4-Difluorobenzene Formaldehyde Acetaldehyde Acetone Borane Trifluoroborane

XX

YY

ZZ

26.1 73.8 28.4 27.9 40.2 39.8 18.0 31.3 45.6 19.4 16.9

22.3 73.8 22.0 21.4 73.8 72.8 13.1 24.2 34.4 14.9 13.6

34.8 40.6 35.6 35.1 75.3 76.9 22.0 35.0 45.8 19.4 16.9

XX

26.1 73.8 28.4 27.9 40.2 39.8 18.0 31.3 45.6 19.4 16.9

YY

22.3 73.8 22.0 21.4 73.3 73.1 13.1 24.2 34.4 14.9 13.6

ZZ

34.8 40.6 35.6 35.1 75.3 76.6 22.0 35.0 45.8 19.4 16.9

3 3 4 4 4 4 3 5 4 2 2

Tabel 4. Common set of bond capacity parameters with Car and Ccar denoting aromatic and carbonyl carbon, respectively, and X denoting a dummy atom. Bond

C-H

Car-H

C-C

C=C

Car=Car

C=O

C-F

C-O

C-X

Car-X

Ccar-X

Length(bohr)

2.08

2.08

2.83

2.55

2.46

2.27

2.55

2.36

1.0

1.0

1.0

0

1.317

0.526

1.072

1.920

1.222

1.891

1.316

1.557

5.480

3.351

5.924


32

Page 33 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Tabel 5. Polarizability components (atomic units) calculated by the BC model using the parameters in Table 4 compared to reference values. The BC results for planar systems are with dummy atoms as shown in Figure 1. The mean absolute deviation for all the components is 4.1. Molecule

BC polarizabilities B97X-D/aug-pcseg-3

xx

yy

zz

xx

yy

zz

Methane

14.7

14.7

14.7

17.2

17.2

17.2

Ethane

28.5

28.5

30.8

27.8

27.8

31.0

Propane Cyclopropane

48.5 41.5

41.3 41.5

46.8 41.9

41.4 37.9

38.1 37.9

42.8 34.5

Butane

63.0

53.6

56.2

51.8

48.1

58.1

Cyclobutane Fluoromethane Fluoroethane Ethene Benzene 1,1-Difluoroethene

67.2 13.9 24.5 28.0 77.8 30.5

67.2 13.9 27.9 21.9 77.8 21.9

51.0 17.7 27.2 39.3 40.2 34.3

50.4 17.0 32.2 26.1 73.9 28.4

50.4 17.0 31.2 22.3 73.8 22.0

44.2 17.8 34.0 34.8 40.6 35.6

E-1,2-Difluoroethene 30.7

21.9

34.6

27.9

21.4

35.2

Fluoro-benzene 1,4-Difluorobenzene Formaldehyde Acetaldehyde

40.2 40.2 14.2 26.0

73.5 69.4 11.9 26.6

75.1 72.3 23.0 33.0

40.2 39.8 18.0 31.3

73.3 73.1 13.1 24.2

75.3 76.6 22.0 35.0

Acetone Dimethylethere Oxirane

37.6 17.1 32.8

40.2 30.3 32.6

42.7 23.6 28.2

45.6 40.3 29.0

34.4 34.7 29.5

45.8 44.4 26.7


33


Supporting Information. Figures corresponding to Figures 2-7 for a similar system with –NH3+ instead of –F and Tables with bond capacities are provided as Supporting Information. The files are available free of charge.

Corresponding Author *[email protected] Funding Sources This work was supported by Grant No. 4181-00030B from the Danish Council for Independent Research.


34

Page 34 of 59

Page 35 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60



35


Page 36 of 59

Appendix A. Eq. (29) ensures charge conservation, which can be written as in eq. (A1). Qtot  1t  Q neut  Qion 

1

 1t Q0

(A1)

The neutral and ionic parts of the polarized charges are given by eqs. (A2) and (A3). Qneut    I  CJ  Qion   I  CJ 

1

1

Cχ  ψ

(A2)

 I  CJ  Q0

(A3)

The (I+CJ)-1 matrix can be written in its Taylor-expanded form shown in eq. (A4).

 I  CJ 1  I  CJ  CJCJ  CJCJCJ  L

(A4)

This allows eq. (A2) to be written as eq. (A5). 1t Q neut  1t  I  CJ  CJCJ  CJCJCJ  L  C  χ  ψ 

(A5)

Equation (A5) can be rearranged to obtain eq. (A6), where all terms vanish due to the property of the C matrix (eq. (22): C1 = 1tC = 0). 1t Q neut  1t C  I  JC  JCJC  JCJCJC  L

 χ  ψ   0

(A6)

Eq. (A3) can similarly be rewritten and rearranged to give eq. (A7), where the property of the C matrix ensures that only the first term survives. 1t Qion  1t Q0  1t C  2JC  2JCJC  2JCJCJC  L  Q0  1t Q0  Qtot

(A7)

The dipole moment of a system with a net charge Qtot is not invariant to a translation. The change by a translation in the  direction (𝑇) of magnitude t is given in eq. (A8), and this we want to prove for a dipole moment generated by the atomic charges in eq. (29). Tµ     t Qtot

(A8)

The dipole component  can be written using the atomic charges in eq. (29) as shown in eq. (A9).


36

Page 37 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


  Rt  Q neut  Qion 

(A9)

The change of origin described by 𝑇 can be written as in eq. (A10). t Tµ    R  t 1  Q neut  Qion 

(A10)

Rearranging eq. (A10) leads to eq. (A11), which can be simplified by using eqs. (A6) and (A7), and this proves eq. (A8). Tµ   Rt  Q neut  Qion   t 1t Q neut  t 1t Qion    t Qtot

(A11)

Appendix B We here show that the polarizability given in eq. (35) scales linearly with system size. We assume a super-system of N identical, non-interacting and uncharged fragments, with a constant translation vector between neighboring fragments. The choice of neutral fragments does not affect the generality of this derivation, as the polarizability tensor does not depend on the net molecular charge, as shown in eq. (35). Eq. (25) can for this super-system be written as in eq. (B1).

 I  CJ  N  Q N   C  χ  ψ  N 

(B1)

An out-of-diagonal element of (I+CJ) can be decomposed by the following summations:

 I  CJ ik N    ik   Cij J jk   Cij J jk jI

(B2)

jI

In eq. (B2), I represents the fragment the atom i belongs to. The last summation in eq. (B2) runs over null out-of-diagonal elements of C since the pair capacity involving two atoms belonging to different fragments is zero (Cij = -ij = 0) by the non-interaction assumption. The first summation in eq. (B2) is also zero if atom k does not belong to I as Jjk is zero for non-interacting fragments. Eq. (B2) can thus rewritten as in eq. (B3).


37


N

 I  CJ ik

  ik   Cij J jk  j   ik 

Page 38 of 59

 j, k   I

(B3)

 j, k   I

Eq. (B3) implies that (I+CJ) is N block-diagonal and its inverse (assuming that it exists) can therefore be written as in eq. (B4).

 I  CJ 1 N    I  CJ 11  L

  I  CJ 

1 I 

 L   I  CJ 

1 N 

(B4)

A similar analysis can be performed on the right-hand side of eq. (B1). Cχ  ψ

N







  Cij  j   j   Cij  j   j jI

jI



(B5)

The last summation is zero by the same arguments as for eq. (B2), and eq. (B5) can be written as direct sum of fragment contribution as shown in eq. (B6). N 1 I  N Cχ  ψ  C  χ  ψ  L  C  χ  ψ  L  C  χ  ψ 

(B6)

Combining eqs. (B4-B6) leads to eq. (B7). N 1 I N Q   Q   L  Q   L  Q 

(B7)

The fragments are related to one another by a translation that preserves atomic distances within a fragment, and the absence of rotations implies that each fragment is polarized in the same direction by an external field. Since all fragments are identical, all blocks in eqs. (B4), (B6) and (B7) are therefore also identical. The  dipole moment component of the super-system is given by a sum over single fragments contributions.

 

t

 N   Rt Q N    R I  Q I      I  I

(B8)

I

The dipole moment of a neutral molecule is invariant under translations, as shown by eq. (A8). Since the fragments are assumed to be uncharged, identical and related to one another by


38

Page 39 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


translations, the  dipole component of the super-system can be written as in eq. (B9), where the

 indicates the  dipole moment of a single fragment.  N   N 

(B9)

Combining eqs. (34) and (B9) allows the induced dipole moment of the super-system to be written as in eq. (B10).

 N ,ind  NRt  I  CJ  CRF 1

(B10)

Differentiating eq. (B10) with respect the  field component gives the  polarizability element shown in eq. (B11).  N   NR t I  CJ   CR  N   1

(B11)

Eq. (B11) shows that the bond capacity model displays a linear scaling of the polarizability tensor for a system of non-interacting fragments. This has been numerically tested for super-systems of up to 30 CO monomers aligned along the z-axis and separated by 3 Å, as shown in Figure A1.


39


References 1.

Frenkel, D.; Smith, B., Understanding Molecular Simulations. Academic Press:

San Diego, USA, 2002. 2.

Allen, M. P.; Tildesley, D. J., Computer Simulations of Liquids. Clarendon Press:

Oxford, UK, 1987. 3.

Schlick, T., Molecular Modeling and Simulation. Springer: 2002.

4.

Vanommeslaeghe, K.; Hatcher, E.; Acharya, C.; Kundu, S.; Zhong, S.; Shim, J.;

Darian, E.; Guvench, O.; Lopes, P.; Vorobyov, I.; MacKerell, A. D., Jr., CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields. Journal of Computational

Chemistry 2010, 31, 671-690. 5.

Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G. M.; Zhang, W.; Yang, R.;

Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J. M.; Kollman, P., A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations. Journal of Computational Chemistry 2003, 24, 19992012. 6.

Plazinski, W.; Lonardi, A.; Huenenberger, P. H., Revision of the GROMOS

56A6(CARBO) force field: Improving the description of ring-conformational equilibria in hexopyranose-based carbohydrates chains. Journal of Computational Chemistry 2016,

37, 354-365. 7.

Cieplak, P.; Caldwell, J.; Kollman, P., Molecular mechanical models for organic

and biological systems going beyond the atom centered two body additive approximation: Aqueous solution free energies of methanol and N-methyl acetamide, nucleic acid base, and amide hydrogen bonding and chloroform/water partition coefficients of the nucleic acid bases. Journal of Computational Chemistry 2001, 22, 1048-1057. 8.

Lopes, P. E. M.; Huang, J.; Shim, J.; Luo, Y.; Li, H.; Roux, B.; MacKerell, A. D.,

Polarizable Force Field for Peptides and Proteins Based on the Classical Drude Oscillator. Journal of Chemical Theory and Computation 2013, 9, 5430-5449.


40

Page 40 of 59

Page 41 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


9.

Patel, S.; Brooks, C. L., CHARMM fluctuating charge force field for proteins: I

parameterization and application to bulk organic liquid simulations. Journal of

Computational Chemistry 2004, 25, 1-15. 10.

Shi, Y.; Xia, Z.; Zhang, J.; Best, R.; Wu, C.; Ponder, J. W.; Ren, P., Polarizable

Atomic Multipole-Based AMOEBA Force Field for Proteins. Journal of Chemical Theory

and Computation 2013, 9, 4046-4063. 11.

Stone, A. J., Distributed multipole analysis: Stability for large basis sets. Journal

of Chemical Theory and Computation 2005, 1, 1128-1132. 12.

Gagliardi, L.; Lindh, R.; Karlstrom, G., Local properties of quantum chemical

systems: The LoProp approach. Journal of Chemical Physics 2004, 121, 4494-4500. 13.

Williams, D. E., Rev. Comp. Chem. 1991, 2, 219.

14.

Breneman, C. M.; Wiberg, K. B., Determining atom-centered monopoles from

molecular electrostatic potentials - the need for high sampling density in formamide conformational-analysisS. Journal of Computational Chemistry 1990, 11, 361-373. 15.

Hirshfeld, F. L., Bonded-atom fragments for describing molecular charge-

densities. Theoretica Chimica Acta 1977, 44, 129-138. 16.

Pérez de la Luz, A.; Aguilar-Pineda, J. A.; Méndez-Bermúdez, J. G.; Alejandre,

J., Force Field Parametrization from the Hirshfeld Molecular Electronic Density. Journal


Verstraelen, T.; Pauwels, E.; De Proft, F.; Van Speybroeck, V.; Geerlings, P.;

Waroquier, M., Assessment of Atomic Charge Models for Gas-Phase Computations on Polypeptides. Journal of Chemical Theory and Computation 2012, 8, 661-676. 18.

Verstraelen, T.; Speybroeck, V. V.; Waroquier, M., The electronegativity

equalization method and the split charge equilibration applied to organic systems: Parametrization, validation, and comparison. The Journal of Chemical Physics 2009,

131, 044127. 19.

Vilseck, J. Z.; Tirado-Rives, J.; Jorgensen, W. L., Evaluation of CM5 Charges for

Condensed-Phase Modeling. Journal of Chemical Theory and Computation 2014, 10, 2802-2812.


41


20.

Sedghamiz, E.; Nagy, B.; Jensen, F., Probing the Importance of Charge Flux in

Force Field Modeling. Journal of Chemical Theory and Computation 2017, 13, 37153721. 21.

Silva, A. F.; Richter, W. E.; Meneses, H. G. C.; Faria, S. H. D. M.; Bruns, R. E.,

How Accessible Is Atomic Charge Information from Infrared Intensities? A QTAIM/CCFDF Interpretation. The Journal of Physical Chemistry A 2012, 116, 82388249. 22.

Silva, A. F.; Richter, W. E.; Bassi, A.; Bruns, R. E., Dynamic atomic contributions

to infrared intensities of fundamental bands. Physical Chemistry Chemical Physics 2015, 17, 30378-30388. 23.

Soderhjelm, P.; Krogh, J. W.; Karlstrom, G.; Ryde, U.; Lindh, R., Accuracy of

distributed multipoles and polarizabilities: Comparison between the LoProp and MpProp models. Journal of Computational Chemistry 2007, 28, 1083-1090. 24.

Jakobsen, S.; Jensen, F., Systematic Improvement of Potential-Derived Atomic

Multipoles and Redundancy of the Electrostatic Parameter Space. Journal of Chemical

Theory and Computation 2014, 10, 5493-5504. 25.

Jakobsen, S.; Jensen, F., Searching the Force Field Electrostatic Multipole

Parameter Space. Journal of Chemical Theory and Computation 2016, 12, 1824-1832, 2999. 26.

Verstraelen, T.; Ayers, P. W.; Van Speybroeck, V.; Waroquier, M., The

conformational sensitivity of iterative stockholder partitioning schemes. Chemical

Physics Letters 2012, 545, 138-143. 27.

Van Damme, S.; Bultinck, P.; Fias, S., Electrostatic Potentials from Self-

Consistent Hirshfeld Atomic Charges. Journal of Chemical Theory and Computation 2009, 5, 334-340. 28.

Bultinck, P.; Van Alsenoy, C.; Ayers, P. W.; Carbo-Dorca, R., Critical analysis

and extension of the Hirshfeld atoms in molecules. Journal of Chemical Physics 2007,

126, 144111. 29.

Lillestolen, T. C.; Wheatley, R. J., Atomic charge densities generated using an

iterative stockholder procedure. Journal of Chemical Physics 2009, 131, 144101.


42

Page 42 of 59

Page 43 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


30.

Verstraelen, T.; Vandenbrande, S.; Heidar-Zadeh, F.; Vanduyfhuys, L.; Van

Speybroeck, V.; Waroquier, M.; Ayers, P. W., Minimal Basis Iterative Stockholder: Atoms in Molecules for Force-Field Development. Journal of Chemical Theory and

Computation 2016, 12, 3894-3912. 31.

Ghillemijn, D.; Bultinck, P.; Van Neck, D.; Ayers, P. W., A Self-Consistent

Hirshfeld Method for the Atom in the Molecule Based on Minimization of Information Loss. Journal of Computational Chemistry 2011, 32, 1561-1567. 32.

Verstraelen, T.; Ayers, P. W.; Van Speybroeck, V.; Waroquier, M., Hirshfeld-E

Partitioning: AIM Charges with an Improved Trade-off between Robustness and Accurate Electrostatics. Journal of Chemical Theory and Computation 2013, 9, 22212225. 33.

Bultinck, P.; Cooper, D. L.; Van Neck, D., Comparison of the Hirshfeld-I and

iterated stockholder atoms in molecules schemes. Physical Chemistry Chemical

Physics 2009, 11, 3424-3429. 34.

Finzel, K.; Pendas, A. M.; Francisco, E., Efficient algorithms for Hirshfeld-I

charges. Journal of Chemical Physics 2015, 143, 084115. 35.

Manz, T. A.; Sholl, D. S., Chemically Meaningful Atomic Charges That

Reproduce the Electrostatic Potential in Periodic and Nonperiodic Materials. Journal of

Chemical Theory and Computation 2010, 6, 2455-2468. 36.

Manz, T. A.; Sholl, D. S., Improved Atoms-in-Molecule Charge Partitioning

Functional for Simultaneously Reproducing the Electrostatic Potential and Chemical States in Periodic and Nonperiodic Materials. Journal of Chemical Theory and

Computation 2012, 8, 2844-2867. 37.

Manz, T. A.; Limas, N. G., Introducing DDEC6 atomic population analysis: part 1.

Charge partitioning theory and methodology. Rsc Advances 2016, 6, 47771-47801. 38.

Marenich, A. V.; Jerome, S. V.; Cramer, C. J.; Truhlar, D. G., Charge Model 5:

An Extension of Hirshfeld Population Analysis for the Accurate Description of Molecular Interactions in Gaseous and Condensed Phases. Journal of Chemical Theory and

Computation 2012, 8, 527-541.


43


39.

Popelier, P. L. A.; Joubert, L.; Kosov, D. S., Convergence of the electrostatic

interaction based on topological atoms. Journal of Physical Chemistry A 2001, 105, 8254-8261. 40.

Dinur, U., On the interpretation of infrared intensitites in planar molecular-

systems. Chemical Physics Letters 1990, 166, 211-216. 41.

Dinur, U., Atomic multipoles and perpendicular electrostatic forces in diatomic

and planar molecules. Journal of Computational Chemistry 1991, 12, 469-486. 42.

Cioslowski, J., A new population analysis based on atomic polar tensors. Journal

of the American Chemical Society 1989, 111, 8333-8336. 43.

Dinur, U., Charge flux and electrostatic forces in planar molecules. Journal of

Physical Chemistry 1991, 95, 6201-6211. 44.

Mei, Y.; Simmonett, A. C.; Pickard, F. C.; DiStasio, R. A.; Brooks, B. R.; Shao, Y.

H., Numerical Study on the Partitioning of the Molecular Polarizability into Fluctuating Charge and Induced Atomic Dipole Contributions. Journal of Physical Chemistry A 2015, 119, 5865-5882. 45.

Unke, O. T.; Devereux, M.; Meuwly, M., Minimal distributed charges: Multipolar

quality at the cost of point charge electrostatics. Journal of Chemical Physics 2017, 147, 161712. 46.

Holt, A.; Bostrom, J.; Karlstrom, G.; Lindh, R., A NEMO Potential that Includes

the Dipole-Quadrupole and Quadrupole-Quadrupole Polarizability. Journal of


El Khoury, L.; Naseem-Khan, S.; Kwapien, K.; Hobaika, Z.; Maroun, R. G.;

Piquemal, J. P.; Gresh, N., Importance of Explicit Smeared Lone-Pairs in Anisotropic Polarizable Molecular Mechanics. Torture Track Angular Tests for Exchange-Repulsion and Charge Transfer Contributions. Journal of Computational Chemistry 2017, 38, 1897-1920. 48.

Gokcan, H.; Kratz, E.; Darden, T. A.; Piquemal, J. P.; Cisneros, G. A., QM/MM

Simulations with the Gaussian Electrostatic Model: A Density-based Polarizable Potential. Journal of Physical Chemistry Letters 2018, 9, 3062-3067.


44

Page 44 of 59

Page 45 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


49.

Gordon, M. S.; Smith, Q. A.; Xu, P.; Slipchenko, L. V., Accurate First Principles

Model Potentials for Intermolecular Interactions. In Annual Review of Physical

Chemistry, Johnson, M. A.; Martinez, T. J., Eds. 2013; Vol. 64, pp 553-578. 50.

Naserifar, S.; Brooks, D. J.; Goddard, W. A.; Cvicek, V., Polarizable charge

equilibration model for predicting accurate electrostatic interactions in molecules and solids. Journal of Chemical Physics 2017, 146, 124117. 51.

Oppenheim, J. J.; Naserifar, S.; Goddard, W. A., Extension of the Polarizable

Charge Equilibration Model to Higher Oxidation States with Applications to Ge, As, Se, Br, Sn, Sb, Te, I, Pb, Bi, Po, and At Elements. Journal of Physical Chemistry A 2018,

122, 639-645. 52.

Rappe, A. K.; Goddard, W. A., Charge equilibration for molecular-dynamics

simulations. Journal of Physical Chemistry 1991, 95, 3358-3363. 53.

Wells, B. A.; De Bruin-Dickason, C.; Chaffee, A. L., Charge Equilibration Based

on Atomic Ionization in Metal-Organic Frameworks. Journal of Physical Chemistry C 2015, 119, 456-466. 54.

Wilmer, C. E.; Kim, K. C.; Snurr, R. Q., An Extended Charge Equilibration

Method. Journal of Physical Chemistry Letters 2012, 3, 2506-2511. 55.

Martin-Noble, G. C.; Reilley, D.; Rivas, L. M.; Smith, M. D.; Schrier, J., EQeq plus

C: An Empirical Bond-Order-Corrected Extended Charge Equilibration Method. Journal


Mortier, W. J.; Ghosh, S. K.; Shankar, S., Electronegativity equalization method

for the calculation of atomic charges in molecules. Journal of the American Chemical

Society 1986, 108, 4315-4320. 57.

Chelli, R.; Procacci, P.; Righini, R.; Califano, S., Electrical response in chemical

potential equalization schemes. Journal of Chemical Physics 1999, 111, 8569-8575. 58.

Nistor, R. A.; Polihronov, J. G.; Muser, M. H.; Mosey, N. J., A generalization of

the charge equilibration method for nonmetallic materials. Journal of Chemical Physics 2006, 125, 094108. 59.

Dapp, W. B.; Muser, M. H., Towards time-dependent, non-equilibrium charge-

transfer force fields. European Physical Journal B 2013, 86, 337.


45


60.

Smirnov, K. S., Assessment of split-charge equilibration model for development

of polarizable force fields. Modelling and Simulation in Materials Science and

Engineering 2015, 23. 61.

Chen, J. H.; Martinez, T. J., QTPIE: Charge transfer with polarization current

equalization. A fluctuating charge model with correct asymptotics. Chemical Physics

Letters 2007, 438, 315-320. 62.

Chen, J. H.; Martinez, T. J., QTPIE: Charge transfer with polarization current

equalization. A fluctuating charge model with correct asymptotics (vol 438, pg 315, 2007). Chemical Physics Letters 2008, 463, 288-288. 63.

York, D. M.; Yang, W. T., A chemical potential equalization method for molecular

simulations. Journal of Chemical Physics 1996, 104, 159-172. 64.

Bret, C.; Field, M. J.; Hemmingsen, L., A chemical potential equalization model

for treating polarization in molecular mechanical force fields. Molecular Physics 2000,

98, 751-763. 65.

Verstraelen, T.; Ayers, P. W.; Van Speybroeck, V.; Waroquier, M., ACKS2: Atom-

condensed Kohn-Sham DFT approximated to second order. Journal of Chemical

Physics 2013, 138, 074108. 66.

Parr, R. G.; Pearson, R. G., Absolute hardness - companion parameters to

absolute electronegativity. Journal of the American Chemical Society 1983, 105, 75127516. 67.

Geerlings, P.; De Proft, F.; Langenaeker, W., Conceptual density functional

theory. Chemical Reviews 2003, 103, 1793-1873. 68.

Stern, H. A.; Kaminski, G. A.; Banks, J. L.; Zhou, R. H.; Berne, B. J.; Friesner, R.

A., Fluctuating charge, polarizable dipole, and combined models: Parameterization from ab initio quantum chemistry. Journal of Physical Chemistry B 1999, 103, 4730-4737. 69.

Davari, N.; Haghdani, S.; Astrand, P. O.; Schatz, G. C., Local electric field factors

by a combined charge-transfer and point-dipole interaction model. Rsc Advances 2015,

5, 31594-31605. 70.

Jensen, L.; Astrand, P. O.; Osted, A.; Kongsted, J.; Mikkelsen, K. V.,

Polarizability of molecular clusters as calculated by a dipole interaction model. Journal

of Chemical Physics 2002, 116, 4001-4010.


46

Page 46 of 59

Page 47 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


71.

Wang, L. P.; Chen, J. H.; Van Voorhis, T., Systematic Parametrization of

Polarizable Force Fields from Quantum Chemistry Data. Journal of Chemical Theory

and Computation 2013, 9, 452-460. 72.

Nalewajski, R. F.; Korchowiec, J.; Zhou, Z. X., Molecular hardness and softness

parameters and their use in chemistry. International Journal of Quantum Chemistry 1988, 349-366. 73.

Oda, A.; Hirono, S., Geometry-dependent atomic charge calculations using

charge equilibration method with empirical two-center Coulombic terms. Journal of

Molecular Structure-Theochem 2003, 634, 159-170. 74.

Chen, J. H.; Martinez, T. J., Charge conservation in electronegativity equalization

and its implications for the electrostatic properties of fluctuating-charge models. Journal

of Chemical Physics 2009, 131, 044114. 75.

Warren, G. L.; Davis, J. E.; Patel, S., Origin and control of superlinear

polarizability scaling in chemical potential equalization methods. Journal of Chemical

Physics 2008, 128, 144110. 76.

Verstraelen, T.; Bultinck, P.; Van Speybroeck, V.; Ayers, P. W.; Van Neck, D.;

Waroquier, M., The Significance of Parameters in Charge Equilibration Models. Journal


Krykunov, M.; Demone, C.; Lo, J. W. H.; Woo, T. K., A New Split Charge

Equilibration Model and REPEAT Electrostatic Potential Fitted Charges for Periodic Frameworks with a Net Charge. Journal of Chemical Theory and Computation 2017, 13, 2858-2869. 78.

Chen, J. H.; Hundertmark, D.; Martinez, T. J., A unified theoretical framework for

fluctuating-charge models in atom-space and in bond-space. Journal of Chemical

Physics 2008, 129, 214113. 79.

Mathieu, D., Split charge equilibration method with correct dissociation limits.

Journal of Chemical Physics 2007, 127, 224103. 80.

Mikulski, P. T.; Knippenberg, M. T.; Harrison, J. A., Merging bond-order

potentials with charge equilibration. Journal of Chemical Physics 2009, 131, 241105.


47


81.

Bauer, B. A.; Patel, S., Recent applications and developments of charge

equilibration force fields for modeling dynamical charges in classical molecular dynamics simulations. Theoretical Chemistry Accounts 2012, 131, 1153. 82.

Wang, Q.; Rackers, J. A.; He, C.; Qi, R.; Narth, C.; Lagardere, L.; Gresh, N.;

Ponder, J. W.; Piquemal, J.-P.; Ren, P., General Model for Treating Short-Range Electrostatic Penetration in a Molecular Mechanics Force Field. Journal of Chemical

Theory and Computation 2015, 11, 2609-2618. 83.

Rackers, J. A.; Wang, Q.; Liu, C.; Piquemal, J.-P.; Ren, P.; Ponder, J. W., An

optimized charge penetration model for use with the AMOEBA force field. Physical

Chemistry Chemical Physics 2017, 19, 276-291. 84.

Gresh, N.; Stevens, W. J.; Krauss, M., Mono-ligated and poly-ligated complexes

of Zn2+ - an ab initio analysis of the metal-lgand interaction energy. Journal of


Maxwell, J., A treatise on electricity and magnetism. Clarendon Press: 1873.

86.

Liang, T.; Shin, Y. K.; Cheng, Y. T.; Yilmaz, D. E.; Vishnu, K. G.; Verners, O.;

Zou, C. Y.; Phillpot, S. R.; Sinnott, S. B.; van Duin, A. C. T., Reactive Potentials for Advanced Atomistic Simulations. In Annual Review of Materials Research, Vol 43, Clarke, D. R., Ed. 2013; Vol. 43, pp 109-129. 87.

Holt, A.; Karlstrom, G.; Roos, B. O., The Charge Capacitance of the Chemical

Bond: Application to Bonds Containing Metals. International Journal of Quantum

Chemistry 2009, 109, 618-628. 88.

Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.;

Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; J. A. Montgomery, J.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.;


48

Page 48 of 59

Page 49 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox;, D. J.; Gaussian Inc. Gaussian-09 rev.D-01, 2009. 89.

Chai, J.-D.; Head-Gordon, M., Long-range corrected hybrid density functionals

with damped atom-atom dispersion corrections. Physical Chemistry Chemical Physics 2008, 10, 6615-6620. 90.

Jensen, F., Unifying General and Segmented Contracted Basis Sets. Segmented

Polarization Consistent Basis Sets. Journal of Chemical Theory and Computation 2014,

10, 1074-1085. 91.

Dunning, T. H., Gaussian-basis sets for use in correlated molecular calculations.

1. The atoms boron through neon and hydrogen. Journal of Chemical Physics 1989, 90, 1007-1023. 92.

Lu, T.; Chen, F. W., Multiwfn: A multifunctional wavefunction analyzer. Journal of


Mahoney, M. W.; Jorgensen, W. L., A five-site model for liquid water and the

reproduction of the density anomaly by rigid, nonpolarizable potential functions. The

Journal of Chemical Physics 2000, 112, 8910-8922. 94.

Holt, A.; Karlstrom, G.; Lindh, R., The charge capacity of the chemical bond.

Chemical Physics Letters 2007, 436, 297-301. 95.

Francl, M. M.; Chirlian, L. E., The pluses and minuses of mapping atomic

charges to electrostatic potentials. Reviews in Computational Chemistry 2000, 14, 1-31. 96.

Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A., A well-behaved

electrostatic potential based method using charge restraints for deriving atomic charges - the RESP model. Journal of Physical Chemistry 1993, 97, 10269-10280. 97.

Sifain, A. E.; Lubbers, N.; Nebgen, B. T.; Smith, J. S.; Lokhov, A. Y.; Isayev, O.;

Roitberg, A. E.; Barros, K.; Tretiak, S., Discovering a Transferable Charge Assignment Model Using Machine Learning. Journal of Physical Chemistry Letters 2018, 9, 44954501. 98.

Nebgen, B.; Lubbers, N.; Smith, J. S.; Sifain, A. E.; Lokhov, A.; Isayev, O.;

Roitberg, A. E.; Barros, K.; Tretiak, S., Transferable Dynamic Molecular Charge


49


Assignment Using Deep Neural Networks. Journal of Chemical Theory and

Computation 2018, 14, 4687-4698.


50

Page 50 of 59

Page 51 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Illustrating the position of dummy atoms for water, ammonia and planar molecules. The distance to the dummy atom is in all cases taken as 1 Bohr. The molecular plane for water is xz, with dummy atoms placed in the zy-plane with an XOX angle of 97.8 while for ammonia the dummy atom is placed along the C3 axis. Structures A, B, C show the placement of dummy atoms for planar molecules/groups as benzene and ethane, carbonyl groups and boranes, respectively. 151x82mm (90 x 90 DPI)



Figure 2. The fluoromethane dimer for probing the intramolecular polarization shown in Figure 3.


Page 52 of 59

Page 53 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Change in the dipole moment for monomer A (Figure 2) as a function of intermolecular distance with results aligned at 13 Å 500x345mm (72 x 72 DPI)



Figure 4. Illustrative CH3F dimer geometries used to probe intermolecular polarization shown in Figure 5 as a function of angle.


Page 54 of 59

Page 55 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. Change in the x- and z-dipole components with zero taken at  = 0 for monomer A with the geometry variation shown in Figure 4. 500x345mm (72 x 72 DPI)



Figure 5. Change in the x- and z-dipole components with zero taken at  = 0 for monomer A with the geometry variation shown in Figure 4. 500x345mm (72 x 72 DPI)


Page 56 of 59

Page 57 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. 1-Fluoro-decane conformations used for probing intramolecular polarization shown in Figure 7 with the carbon backbone laying in the xy-plane. 380x169mm (90 x 90 DPI)



Figure 7. Electrostatic potential in atomic units in the xy-plane arising from the atomic charges difference for the C5-C10 atoms and bonded hydrogens for 1-fluoro-decane (Figure 6) for BC charges with Jii = 0, Jii = 1 and HI.


Page 58 of 59

Page 59 of 59 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure A1. Principal polarizability component calculated with the BC model for a system of up to 30 CO molecules aligned along the z-axis, with an interfragment distance of 3 Å and a bond capacity of 3.273 a.u. 500x345mm (72 x 72 DPI)


Describing Molecular Polarizability by a Bond Capacity Model

Recommend Documents