Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2922695

Formats

Article sections

- Abstract
- 1. Introduction
- 2. Method
- 3. Computational Details
- 4. Results and Discussion
- 5. Conclusions
- Supplementary Material
- References

Authors

Related links

J Chem Theory Comput. Author manuscript; available in PMC 2011 January 1.

Published in final edited form as:

J Chem Theory Comput. 2010; 6(8): 2469–2476.

doi: 10.1021/ct100268pPMCID: PMC2922695

NIHMSID: NIHMS223806

Alessandro Cembran, Peng Bao, Yingjie Wang, Lingchun Song,^{*} Donald G. Truhlar,^{*} and Jiali Gao^{*}

Department of Chemistry and Minnesota Supercomputing Institute University of Minnesota, Minneapolis, MN 55455

See other articles in PMC that cite the published article.

The inclusion of exchange repulsion terms in the explicit polarization (X-Pol) model is examined by antisymmetrizing the X-Pol Hartree-product wave function; this yields XPol with full eXchange, called X-Pol-X. When the monomers are treated by Hartree-Fock theory, this calculation can be accomplished by using the formalism of block-localized wave functions (BLW) that has been used in a variety of applications. In this case the block-localized structure in the X-Pol-X wave function allows for decomposition of the full Fock matrix of a dimension of *M* blocks into *M* smaller Fock matrices. The method is illustrated by considering two trimer structures of water clusters, and it is found that the total exchange repulsion energies in these hydrogen-bonding test cases are adequately treated and—to a good approximation— are pairwise additive. We also present a formalism to yield a simplified Fock matrix by making use of the neglect of interfragment differential overlap (NIDO) approximation, which is less severe than the neglect of diatomic differential overlap (NDDO) approximation.

The explicit polarization (X-Pol) method,^{1} initially called molecular orbital derived potentials for liquids (MODEL), was introduced in 1997^{2} for the treatment of macromolecular systems in statistical mechanical Monte Carlo and molecular dynamics simulations using electronic structure theory.^{3}^{–}^{5} It is based on the partition of a whole system into subsystems, which are called fragments or blocks, each of which is represented by an antisymmetrized wave function.^{2}^{–}^{3} The system wave function is then a Hartree product of fragment wave functions, which are optimized by the self-consistent field (SCF) method in the presence of the external field of all other blocks until the energy or electron density of the entire system is converged.^{2}^{–}^{3}^{,}^{6}^{– }^{8} The individual wave functions of the subsystems can be obtained at any level of theory—ab initio Hartree-Fock (HF), semiempirical molecular orbital theory, correlated wave function theory, or Kohn-Sham (KS) density functional theory (DFT).^{9} In general, the X-Pol method can be applied to post-SCF calculations such MP2 and multiconfigurational methods, and a post-SCF treatment can be used to estimate dispersion interaction energies between different fragments. In contrast to the methods of Stoll and Preuss^{10}^{–}^{11} (whose emphasis was on the localization of orbitals on atomic centers with the neighboring bonding partners and lone-pairs) and Yang^{12} (whose emphasis was on efficient full calculations for large systems by a divide-and-conquer approach), X-Pol is not designed to reproduce the full HF or KS energy; rather it is intended to establish a theoretical framework for economical and accurate applications to macromolecular systems.^{3}^{,}^{5}

The Hartree-product approximation of the X-Pol wave function implies that interfragment exchange and correlation interactions are neglected.^{1}^{–}^{3}^{,}^{6}^{,}^{8} The most noticeable contributions to the interfragment exchange term are from the repulsion interactions at geometries where orbitals on different blocks have significant overlap (short-range), whereas the dominant contributions to interfragment correlation energy are dispersion interactions at geometries where they have negligible overlap (long-range). To achieve computational efficiency, we have proposed to use empirical Lennard-Jones terms to account for these interactions in the spirit of developing a quantal force field for macromolecular simulations.^{1}^{–}^{3}^{,}^{5}^{–}^{6} Alternatively, interfragment exchange can be included in the X-Pol method by fully antisymmetrizing the block-localized molecular orbitals (BLMOs) to form a determinant wave function, and this is equivalent to the block-localized wave function (BLW) method described by Mo,^{13} which has been applied to a variety of problems.^{9}^{,}^{14}^{–}^{25} In the present context, this approach will be called X-Pol with full eXchange or X-Pol-X. In such a wave function the orbital coefficient matrix is block-diagonal, which has the property that orbitals within each block are orthogonal, and orbitals between different subsystems are non-orthogonal.^{11}^{,}^{13}^{–}^{14}^{,}^{26} The latter is characteristic of valence bond theory.^{9}^{,}^{16}^{,}^{23}^{,}^{27} Methods for treating non-orthogonal determinantal wave functions have been described a long time ago and used in various applications, including the work of King et al.,^{28} Gallup,^{29} Stoll,^{11}^{,}^{30} Gianinetti and Raimondi,^{31}^{–}^{32} Mo et al.,^{13}^{–}^{16}^{,}^{27} Head- Gordon and coworkers,^{33} and Cembran et al.^{9} The BLW method has been used to define valence bond-like states in the mixed molecular orbital and valence bond (MOVB) theory^{16}^{,}^{22}^{–}^{23}^{,}^{27} and in the multistate density functional theory (MSDFT) based on valence bond or VBDFT.^{9} Here, we use the X-Pol-X (*i.e*., BLW) method to elucidate the magnitude of the exchange repulsion term that is required in the X-Pol method.^{15} In the Appendix we derive a simplified X-Pol-X method by invoking neglect of interfragment differential overlap (NIDO).

It should be noted that in 1999, Kitaura et al.^{34}^{–}^{35} used the wave functions of dimeric fragments in a post SCF correction to estimate the pair-wise exchange repulsion. The expansion of the monomeric subspace into two fragments also includes pair-wise charge transfer effects which are important for describing short-range interactions.^{10} This approach was called the fragment molecular orbital (FMO) method, and it has been used in numerous applications with various electronic structural models.^{36} However, aside from the dimer correction term, first described by Stoll and Preuss,^{10} the construction of block-localized (fragment) orbitals and the SCF computational algorithm (see page 321 of ref. ^{34}) of the FMO method^{34}^{–}^{35} are essentially the same as that (see page 660 of ref ^{2}) in X-Pol theory.^{2}^{–}^{3} Nevertheless, these studies illustrate the usefulness of block localization of molecular orbitals for macromolecular systems.

Section 2 presents the X-Pol method in a general way and then describes the X-Pol-X extension of the X-Pol method to include the exchange repulsion energy in SCF calculations. Section 3 describes an illustrative application at the ab initio Hartree-Fock level. Section 4 gives the results and discusses them, with special emphasis on comparison of X-Pol to X-Pol-X. Section 5 summarizes the conclusions. The NIDO simplification procedure, which may be very useful in future work, is presented in an Appendix.

We first describe the X-Pol method. Then, we present the SCF procedure for optimization of an X-Pol-X wave function. It is especially interesting to consider the case where the fragment wave functions are Hartree-Fock wave functions, obtained either by ab initio theory or by semiempirical molecular orbital theory. Then the fragment wave functions are optimized by the self-consistent field (SCF) method in the presence of the external field of all other blocks.^{2}^{–}^{3}^{,}^{6}^{–}^{8} This is sometimes called the double SCF (DSCF) procedure because one must achieve self-consistency both within each fragment and among all blocks in the system. As a consequence of the interfragment self-consistency, the fragment wave functions are polarized.

We partition a large molecular system into *M* blocks (which have been equivalently and interchangeably called subgroups, subsystems, molecules, residues, monomers, fragments, or bodies). The term block is more general than some other terms, such as fragment and body, in that orbitals on the same atom can be divided into different subsystems without the physical constraint of atomic structure fragmentation.^{14} Block *a* contains *k _{a}* basis functions and

$$K=\sum _{a=1}^{M}{k}_{a}\phantom{\rule{0.38889em}{0ex}}\text{and}\phantom{\rule{0.38889em}{0ex}}N=\sum _{a=1}^{M}{n}_{a}$$

(1)

The molecular orbitals
${\varphi}_{j}^{a}$ in each block *a* are written as linear combinations of the primitive basis functions {
${\chi}_{\mu}^{a}$; *μ* = 1, ···, *k _{a}*} in that specific subspace:

$${\varphi}_{j}^{a}=\sum _{\mu =1}^{{k}_{a}}{c}_{j\mu}^{a}{\chi}_{\mu}^{a}$$

(2)

Thus, by construction, the molecular orbitals are strictly localized within the subspace of each block.

The X-Pol wave function is approximated as a product of the individual, antisymmetric fragment wave functions of all *M* blocks:

$${\mathrm{\Psi}}^{\text{X}-\text{Pol}}=\prod _{a=1}^{M}\widehat{A}{\mathrm{\Phi}}_{a}$$

(3)

where *Â* is an antisymmetrization operator, and Φ* _{a}* is a product of the occupied spin-orbitals in block

$${\mathrm{\Phi}}_{a}={\varphi}_{1}^{a}{\varphi}_{2}^{a}\cdots {\varphi}_{{n}_{a}}^{a}$$

(4)

The Hartree-product wave function in eq 3 includes Coulomb interactions, but the exchange interactions are ignored between different blocks. If we make a further approximation by representing the Coulomb potentials by a classical multipole expansion or by the potential due to distributed monopoles calculated by Mulliken population analysis, the computation of the Fock matrix for block *a*, in the presence of the external field of the rest of the system, is simplified to

$${\mathbf{F}}^{a}={\mathbf{H}}^{a}+\sum _{i}({\mathbf{J}}_{i}^{a}-\frac{1}{2}{\mathbf{K}}_{i}^{a})+\sum _{b\ne a}{\mathbf{V}}_{b}^{a}$$

(5)

where **H*** ^{a}* is the one-electron Hamiltonian matrix,
${\mathbf{J}}_{i}^{a}$ and
${\mathbf{K}}_{i}^{a}$ are the Coulomb and exchange integral matrices associate with molecular orbital

We note that the electrostatic potential
${\mathbf{V}}_{b}^{a}$, which is derived from the wave function of block *b*, can be calculated in more than one way. For example, it can be determined by integration over the BLMOs in that block, it can be represented by fitted atomic densities, or it can be approximated by electrostatic-potential-fitted charges or by a multipole expansion.^{2}^{,}^{6}^{,}^{37} These choices may be dictated by the context in which X-Pol is used. In previous applications, we have used charges obtained by Mulliken population analysis to approximate the Coulomb potentials of each of the monomers.^{1}^{–}^{2}^{,}^{5}^{–}^{6}^{,}^{8} In the present work, we use the explicit evaluation of all two-electron Coulomb integrals to account for inter-block electrostatic interactions. Therefore, the calculation of Coulomb interactions between different monomer blocks does not introduce approximations here, and this makes it easier to evaluate the new theoretical approach to exchange repulsion.

The short-range exchange repulsion energies between different blocks can be approximated empirically by Lennard-Jones terms,
${\epsilon}_{ab}^{LJ}$, which also account for long-range dispersion interactions.^{2}^{–}^{3} Thus, the total X-Pol energy of the system is

$${E}^{\text{X}-\text{Pol}}={E}^{\text{SCF}}+\frac{1}{2}\sum _{a}\sum _{b\ne a}{\epsilon}_{ab}^{LJ}$$

(6)

The use of the *R*^{−12} repulsive terms in the Lennard-Jones potential (where *R* is a distance between atoms in different monomers) to approximate the short-range exchange repulsion significantly reduces the computational cost needed to solve the Hartree-Fock- Roothaan equations for the entire block-localized (*i.e*., fragmental) system. Alternatively, the exchange energies can be determined without the use of an empirical potential by antisymmetrizing the X-Pol wave function, keeping the block-localized structure of the molecular orbitals:^{13}^{–}^{16}^{,}^{27}

$${\mathrm{\Psi}}^{\text{X}-\text{Pol}-\text{X}}=\widehat{A}\prod _{a=1}^{M}{\mathrm{\Phi}}_{a}$$

(7)

The coefficient matrix of the corresponding orbitals has the following block-diagonal form:

$$\mathbf{C}=\left(\begin{array}{cccc}{\mathbf{C}}^{\mathbf{1}}& 0& \dots & 0\\ 0& {\mathbf{C}}^{\mathbf{2}}& \dots & 0\\ \dots & \dots & \dots & \dots \\ 0& 0& \dots & {\mathbf{C}}^{\mathbf{k}}\end{array}\right)$$

(8)

where **C*** ^{a}* is an

An alternative to the expression of eq 7 is to use the antisymmetrized function of antisymmetrized block-localized wave functions of eq 3. However, in the present case, the molecular orbitals are strictly block-localized within each fragment space, and the difference from eq 7 is likely to be minimal (antisymmetrying determinant states is useful if each of the block-localized fragment is used as a reference to form multiconfigurational states).^{39} In addition, the present test results show that the exchange repulsion energy is adequately treated, suggesting that eq 7 in the current form is a good choice for the purpose of this study.

The molecular orbitals in eq 7 are subject to the constraint that orbitals within the same subgroup are orthogonal, whereas orbitals in different subgroups are non-orthogonal:

$$<{\varphi}_{i}^{a}\mid {\varphi}_{j}^{b}>=\{\begin{array}{ll}{\delta}_{ij}^{aa},\hfill & a=b\hfill \\ {S}_{ij}^{ab},\hfill & a\ne b\hfill \end{array}$$

(9)

where
${S}_{ij}^{ab}$ is the overlap integral between molecular orbitals *i* and *j* from blocks *a* and *b*, respectively, and is given by

$${S}_{ij}^{ab}={({\mathbf{c}}_{i}^{a})}^{\top}{\mathbf{R}}^{ab}{\mathbf{c}}_{j}^{b}$$

(10)

where {
${\mathbf{c}}_{i}^{a}$; *i* = 1, ···, *n _{a}*,

$${R}_{\mu \nu}^{ab}=\langle {\chi}_{\mu}^{a}|{\chi}_{\nu}^{b}\rangle $$

(11)

The orthogonality constraint on orbitals in the same block does not affect the total energy of the system because the energy is invariant to a unitary transformation of the orbitals within a block.

The density matrix is defined, using the occupied orbitals, which include the column vectors {
${\mathbf{c}}_{i}^{a}$} for block **C*** ^{a}* in eq 8. by

$$\mathbf{D}=\mathbf{C}{({\mathbf{C}}^{\top}\mathbf{RC})}^{-1}{\mathbf{C}}^{\top}$$

(12)

The Fock matrix, **F**, has an identical form to that in Hartree-Fock theory, in which an element, in the atomic basis {*χ _{μ}*}, is given for a closed-shell system as

$${F}_{\mu \nu}={H}_{\mu \nu}+\sum _{\lambda \sigma}^{K}{D}_{\lambda \sigma}\{({\chi}_{\mu}{\chi}_{\nu}\mid {\chi}_{\lambda}{\chi}_{\sigma})-\frac{1}{2}({\chi}_{\mu}{\chi}_{\lambda}\mid {\chi}_{\nu}{\chi}_{\sigma})\}$$

(13)

where *H _{μν}* is an element of the one-electron matrix representing kinetic energy and nuclear attraction, and the brackets are the Coulomb and exchange integrals in the Mulliken notation. The total X-Pol-X energy is

$${E}^{\text{X}-\text{Pol}-\text{X}}=\sum _{\mu \nu}^{K}{D}_{\mu \nu}\{{H}_{\mu \nu}+{F}_{\mu \nu}\}$$

(14)

In comparison with the X-Pol SCF energy *E*^{SCF} (eq 6), the exchange energy included in the XPol- X expression (eq 14) originates from the orthogonalization constraint **(C**^{T}**RC)**^{−1} in eq 12.^{15}^{,}^{40}

Since the orbital coefficient matrix is block diagonal (eq 8), the *K × K* Fock matrix can be reduced to smaller sizes, one for each block (*k _{a} × k_{a}*) by transforming the original Fock matrix into a block-diagonal form using a projection operator.

$${\mathfrak{J}}^{aa}={\mathbf{F}}^{aa}+{V}_{2}+{V}_{3}+{V}_{4}$$

(15)

where

$${V}_{2}=-\sum _{b\ne a}^{k}{\mathbf{F}}^{ab}\left(\sum _{c\ne a}^{k}{\mathbf{D}}^{bc}{\mathbf{R}}^{ca}\right)$$

(16)

$${V}_{3}=-\sum _{b\ne a}^{k}{\left(\sum _{c\ne a}^{k}{\mathbf{D}}^{bc}{\mathbf{R}}^{ca}\right)}^{\top}{\mathbf{F}}^{ba}$$

(17)

and

$${V}_{4}=\sum _{d\ne a}^{k}\sum _{b\ne a}^{k}{\left(\sum _{c\ne a}^{k}{\mathbf{D}}^{bc}{\mathbf{R}}^{ca}\right)}^{\top}{\mathbf{F}}^{bd}\left(\sum _{c\ne a}^{k}{\mathbf{D}}^{dc}{\mathbf{R}}^{ca}\right)$$

(18)

The upper indices in eqs 15–18 are used to indicate that the dimension of each matrix with such indices is determined by the number of basis functions in that block (see eq 1). Note that eq 15 is general in that it can be used both in Hartree-Fock^{13}^{–}^{16}^{,}^{27} and DFT calculations.^{9} In the latter case, the exchange potential in Hartree-Fock theory is replaced by the exchangecorrelation potential.^{9}

All computations have been performed using a locally modified GAMESS program^{41} and the Xiamen University Valence Bond (XMVB) program.^{42} The minimally augmented polarized valence double-zeta basis set 6-31+G(d) was used for all calculations. All X-Pol-X calculations are fully self-consistent, that is, the orbitals are optimized in the presence of interfragment exchange.

We present illustrative calculations for two examples, both of which are structures of water trimer. First, we chose the minimum energy configuration of a cyclic water trimer structure, **c-W _{3}**, in which each water molecule donates and accepts one hydrogen bond from one of the other two water monomers.

Illustration of the minimum water trimer structure (**c-W**_{3}) in (a) and a symmetric configuration (**s-W**_{3}) in (b).

In the following discussion, each water monomer in the trimer complexes is treated as a single block or fragment, and the internal geometries of the monomers are kept the same as in the fully optimized configuration at the HF/6-31+G(d) level for **c-W _{3}** and the same as constructed above for

The energy components we consider are the Coulomb energy, the exchange energy, the polarization energy, and the charge transfer energy. The Coulomb and exchange components are the same in all energy decomposition analysis (EDA) methods, including the approach of Kitaura and Morokuma,^{44} whereas the polarization and charge transfer terms are defined in the block-localized wave-function energy decomposition (BLW-ED) method developed by Mo et al.^{15}

The main goal of the present article is to show that the theory presented in Section 2 can provide a reasonable estimate of the exchange repulsion energy in the X-Pol method.

Table S1 (see supporting information) lists the total energies and energy components determined for the fully separated monomers and for the trimer complexes as calculated from *Â ψ _{a} ψ_{b} ψ_{c}* and from
$\widehat{A}{\psi}_{a}^{0}{\psi}_{b}^{0}{\psi}_{c}^{0}$ and also by full Hartree-Fock calculations. In this discussion we shall focus on interaction energies, defined as the energy of the trimer minus the sum of the energies of the three separated monomers, each with the same geometry that it has in the trimer. We shall also discuss the components of the interactions. The interaction energies and their components are given in Table 1; all values in Table 1 were computed from the absolute energies in Table S1.

Computed relative energies for the cyclic water trimer minimum structure (**c-W**_{3}) and for a symmetric trimer geometry (**s-W**_{3}) (kcal/mol). All energies are determined using the 6-31+G(d) basis set.

Row 1 of Table 1 gives the X-Pol-X interaction energy, and the following rows of the table report quantities that help us dissect this interaction energy. Row 2 is calculated using
${\psi}_{a}^{0}{\psi}_{b}^{0}{\psi}_{c}^{0}$, and row 3 is calculated using
$\widehat{A}{\psi}_{a}^{0}{\psi}_{b}^{0}{\psi}_{c}^{0}$. Row 2 gives the unpolarized Coulomb contribution to the X-Pol or X-Pol-X binding energy. The results in Table 2 show that the Coulomb interaction to the total binding energy is quite large in both water trimer complexes. In particular, with the monomers unpolarized, the total Coulomb interaction energies are −25.6 and −16.0 kcal/mol for **c-W _{3}** and

Additivity of dimer contributions to the Coulomb and exchange-repulsion interaction energies in the cyclic water trimer minimum structure (**c-W**_{3}) and in a symmetric trimer geometry (**s-W**_{3}) (kcal/mol). All energies are determined using the 6-31+G(d) basis **...**

Next we consider the analogous breakdown of the X-Pol-X interaction energies and interaction energy components. For these calculations, the wave function is *Â ψ _{a} ψ_{b} ψ_{c}* where the antisymmetrization is performed before the orbitals are optimized, so that they are optimized in the presence of interfragment exchange. Now, as compared to rows 2–4, there is an additional component of the total interaction energy, namely the increase in the internal energy of the monomers due to the polarization of the orbitals. This is given in row 5, the Coulomb energy is in row 6, and the exchange energy is in row 7. The sum of rows 5–7 gives the X-Pol-X interaction energy in row 1. Note that row 4 is the total energy of the antisymmetrized fragment product wave function where the orbitals are the unpolarized monomer orbitals, whereas row 1 is the total energy for the X-Pol-X wave function, in which the mutual charge and exchange polarizations of the monomer charge density are included.

Mutual polarization of the block-localized monomer wave functions in the presence of other monomers enhances Coulomb interactions by −6.2 and −2.3 kcal/mol, but reduces the exchange repulsion slightly. Taking into account of the energy costs for polarizing the monomer wave functions in the trimer complex and the small decrease in exchange repulsion due to polarization effects, we see that the overall polarization, calculated at this level, lowers the **c-W _{3}** interaction energy from −9.3 to −12.5 kcal/mol and the

Electron density difference contour between the polarized and unpolarized systems,
$\mathrm{\Delta}\rho =\rho (\widehat{A}{\psi}_{a}{\psi}_{b}{\psi}_{c})-{\rho}^{o}({\psi}_{a}^{o}{\psi}_{b}^{o}{\psi}_{c}^{o})$, for (a) the minimum-energy structure (**c-W3**) and (b) **...**

Row 10 depicts the interaction energy from a full (delocalized) Hartree-Fock calculation. The energy difference between the Hartree-Fock SCF result and the X-Pol-X energy for each trimer complex is defined as the charge-transfer energy, which is not included in the X-Pol-X method and will be addressed in a separate study. This charge transfer energy (row 9) is −3.0 and −2.1 kcal/mol using the 6-31+G(d) basis functions in the cyclic and symmetric complexes, respectively. For comparison, Chen and Gordon found a charge transfer energy of −1.7 kcal/mol using the RVS approach, smaller than the present value.^{46} The BLW energy decomposition analysis (BLW-ED) method has been shown to be relatively stable in computed polarization and charge transfer energies with increased size of basis set.^{15}^{,}^{18}^{,}^{47}

It is interesting to further analyze the exchange repulsion energy of the X-Pol-X method. The numerical results in Table 1 show that the exchange repulsion is significant, reducing the Coulomb interaction energies by more than 50% in both cases. Table 2 shows the pair-wise contributions to the total Coulomb and exchange energies. The variation of Coulomb energies for the dimer interactions depicted in the first three rows in Table 2 suggests that stronger Coulomb attractions are accompanied by greater charge overlap and exchange repulsions. The Coulomb interaction energy is additive by definition in the energy decomposition analysis,^{15}^{,}^{40}^{,}^{44} but the exchange repulsion term is not necessarily additive due to orthogonalization of the product of nonorthogonal block localized orbitals. Nevertheless, Table 2 reveals that the exchange repulsion energies are additive to within 0.03 kcal/mol in both trimer complexes. This result holds even though one complex (**c-W _{3}**) has three significant pairwise exchange repulsion interactions, while the other (

The dependence of exchange repulsion on intermolecular separation is depicted in Figure 3 for a water dimer as a function of the hydrogen-bond distance between the donor hydrogen and acceptor oxygen atoms. As is well-known, the exchange repulsion decreases exponentially to a diminishing value as the hydrogen bond distance is elongated by more than 1 Å from the optimal geometry, well within the first solvation layer in liquid water. Figure 3 shows that for all practical purposes, one only needs to consider the explicit exchange repulsion for monomers within the first solvation shell, and there is no need to construct a fully antisymmetric X-Pol determinant wave function for the entire system.^{47} Further simplifications of the X-Pol-X method can be made by imposing the NIDO approximation presented in the Appendix.

Computed Coulomb (X-Pol) and exchange repulsion (X-Pol-X) interaction energy as a function of the hydrogen bond distance between the donor hydrogen and the acceptor oxygen atoms in the water dimer complex with C_{s} symmetry. Energies are computed using **...**

Several topics suggest themselves for further work. First, it would be useful to examine more systems. Second, the method can be used to account for charge transfer effects.^{9}^{,}^{16}^{,}^{23}^{,}^{27}^{,}^{47} Third, the method may be examined by treating the monomers with density functional theory.^{6}^{,}^{9} Fourth, it would be interesting to examine the effect of computing the interfragment exchange terms of X-Pol-X in a post-SCF fashion, that is, optimizing the orbitals by X-Pol using a constant interfragment exchange based on the unpolarized fragmental wave functions and then using the optimized BLMOs to compute the exchange repulsion. Fifth, the NIDO approximation presented in the Appendix can be used with either full X-Pol-X calculations, including the use of full NDDO approximation in effective-Hamiltonian MOVB method,^{25} or with post-SCF interfragment exchange. Finally, it would be interesting to use the interfragment exchange terms of X-Pol-X calculations as a basis for parameterizing more accurate analytical pairwise exchange repulsion terms.

In this article, we have examined the inclusion of exchange repulsion energies directly into the explicit polarization (X-Pol) model by antisymmetrizing the Hartree-product wave function in a way that is equivalent to the block-localized wave function (BLW) approach. In general, this may be called X-POL with full exchange or X-Pol-X. The block-localized structure in the X-Pol-X wave function allows for decomposition of the full Fock matrix of *M* blocks into *M* smaller Fock matrices, reducing the computational costs from *Q*[(*KM*)^{3}] scaling to *Q*[*M*(*K*)^{3}]. Furthermore, the complexity of the Fock matrix can be reduced by making use of the neglect fragmental diatomic overlap approximation, a modification of the traditional neglect diatomic differential overlap (NDDO) method. Another way to reduce the cost would be to compute the exchange repulsion by a set of pairwise X-Pol-X dimer calculations, including the contributions of the nearby pairs. The computation of the exchange repulsion energy is illustrated by considering two trimer structures of water, and it was found that the total exchange repulsion energies in these simple, but general and strongly hydrogen-bonding cases are fully pair-wise additive.

We are grateful to Peng Zhang, Hannah R. Leverentz, and Carlos Sosa for valuable discussions. This work has been partially supported by the National Institutes of Health (grant no. GM46736 for the developments of effective Hamiltonian methods for enzymatic processes and grant no. RC1-GM091445 for the development of a quantal force field for biomolecular simulations), by the National Science Foundation (grant no. CHE09-57162 for the development of mixed molecular orbital and valence bond theory and grant no. CHE09-56776 for incorporating quantum mechanics in the treatment of complex reactive dynamical systems).

The neglect of diatomic differential overlap (NDDO)^{48}^{–}^{49} approximation is adopted in most modern semiempirical quantum mechanical methods such as MNDO,^{50} AM1,^{51} PM3,^{52} and RM1.^{53} In this method, electronic integrals involving the charge density
${\chi}_{\mu}^{A}{\chi}_{\nu}^{B}$ are assumed to be zero when the basis functions
${\chi}_{\mu}^{A}$ and
${\chi}_{\nu}^{B}$ are located on different atoms, *A* ≠ *B*, with the exception that the two-center one-electron nuclear attraction integrals are retained due to their importance for chemical bonding. Here we introduce a less severe approximation called the neglect of interfragment differential overlap (NIDO). In this approximation, we retain all electronic integrals within the same monomer, but electronic integrals involving differential overlap between basis functions belonging to different fragments are neglected. Applying the NIDO approximation greatly simplifies eq 15, as explained in the rest of this appendix. Below, we present the explicit expression for each term in eqs 16–18 by applying the NIDO approximation.

First, the block diagonal terms of the un-projected matrix are reduced to

$${F}_{\mu \nu}^{aa}=\left[{H}_{\mu \nu}^{aa}+\sum _{\lambda \sigma}^{{m}_{a}}{D}_{\lambda \sigma}^{aa}\{({\chi}_{\mu}^{a}{\chi}_{\nu}^{a}\mid {\chi}_{\lambda}^{a}{\chi}_{\sigma}^{a})-\frac{1}{2}({\chi}_{\mu}^{a}{\chi}_{\lambda}^{a}\mid {\chi}_{\nu}^{a}{\chi}_{\sigma}^{a})\}\right]+\left[\sum _{b\ne a}\sum _{\lambda \sigma}^{{m}_{b}}{D}_{\lambda \sigma}^{bb}({\chi}_{\mu}^{a}{\chi}_{\nu}^{a}\mid {\chi}_{\lambda}^{b}{\chi}_{\sigma}^{b})\right]$$

(A1)

Note that the “exchange” terms from orbitals on other fragments,
${D}_{\lambda \sigma}^{ab}({\chi}_{\mu}^{a}{\chi}_{\lambda}^{b}\mid {\chi}_{\nu}^{a}{\chi}_{\sigma}^{c})$, are all zero because of the NIDO approximation. Thus, the un-projected Fock matrix element for orbitals on fragment *a* can be written as:

$${F}_{\mu \nu}^{aa}={({F}_{\mu \nu}^{aa})}^{o}+<{\chi}_{\mu}^{a}\mid \sum _{b\ne a}{\mathbf{V}}_{b}^{a}\mid {\chi}_{\nu}^{a}>$$

(A2)

Eq A2 includes contributions only from the Coulomb integrals between different fragments, which is similar to a combined quantum mechanical and molecular mechanical (QM/MM) method. In the original X-Pol potential,^{2}^{–}^{3} in which the Hartree-product wave function is used, the rest of the density projection terms in eq 15 are not needed.

The *V*_{2} term of eq 16 simplifies under the NIDO approximation to

$$\begin{array}{l}{({V}_{2})}_{\mu \nu}^{aa}=\sum _{b\ne a}^{k}\sum _{\beta}^{{m}_{b}}\left[{H}_{\mu \beta}^{ab}+\sum _{c,d}^{k}\sum _{\lambda \sigma}^{{m}_{c},{m}_{d}}{D}_{\lambda \sigma}^{cd}\{({\chi}_{\mu}^{a}{\chi}_{\beta}^{b}\mid {\chi}_{\lambda}^{c}{\chi}_{\sigma}^{d})-\frac{1}{2}({\chi}_{\mu}^{a}{\chi}_{\lambda}^{c}\mid {\chi}_{\beta}^{b}{\chi}_{\sigma}^{d})\right]{P}_{\beta \nu}^{ba}\\ =-\frac{1}{2}\sum _{b\ne a}^{k}\sum _{\beta}^{{m}_{b}}\left[\sum _{\lambda}^{{m}_{a}}\sum _{\sigma}^{{m}_{b}}{D}_{\lambda \sigma}^{ab}({\chi}_{\mu}^{a}{\chi}_{\lambda}^{a}\mid {\chi}_{\beta}^{b}{\chi}_{\sigma}^{b})\right]{P}_{\beta \nu}^{ba}\end{array}$$

(A3)

in which, only the exchange terms having *a* = *c* and *b* = *d* remain to be non-zero.

The *V*_{3} term becomes

$$\begin{array}{l}{({V}_{3})}_{\mu \nu}^{aa}=\sum _{b\ne a}^{k}\sum _{\beta}^{{m}_{b}}{({P}_{\beta \nu}^{ba})}^{T}\lfloor {H}_{\beta \nu}^{ba}+\sum _{c,d}^{k}\sum _{\lambda \sigma}^{{m}_{c},{m}_{d}}{D}_{\lambda \sigma}^{cd}\{({\chi}_{\beta}^{b}{\chi}_{\nu}^{a}\mid {\chi}_{\lambda}^{c}{\chi}_{\sigma}^{d})-\frac{1}{2}({\chi}_{\beta}^{b}{\chi}_{\lambda}^{c}\mid {\chi}_{\nu}^{a}{\chi}_{\sigma}^{d})\rfloor \\ =\sum _{b\ne a}^{k}\left[-\frac{1}{2}\sum _{\beta}^{{m}_{b}}{({P}_{\beta \nu}^{ba})}^{T}\sum _{\lambda}^{{m}_{b}}\sum _{\sigma}^{{m}_{a}}{D}_{\lambda \sigma}^{ba}({\chi}_{\beta}^{b}{\chi}_{\lambda}^{b}\mid {\chi}_{\nu}^{a}{\chi}_{\sigma}^{a})\right]\end{array}$$

(A4)

Similarly, the non-zero terms are in the exchange part for orbitals that are on fragments such that *b* = *c* and *a* = *d* in the double summation over fragments in parentheses.

Finally, the last term of of eq 15 becomes

$$\begin{array}{l}{({V}_{4})}_{\mu \nu}^{aa}=\sum _{d\ne a}^{k}\sum _{\delta}^{{m}_{d}}\left\{\sum _{b\ne a}^{k}\sum _{\beta}^{{m}_{b}}{\left({P}_{\beta \mu}^{ba}\right)}^{T}\left[{H}_{\beta \delta}^{bd}+\sum _{c,f}^{k}\sum _{\lambda \sigma}^{{m}_{c},{m}_{f}}{D}_{\lambda \sigma}^{cf}\{({\chi}_{\beta}^{b}{\chi}_{\delta}^{d}\mid {\chi}_{\lambda}^{c}{\chi}_{\sigma}^{f})-\frac{1}{2}({\chi}_{\beta}^{b}{\chi}_{\lambda}^{c}\mid {\chi}_{\delta}^{d}{\chi}_{\sigma}^{f})\}\right]\right\}{P}_{\delta \nu}^{da}\\ =\sum _{b\ne a}^{k}\sum _{\beta}^{{m}_{b}}\sum _{\delta}^{{m}_{b}}{({P}_{\beta \nu}^{ba})}^{T}\left[{H}_{\beta \delta}^{bd}+\sum _{c}^{k}\sum _{\lambda \sigma}^{{m}_{c},{m}_{c}}{D}_{\lambda \sigma}^{cc}({\chi}_{\beta}^{b}{\chi}_{\delta}^{b}\mid {\chi}_{\lambda}^{c}{\chi}_{\sigma}^{c})\right]{P}_{\delta \nu}^{ba}-\frac{1}{2}\sum _{d\ne a}^{k}\sum _{\delta}^{{m}_{d}}\sum _{b\ne a}^{k}\sum _{\beta}^{{m}_{b}}{\left({P}_{\beta \mu}^{ba}\right)}^{T}\left[\sum _{\lambda}^{{m}_{b}}\sum _{\sigma}^{{m}_{d}}{D}_{\lambda \sigma}^{bd}({\chi}_{\beta}^{b}{\chi}_{\lambda}^{b}\mid {\chi}_{\delta}^{d}{\chi}_{\sigma}^{d})\right]{P}_{\delta \nu}^{da}\end{array}$$

(A5)

We presented the NIDO approximation here for a general case of X-Pol-X, that is the monomers can be treated by ab initio Hartree-Fock, by post-Hartree-Fock correlation methods, or by density functional theory. However, if desired, the monomers could also be treated by semiempirical molecular theory employing the NDDO approximation within the monomers. This procedure has been implemented into the program CHARMM for defining diabatic states in the effective Hamiltonian mixed molecular orbital and valence bond (EH-MOVB) method.^{25}

Supporting Information Available: Absolute energies and energy components in hartrees for trimer calculations. This material is available free of charge via the Internet at http://pubs.acs.org.

1. Xie W, Gao J. J Chem Theory Comput. 2007;3:1890. [PMC free article] [PubMed]

2. Gao J. J Phys Chem B. 1997;101:657.

3. Gao J. J Chem Phys. 1998;109:2346.

4. Wierzchowski SJ, Kofke DA, Gao J. J Chem Phys. 2003;119:7365.

5. Xie W, Orozco M, Truhlar DG, Gao J. J Chem Theory Comput. 2009;5:459. [PMC free article] [PubMed]

6. Song L, Han J, Lin YL, Xie W, Gao J. J Phys Chem A. 2009;113:11656. [PMC free article] [PubMed]

7. Xie W, Song L, Truhlar DG, Gao J. J Phys Chem B. 2008;112:14124. [PMC free article] [PubMed]

8. Xie W, Song L, Truhlar DG, Gao J. J Chem Phys. 2008;128:234108/1. [PubMed]

9. Cembran A, Song L, Mo Y, Gao J. J Chem Theory Comput. 2009;5:2702. [PMC free article] [PubMed]

10. Stoll H, Preuss H. Theor Chem Acc. 1977;46:12.

11. Stoll H, Wagenblast G, Preuss H. Theor Chim Acta. 1980;57:169.

12. Yang W. Phys Rev Lett. 1991;66:1438. [PubMed]

13. Mo Y, Peyerimhoff SD. J Chem Phys. 1998;109:1687.

14. Mo Y, Zhang Y, Gao J. J Am Chem Soc. 1999;121:5737.

15. Mo Y, Gao J, Peyerimhoff SD. J Chem Phys. 2000;112:5530.

16. Mo Y, Gao J. J Comput Chem. 2000;21:1458.

17. Gao J, Mo Y. Prog Theor Chem Phys. 2000;5:247.

18. Mo Y, Gao J. J Phys Chem A. 2001;105:6530.

19. Gao J, Garcia-Viloca M, Poulsen TD, Mo Y. Adv Phys Org Chem. 2003;38:161.

20. Brauer CS, Craddock MB, Kilian J, Grumstrup EM, Orilall MC, Mo Y, Gao J, Leopold KR. J Phys Chem A. 2006;110:10025. [PubMed]

21. Mo Y, Gao J. Acc Chem Res. 2007;40:113. [PubMed]

22. Mo Y. J Chem Phys. 2007;126:224104. [PubMed]

23. Song L, Gao J. J Phys Chem A. 2008;112:12925. [PMC free article] [PubMed]

24. Valero R, Song L, Gao J, Truhlar DG. J Chem Theory Comput. 2009;5:1. [PMC free article] [PubMed]

25. Cembran A, Payaka A, Lin YL, Xie W, Song L, Mo Y, Gao J. J Chem Theory Comput. 2010 ASAP. [PMC free article] [PubMed]

26. Hunt WJ, Hay PJ, Goddard WAI. J Chem Phys. 1972;57:749.

27. Mo Y, Gao J. J Phys Chem A. 2000;104:3012.

28. King HF, Staton RE, Kim H, Wyatt RE, Parr RG. J Chem Phys. 1967;47:1936.

29. Gallup GA. Int J Quantum Chem. 1972;6:899.

30. Shukla A, Dolg M, Stoll H, Fulde P. Chem Phys Lett. 1996;262:213.

31. Raimondi M, Gianinetti E. J Phys Chem. 1988;92:899.

32. Gianinetti E, Raimondi M, Tornaghi E. Int J Quantum Chem. 1996;60:157.

33. Khaliullin RZ, Head-Gordon M, Bell AT. J Chem Phys. 2006;124:204105/1. [PubMed]

34. Kitaura K, Ikeo E, Asada T, Nakano T, Uebayasi M. Chem Phys Lett. 1999;313:701.

35. Kitaura K, Sawai T, Asada T, Nakano T, Uebayasi M. Chem Phys Lett. 1999;312:319.

36. Fedorov DG, Kitaura K. J Phys Chem A. 2007;111:6904. [PubMed]

37. Leverentz HR, Gao J, Truhlar D. G Theor Chem Acc. 2010 submitted.

38. Nagata T, Takahashi O, Saito K, Iwata S. J Chem Phys. 2001;115:3553.

39. McWeeny R. Proc R Soc Lond A. 1959;253:242.

40. Su PF, Li H. J Chem Phys. 2009:131. [PubMed]

41. Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, Jensen JH, Koseki S, Matsunaga N, Nguyen KA, Su SJ, Windus TL, Dupuis M, Montgomery JS. J Comput Chem. 1993;14:1347.

42. Song L, Mo Y, Zhang Q, Wu W. J Comput Chem. 2005;26:514. [PubMed]

43. Schutz M, Burgi T, Leutwyler S, Burgi HB. J Chem Phys. 1993;99:5228.

44. Kitaura K, Morokuma K. Int J Quantum Chem. 1976;10:325.

45. Stevens WJ, Fink WH. Chem Phys Lett. 1987;139:15.

46. Chen W, Gordon MS. J Phys Chem. 1996;100:14316.

47. Mo Y, Gao J. J Phys Chem B. 2006;110:2976. [PubMed]

48. Pople JA, Santry DP, Segal GA. J Chem Phys. 1965;43:S129.

49. Pople JA, Segal GA. J Chem Phys. 1965;43:S136.

50. Dewar MJS, Thiel W. J Am Chem Soc. 1977;99:4899.

51. Dewar MJS, Zoebisch EG, Healy EF, Stewart JJP. J Am Chem Soc. 1985;107:3902.

52. Stewart JJP. J Comp Chem. 1989;10:209.

53. Rocha GB, Freire RO, Simas AM, Stewart JJP. J Comput Chem. 2006;27:1101. [PubMed]