PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
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/ct100268p
PMCID: PMC2922695
NIHMSID: NIHMS223806

On the Interfragment Exchange in the X-Pol Method

Abstract

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.

1. Introduction

The explicit polarization (X-Pol) method,1 initially called molecular orbital derived potentials for liquids (MODEL), was introduced in 19972 for the treatment of macromolecular systems in statistical mechanical Monte Carlo and molecular dynamics simulations using electronic structure theory.35 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.23 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.23,68 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 Preuss1011 (whose emphasis was on the localization of orbitals on atomic centers with the neighboring bonding partners and lone-pairs) and Yang12 (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.13,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.13,56 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,1425 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,1314,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,3132 Mo et al.,1316,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) theory16,2223,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.3435 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 method3435 are essentially the same as that (see page 660 of ref 2) in X-Pol theory.23 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.

2. Method

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.23,68 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.

2.1. Block localization and the X-Pol wave function

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 ka basis functions and na electrons. Thus, the total number of primitive basis functions, K, and the total number of electrons, N, in the system are:

K=a=1MkaandN=a=1Mna
(1)

The molecular orbitals ϕja in each block a are written as linear combinations of the primitive basis functions { χμa; μ = 1, ···, ka} in that specific subspace:

ϕja=μ=1kacjμaχμ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:

ΨXPol=a=1MA^Φa
(3)

where  is an antisymmetrization operator, and Φa is a product of the occupied spin-orbitals in block a (eq 2):

Φa=ϕ1aϕ2aϕnaa
(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

Fa=Ha+i(Jia12Kia)+baVba
(5)

where Ha is the one-electron Hamiltonian matrix, Jia and Kia are the Coulomb and exchange integral matrices associate with molecular orbital i in block a, and the last summation term is the Coulomb interaction matrix due to the nuclei and wave functions of all other blocks in the system. This is a level comparable to a force field for macromolecular systems.12,6,8 The Fock matrix Fa depends on the instantaneous polarization of the wave functions of other blocks through Vba, which, in turn, depends on Fa. The electronic energy of the full system is obtained by the DSCF method. We have presented a nonvariational version of the Fock operator to illustrate the main ideas as simply as possible, which is sufficiently accurate for use in statistical mechanical Monte Carlo simulations without the need for computing analytical gradients,24 and this approach has been adopted in other applications.3436 A variational version that can be used to calculate analytic gradients, suitable for geometry optimization and molecular dynamics simulations,5 is presented elsewhere.78 The final SCF energy is called the BL term or E SCF.

We note that the electrostatic potential Vba, 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.12,56,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, εabLJ, which also account for long-range dispersion interactions.23 Thus, the total X-Pol energy of the system is

EXPol=ESCF+12abaεabLJ
(6)

2.2. X-Pol-X: Block-localized wave function that accounts for exchange interactions

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:1316,27

ΨXPolX=A^a=1MΦa
(7)

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

C=(C1000C2000Ck)
(8)

where Ca is an ka×ka matrix, with the orbital coefficients of eq 2 arranged as column vectors. As mentioned in the introduction, this is equivalent to constructing a block-localized wave function (BLW) using the block-localized molecular orbitals (BLMO),1315 and the computational procedure has been described by several authors in various contexts,9,1316,27,3233,38 dating back to the work of Stoll.11

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:

<ϕiaϕjb>={δijaa,a=bSijab,ab
(9)

where Sijab is the overlap integral between molecular orbitals i and j from blocks a and b, respectively, and is given by

Sijab=(cia)Rabcjb
(10)

where { cia; i = 1, ···, na, a = 1, ···, ka} are the columns of Ca in eq 8 (the component of these vectors are { ciμa; μ = 1, ···, ka}), where T denotes a transpose, and where Rab is the block of the overlap matrix R corresponding to basis functions in fragments a and b:

Rμνab=χμa|χνb
(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 { cia} for block Ca in eq 8. by

D=C(CRC)1C
(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 as9,16,27

Fμν=Hμν+λσKDλσ{(χμχνχλχσ)12(χμχλχνχσ)}
(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

EXPolX=μνKDμν{Hμν+Fμν}
(14)

In comparison with the X-Pol SCF energy ESCF (eq 6), the exchange energy included in the XPol- X expression (eq 14) originates from the orthogonalization constraint (CTRC)−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 (ka × ka) by transforming the original Fock matrix into a block-diagonal form using a projection operator.9,11 Thus, for block a,

Jaa=Faa+V2+V3+V4
(15)

where

V2=bakFab(cakDbcRca)
(16)

V3=bak(cakDbcRca)Fba
(17)

and

V4=dakbak(cakDbcRca)Fbd(cakDdcRca)
(18)

The upper indices in eqs 1518 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-Fock1316,27 and DFT calculations.9 In the latter case, the exchange potential in Hartree-Fock theory is replaced by the exchangecorrelation potential.9

3. Computational Details

All computations have been performed using a locally modified GAMESS program41 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-W3, in which each water molecule donates and accepts one hydrogen bond from one of the other two water monomers.43 All three pairwise interactions are stabilizing. For the second structure, we constructed a trimer complex structure that has a significantly Coulomb repulsive interaction as follows. We first optimized the dimer complex structure, and then we placed the third water molecule at the C2 image position of the donor water molecule about the bisection axis of the acceptor water molecule. This complex is denoted as the symmetric trimer, s-W3. Both structures are depicted in Figure 1.

Figure 1
Illustration of the minimum water trimer structure (c-W3) in (a) and a symmetric configuration (s-W3) 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-W3 and the same as constructed above for s-W3. The X-Pol-X wave function is the three-block localized system  ψa ψb ψc. The energy and energy components of the X-Pol-X wave function will be compared to those for the unrelaxed wave function A^ψa0ψb0ψc0, where ψw0 is the wave function of a monomer in the gas phase (unpolarized wave function) at the geometry in the trimer complex.

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

4. Results and Discussion

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 A^ψa0ψb0ψc0 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.

Table 1
Computed relative energies for the cyclic water trimer minimum structure (c-W3) and for a symmetric trimer geometry (s-W3) (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 ψa0ψb0ψc0, and row 3 is calculated using A^ψa0ψb0ψc0. 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-W3 and s-W3 complexes, whereas row 3 shows that the exchange-repulsion energies are 16.3 and 10.8 kcal/mol, respectively. The sum of Coulomb and exchange energies is often called the total unpolarized electrostatic energy due to the charge density of each monomer in the gas phase alone; this is given in row 4, and it amounts to −9.3 and −5.2 kcal/mol, respectively, for c-W3 and s-W3. The water trimer has been extensively studied, but they are not listed here since our goal is not to investigate the water clusters. In one energy decomposition calculation using the reduced variational space (RVS) method,45 which is analogous to the BLW-ED process, Chen and Gordon found that the total electrostatic interaction energy is −9.0 kcal/mol using the 6-31++G(d,p) basis set at the minimum configuration.46

Table 2
Additivity of dimer contributions to the Coulomb and exchange-repulsion interaction energies in the cyclic water trimer minimum structure (c-W3) and in a symmetric trimer geometry (s-W3) (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-W3 interaction energy from −9.3 to −12.5 kcal/mol and the s-W3 interaction energy from −5.2 to −6.7 kcal/mol. The differences of these quantities are listed in row 8, which shows a net X-Pol polarization energy of −3.2 and −1.5 kcal/mol, or 26% and 22%, for the two complexes, respectively. The electronic polarization effects in these two water trimer complexes are illustrated in Figure 2, which depict electron density shifting away for the hydrogen atoms that donate hydrogen bonds to the acceptor oxygen atoms. Concomitantly, electron density is attracted towards regions where oxygen atoms bind the proton donor. It is of interest to notice that in the s-W3 structure (Figure 2b), electron density is reduced from both hydrogen atoms in the W1 water molecule thanks to block-localization that disallows charge transfer between different water molecules, although they are not directly donating any hydrogen bonds. In the study of Chen and Gordon, polarization effects were found to contribute −2.7 kcal/mol with the slightly larger 6-31++G(d,p) basis set for c-W3.46

Figure 2
Electron density difference contour between the polarized and unpolarized systems, Δρ=ρ(A^ψaψbψc)ρo(ψaoψboψco), 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-W3) has three significant pairwise exchange repulsion interactions, while the other (s-W3) has only two. The finding of high additivity for exchange repulsion interactions suggests that it is sufficient to construct an X-Pol-X wave function that involves only pair-wise antisymmetrized blocks, which can greatly reduce the computational cost of the fully antisymmetrized system.

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.

Figure 3
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 Cs 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.

5. Conclusions

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.

Supplementary Material

1_si_001

Acknowledgments

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).

Appendix. X-Pol-X with the neglect of interfragment differential overlap approximation

The neglect of diatomic differential overlap (NDDO)4849 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 χμAχνB are assumed to be zero when the basis functions χμA and χνB are located on different atoms, AB, 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 1618 by applying the NIDO approximation.

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

Fμνaa=[Hμνaa+λσmaDλσaa{(χμaχνaχλaχσa)12(χμaχλaχνaχσa)}]+[baλσmbDλσbb(χμaχνaχλbχσb)]
(A1)

Note that the “exchange” terms from orbitals on other fragments, Dλσab(χμaχλbχνaχσ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μνaa=(Fμνaa)o+<χμabaVbaχν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,23 in which the Hartree-product wave function is used, the rest of the density projection terms in eq 15 are not needed.

The V2 term of eq 16 simplifies under the NIDO approximation to

(V2)μνaa=bakβmb[Hμβab+c,dkλσmc,mdDλσcd{(χμaχβbχλcχσd)12(χμaχλcχβbχσd)]Pβνba=12bakβmb[λmaσmbDλσab(χμaχλaχβbχσb)]Pβνba
(A3)

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

The V3 term becomes

(V3)μνaa=bakβmb(Pβνba)THβνba+c,dkλσmc,mdDλσcd{(χβbχνaχλcχσd)12(χβbχλcχνaχσd)=bak[12βmb(Pβνba)TλmbσmaDλσba(χβbχλbχνaχσa)]
(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

(V4)μνaa=dakδmd{bakβmb(Pβμba)T[Hβδbd+c,fkλσmc,mfDλσcf{(χβbχδdχλcχσf)12(χβbχλcχδdχσf)}]}Pδνda=bakβmbδmb(Pβνba)T[Hβδbd+ckλσmc,mcDλσcc(χβbχδbχλcχσc)]Pδνba12dakδmdbakβmb(Pβμba)T[λmbσmdDλσbd(χβbχλbχδdχσd)]Pδνda
(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

Footnotes

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.

References

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]