Christophe J M Chipot
email: christophe.chipot@edam.uhp-nancy.fr

François Dehez
email: francois.dehez@edam.uhp-nancy.fr

Fabien Archambault

Ignacio Soteras Gutiérrez

Javier Luque

Franc ¸ois Dehez

F Javier Luque

An ab initio strategy for handling induction phenomena in metal ion complexes

Keywords: Polarization phenomena, Ab inito calculations, Macromolecular force field, Symmetry-adapted perturbation theory, Metal ion complexes

come   L'archive ouverte pluridisciplinaire

Introduction

Over the past years, the development of polarizable potential energy functions targeted at numerical simulations of molecular systems has focused tremendous effort. As access to massively parallel computer architectures opened new perspectives for the investigation of larger assemblies of atoms over time scales compatible with the experimentally observed phenomena, it was also realized that more elaborate force fields were needed for an accurate description of intermolecular interactions. Although nonpolarizable macromolecular potential energy functions, e.g. AMBER [1], CHARMM [2], GROMOS [3], or OPLS-AA [4],

have proven to behave reasonably well for biologically relevant molecular systems in which through-space polarization effects can be accounted for in an average sense, the pairwise additive approximation remains an intrinsic limitation for the study of biophysical processes where the induction effects at play cannot be neglected without proper justification. Among these processes, ions permeating a membrane channel have been shown to perturb the charge distribution of those residues lining the pore through which they diffuse, thereby modulating the interplay of the permeant and its surroundings. Similarly, ions embedded in metalloproteins polarize their environment, altering accordingly the charge distribution of the chelating agents [5], which suggests that induction effects may have a pronounced influence on the binding hierarchy of these ions.

A pioneering strategy devised for modeling explicitly through-space induction effects in numerical simulations consists in parameterizing the induction forces in terms of atomic quantities [6]. From a conceptual perspective, this solution assumes that the electron density response can be partitioned into regions of the Cartesian space corresponding to atoms. It further assumes a truncation of the multipole expansion and a selection of leading terms in the classical expression of the forces exerted between atomic distributions. Much effort has been invested in the task of partitioning the electron density response into distributed polarizabilities [7][8][9][10][11][12][13]. Transformation of cumbersome, large ensembles of components into finite sets of physically sound, optimally partitioned electric properties [14] (OPEP) constitutes a promising alternative to conventional distributed polarizability analysis. In a nutshell, the strategy embodied in OPEP relies upon a least-squares fit of atomic polarizabilities to the reference quantum-mechanical induction energy mapped around the molecule. Combining isotropic dipole polarizabilities with a zeroth-order charge-flow term between vicinal atoms [15], the anisotropy of induction effects can be essentially recovered, thereby obviating the requirement of an explicit inclusion of anisotropic contributions [16]. An ideal route for the design of polarizable force fields would consist in plugging directly into an existing pairwise additive potential energy function the polarizability parameters determined quantum mechanically. A number of obstacles preclude, however, the direct application of this strategy. First, the electrostatic term of nonpolarizable force fields markedly exaggerates the polarity of the constituent molecules, thereby compensating in an average sense for missing induction effects. In the context of an explicit treatment of polarization phenomena, atomic charges ought to be either scaled down [17], or consistently optimized to be representative of the true gas phase [18]. Furthermore, polarizability parameters are expected to depend inherently on the nature of the environment and, hence, should be adapted accordingly [19].

Here, description of both electrostatic and polarization terms of the potential energy function is consistent with a low-pressure gaseous state. In principle, these contributions should be tailored to reflect the characteristics of the surroundings.

Perhaps most critically, potential energy functions are very complex constructs, the building blocks of which are intimately related. Minute modifications of their parameters can perturb significantly the delicate, almost precarious balance between the different contributions. In addition to the scaling, or the novel derivation of atomic charges, van der Waals parameters ought to be crafted accordingly to mirror the gas phase electrostatics and the explicit inclusion of induction effects. Comparing the relative performances of potential energy functions to which polarizability parameters have been incorporated, yet without appropriate adjustment of the other terms, would obviously be unfair to the original pairwise additive force field.

Besides, it should be acknowledged that macromolecular, multi-purpose potential energy functions have not been designed for gas phase calculations. To illustrate the critical role played by polarization phenomena in intermolecular interactions, the complexes formed by a calcium ion and three possible electron-donor chelating agents, namely water, formaldehyde and formate, will be considered (see Figure 1). These test cases represent a paradigm for reproducing accurately the hierarchy of binding affinities in less rudimentary, possibly biologically relevant molecular systems -e.g. metalloproteins, in which a misdescription of the non-covalent interactions involving the ion may be liable for altering irreversibly the three-dimensional structure of the catalytic site and, hence, its function.

In the following section, the theoretical background is outlined, together with the computational details for the determination of atomic charges and polarizabilities. Relying on the functional utilized in macromolecular force fields like CHARMM [2], de novo parametrization of the Lennard-Jones potential is carried out based on symmetry-adapted perturbation theory [20,[START_REF] Bukowski | An ab initio program for manybody symmetry-adapted perturbation theory calculations of intermolecular interaction energies[END_REF] (SAPT) calculations, from which the exact individual contributions to the interaction energy can be readily inferred. Performance of the resulting ab initio polarizable force field are subsequently appraised in terms of reproducing the target quantum mechanical interaction energies for the three ionic complexes. Last, concluding remarks are drawn, highlighting the advantages and shortcomings of the proposed parametrization scheme.

Methods

Theoretical background

Interaction energies, ∆U MM tot , determined from classical polarizable force fields can be expressed as a sum of individual contributions:

∆U MM tot = ∆U ele + ∆U ind + ∆U damp + ∆U vdW (1)
where ∆U ele , ∆U ind , ∆U damp and ∆U vdW stand for the electrostatic, the induction, the damping and the van der Waals energy terms, respectively. SAPT calculations provide a rigorous partitioning of the total interaction energy, ∆U SAPT tot into physically meaningful components: Here, ∆U exch , ∆U disp , ∆U exch-ind and ∆U exch-disp denote the exchange, dispersion, exchange-induction and exchange-dispersion contributions to the interaction energy. δHF embraces induction and exchangeinduction higher order terms of the perturbative expansion.

∆U SAPT tot = ∆U ele + ∆U ind + ∆U exch + ∆U disp + ∆U exch-ind + ∆U exch-disp + δHF (2)
Armed with accurate models of net atomic charges and implicitly interacting distributed polarizabilities [16], the strategy pursued in the present work consists in optimizing the van der Waals part of the classical potential energy function in the light of SAPT calculations [START_REF] Dehez | [END_REF]. This parametrization scheme implies a number of underlying assumptions requiring particular attention, chief among which the ability of a simplistic Coulomb sum truncated at the monopole level to reproduce faithfully the quantum-mechanical electrostatic contribution ∆U ele . Brute incorporation of induction phenomena in molecular mechanical force fields raises additional concerns due to the physically unrealistic forces that can thrust the polarizing charge towards the polarizable center. In principle, the classical and the quantum-mechanical polarization terms are comparable, provided that spatial mapping of the induction energy by the model of atomic polarizabilities has the required precision. The SAPT expansion involves, however, additional contributions that ought to be incorporated in the classical description of ∆U tot . Such is the case of the exchangeinduction term, ∆U exch-ind , which is usually modeled by damping the interaction of the electric field with the polarizable sites, using a surrogate empirical function, ∆U damp . Whereas the latter contribution does not compare directly with its SAPT exchange-induction counterpart, it can be contended that the sum of ∆U ind and ∆U damp , in the classical description, can be related to the sum of ∆U ind and ∆U exch-ind , at the quantum-mechanical level [START_REF] Dehez | [END_REF].

Dissecting the classical van der Waals interaction, ∆U vdW , into meaningful contributions is unfortunately less straightforward. The questionable form of the Lennard-Jones potential utilized in a variety of macromolecular force fields like AMBER [1], CHARMM [2], GROMOS [3], or OPLS-AA [4], precludes any physical interpretation in terms of true repulsion and dispersion energies. In practice, Lennard-Jones potentials encompass contributions that cannot be captured in either the electrostatic or the polarization terms of the classical description, thus, going beyond the components that can be recovered rigorously from an SAPT expansion -viz. the exact dispersion, ∆U disp , exchange, ∆U exch , and exchange-dispersion, ∆U exch-disp , contributions. Noteworthy, a number of classical approaches have been devised to represent the latter sum of contributions, going beyond a mere Lennard-Jones term, among which the SIBFA [23],

the ORIENT [START_REF] Stone | Orient: a program for studying interactions between molecules[END_REF] and the EFP [25] potential. Classical equivalence of the SAPT δHF term is even less obvious, as most of the contributions embraced in the latter do not necessarily appear in the parametrization of the classical polarizability models. For instance, hyper-polarizability effects are generally neglected by computing the target quantum-mechanical induction energy sufficiently far from the nuclei.

Computational details

Models of net atomic charges were derived using OPEP [14] from the electrostatic potential computed quantum mechanically at the MP2/6-311++G(d, p) level of approximation, including inner-shell orbitals [26,27] (see Table I). The grids over which the electrostatic potential created by water, formaldehyde and formate was mapped contained 2,106, 3,079 and 3,056 points, respectively. Models of implicitly interacting distributed polarizabilities were also derived with OPEP, based on maps of induction energies generated using a finite-perturbation approach [28] and the formalism of topologically partitioned electric properties (TPEP) [10,12]. The induction-energy maps of water, formaldehyde and formate were obtained using grids of 3,905, 19,152 and 19,099 points, respectively, at the MP2 level of theory with the Sadlej basis set [29],

which has proven to be well suited for the determination of polarizability parameters [28]. The calcium cation was assumed to be nonpolarizable. For all chelating agents considered in this work, preliminary optimization of the molecular geometries was performed at the MP2/6-311++G(2d, 2p) level, employing the Gaussian98 suite of programs [START_REF] Frisch | GAUSSIAN 98 Revision A.7[END_REF].

The inherent anisotropy of induction phenomena can be recovered without the explicit incorporation of anisotropic dipole components in the distributed polarizability models. As can be seen in Table II, zero-order charge-flow polarizabilities and first-order isotropic dipole polarizabilities offer a reasonable compromise between accuracy and tractability.

Electrostatic and polarization models were associated to the CHARMM [2] macromolecular force field for the computation of the classical potential energy surfaces that delineate the interaction of the calcium ion with the three electron donors along their C 2 axis. Cognizant of the imbalance entailed by the inclusion of induction phenomena in the potential energy function, de novo optimization of the constituent Lennard-Jones parameters was carried out using a Levenberg-Marquardt nonlinear fitting algorithm [START_REF] Marquardt | [END_REF]. Assuming that (i) the pure electrostatic component of equation ( 2 RMSD (in 10 -3 a.u.) between the reference quantum-mechanical electrostatic potential mapped on a grid of points and that regenerated from the model of distributed multipoles; ∆ε (in percents) is the corresponding mean error

12" potential, based on the remaining terms of the perturbative expansion. Last, for comparison purposes, the quantum-mechanical potential energy surface characterizing the three complexes was determined at the MP2/6-311++G(d, p) level of theory, varying the intermolecular distance by increments of 0.1 Å and evaluating the basis set superposition error (BSSE) [32] at each step.

Results and discussion

Metalloproteins represent approximately one third of all known proteins that serve catalytic, regulatory, or structural roles [33]. Yet, modeling of these entities remains a difficult undertaking due to the intrinsic non- ent proteins, they were able to preserve the integrity of the active site over multi-nanosecond molecular dynamics simulations in a condensed phase. This static approximation of both induction phenomena and charge transfer cannot, however, mimic with a proper accuracy processes involving a modification of the surroundings -e.g. insertion, substitution and permutation of a given ligand. Sophisticated potentials including an explicit treatment of these nonadditive interactions have been successfully derived for a variety of metalloproteins [START_REF] Gresh | [END_REF]34,[39][40][41][42].

(a) (b) Here, a similar route is followed, attempting to plug simple, yet accurate polarizable models into a CHARMM-like potential energy function. The chosen prototypical interactions of a calcium cation with different chelating agents -viz. water, formaldehyde and formate are commonly found at the active site of calcium metalloproteins, where water or amino acids like aspartate or glutamate form recurrent motifs in the first coordination shell of the cation [43]. The main focus of the present work is on the ability of the various models to reproduce correctly the hierarchy of interaction energies for a calcium ion and a number of its usual partners at the catalytic site of a metalloprotein. The performance of the derived ab mechanical gas phase interaction energy characterizing the binding of a calcium ion by formaldehyde. In the case of formate, CHARMM performs quite nicely and predicts a binding energy minimum close to the MP2/6-31++G(d, p) reference, albeit slightly shifted towards shorter separations. The MP2 interaction energy estimated for the calcium-formaldehyde complex is about 7.8 kcal/mol lower than the value computed for the calcium-water pair, as highlighted in Table III of the CHARMM force field glaringly exaggerates the binding energy for all three ion metal complexes e.g. by up to 11.2 kcal/mol in the case of the water-calcium dimer (see Figure 2 and Table III). The van der Waals part of the classical CHARMM potential energy function has been optimized for a given set of net atomic charges, assumed to reflect the interaction of the parameterized chemical moieties with an aqueous environment. Alteration of the electrostatic contribution of the nonpolarizable force field creates an imbalance in the construct, which ought to be corrected by de novo optimization of its Lennard-Jones parameters. This rationalizes the disagreement between the profiles of Figures The electrostatic contributions estimated from an SAPT expansion and from a model of simple atomic charges coincide both at the equilibrium geometry and at longer separation distances for all three metal ion complexes. At the quantum-mechanical minimum of the binding energy reported in Table IV, the predicted molecular-mechanical electrostatic energy of -51.8 kcal/mol matches perfectly the SAPT value of -51.7 kcal/mol for the water-calcium ion dimer. In contrast, a deviation of less than 7 % is observed between the contributions from the polarizable models and the SAPT reference for formaldehyde and formate (see Table III). At shorter intermolecular separations, the agreement deteriorates, the SAPT electrostatic contributions becoming more attractive than the classical model as the intermolecular distance decreases.

This deviation is believed to be rooted in charge penetration effects that a classical formulation of the charge interactions clearly cannot capture [44,[START_REF] Stone | The theory of intermolecular forces[END_REF] without an appropriate treatment. These effects arise when two interacting chemical moieties are brought at a close enough distance, such that their respective charge den- expected to undergo a greater attraction for the electron density of the other compound. One possible road to account for electron-cloud penetration effects consists in introducing an appropriate correction function into the classical formulation of the interaction of two or more multipole expansions [START_REF] Piquemal | [END_REF]47].

Considering the polarization term of the interaction energy in Figures 2, 3 and4, the pure induction contribution, ∆U ind , to the SAPT expansion is considerably stronger than the value derived from distributed polarizability models. At the equilibrium geometry, ∆U ind inferred from the SAPT expansion is almost an order of magnitude larger than the classical estimate, and even larger at shorter distances. As has been shown previously [48], ∆U ind in SAPT calculations is excessively large on account of Coulomb singularities in the interaction operator. This shortcoming has already been highlighted by Claverie in related perturbational approaches [START_REF] Claverie | Intermolecular Interactions: From Diatomics to Biopolymers[END_REF]. These singularities also yield very large positive estimates of the exchangeinduction term, ∆U exch-ind . Assignment of a physical meaning to these individual contributions, therefore, constitutes a daunting task. As put forth by Misquitta and Stone [START_REF] Misquitta | [END_REF], the vocabulary SAPT induction energy denotes the sum of the two contributions, ∆U ind and ∆U exch-ind . For polar molecules, nonlinear or many-body, induction and exchange-induction contributions are often large and cannot be neglected.

These third-and higher-order induction and exchange-induction terms are referred to as δHF in the SAPT expansion. They do not appear explicitly in the molecular-mechanical derivation of the interaction energy.

Misquitta and Stone [START_REF] Misquitta | [END_REF] have, however, shown that the higher-order energies could be partially recovered through an appropriate damping function associated to a polarizable model. Admittedly, the damping function adopted in this work [51,[START_REF]IM-SQRT damping function proposed by Jensen et al in reference 51[END_REF] has not been designed to account for these higher-order effects. It may be contended, however, that a comparison of the induction energies computed from the damped classical polarizable model with the total induction of the SAPT calculation supplemented by δHF is legitimate.

As can be observed in Figure 2 and Table III, an almost perfect agreement for the water-calcium ion complex at any separation distances is attained. In the case of formaldehyde depicted in Figure 3, the induction energy computed from the ab initio polarizable force field departs from the reference SAPT expansion for distances shorter than 3.6 Å. This discrepancy probably stems from the limited ability of the distributed polarizability model to reproduce correctly the molecular polarizability tensor. Moreover, the damping function developed by Jensen et al. [51] has been derived from a series of molecules only remotely related to the formaldehyde-calcium ion complex, thereby suggesting the need to design a better appropriate damping function, or, at the very least, an improved parametrization thereof. In the challenging example of formate associated to the calcium ion shown in Figure 4, the induction energy obtained from the damped polarizable model matches satisfactorily the estimate of the SAPT calculation. Yet, the classical curve is slightly shifted, compared to the SAPT profile. As can be seen in Table III, this discrepancy on the order of a few kcal/mol is in line with the 9.2 kcal/mol deviation observed between the total energy of interaction computed at the SAPT and at the MP2/6-311++G(d,p) level. It is also noteworthy that in this study the calcium ion has been considered non-polarizable. However weak, the polarizability of the cation should be included if an extremely precise interaction energy is sought.

The van der Waals contribution is expressed in the SAPT expansion as the sum of exchange, dispersion, and exchange-dispersion terms whereas it is described by means of a Lennard-Jones potential in the present ab initio polarizable force field. For metal ion complexes, brute-force comparison of the two descriptions delineated in see Figures 2,3 and 4 reveals that the ab initio polarizable force field strongly exaggerates the repulsive part of the SAPT reference curve. This illustrates the incapacity of a rudimentary"6-12" potential to reproduce simultaneously the repulsive and the attractive contributions to the SAPT van der Waals energy, thus, calling for the development of alternate mathematical forms of this potential. Such a strategy has been successfully endeavored by Piquemal et al. in the case of divalent cations interacting with the sophisticated polarizable AMOEBA water model [START_REF] Piquemal | [END_REF]. The use of the function devised by Halgren [54] to represent the van der Waals interactions, together with an ad hoc parametrization of the Thole damping function [55] proved to reproduce accurately the overall interaction energy in these metal ion complexes.

Conclusion

In the present contribution, an ab initio strategy based on accurate quantum mechanical calculations is introduced for handling induction phenomena in metal ion complexes. The proposed polarizable force field combines potential derived net atomic charges with atom-centered charge-flow plus isotropic dipole polarizabilities determined numerically from induction energy maps. Investigation of the complexes formed by a calcium ion and three different electron-donor chelating agents, viz. water, formaldehyde and formate, suggests, however, that concomitant reproduction of the electrostatic potential and the polarization energy of the isolated chemical species, does not guarantee an accurate description of intermolecular interactions upon association. The results reported here emphasize that electrostatic and induction models ought to be distinguished from interaction models. Precise description of the target electric properties embodied in low standard deviations with respect to the reference observables for the isolated species does not necessarily imply that the classical electrostatic and induction contributions to the interaction energy match their quantum-mechanical counterparts inferred from an SAPT expansion. As exemplified here, the naive choice of turning to a functional common to many non-polarizable macromolecular force fields, with minimal alteration of the latter for the sake of adaptability, suffers from severe shortcomings, as physically meaningful terms are missing in the proposed description where induction effects are introduced explicitly.

Chief amongst these missing contributions are the short-range penetration of electron clouds, anticipated to counterbalance the repulsive behavior witnessed at small distances separating the ion from the chelating agent -e.g. formate probably constitutes the most striking illustration of this shortcoming. Not too surprisingly, penetration effects are also present in the induction term. Mathematical surrogate functions aimed at damping the interaction of the electric field with the participating polarizable sites generally prove to be suboptimal at short separations -i.e. where the consequences of electron-cloud penetration are expected to be pronounced. Furthermore, ad hoc damping functions are seldom designed to account for higher-order induction and exchange-induction contributions, which are embraced in the so-called δHF term of the perturbative expansion. It still remains that, however perfectible, damping functions represent a reasonable compromise to recover physically meaningful induction component of the classical interaction energy. Most importantly, the form of the energy profile that encompasses the SAPT exchange, dispersion and exchange-dispersion contributions cannot be modeled in a satisfactory fashion by means of the mathematically questionable "6-12" form of a simplistic Lennard-Jones potential. If the attractive part of the curve is favored in the course of the nonlinear fitting procedure, its repulsive part will then inevitably be poorly described, and vice versa. By and large, a strategy that consists in assigning to the van der Waals

Anisotropy of the interaction

To illustrate the ability of the proposed strategy to model the anisotropy of the interaction, binding of the divalent calcium cation to the formate ion was investigated, considering the two additional approaches depicted in Figure 1. Ion pairing was examined as a function of the distance separating Ca 2+ from the central carbon atom of the anion (i) following a monodentate chelation pattern, and, (ii) in the direction orthogonal to the plane of the latter.

Figure 1 :

 1 Figure 1: Interaction of a divalent calcium ion with (a) water, (b) formaldehyde, and (c) a formate ion. Approach of the cation towards the three ligands is considered along the C 2 axis of the latter.

Figure 2 :

 2 Figure 2: Interaction of a divalent calcium ion with water. (a) Comparison of the binding energies determined from BSSE-corrected MP2/6-311++G(d, p) calculations (dark solid line), the classical, nonpolarizable CHARMM force field (black lines), and the latter supplemented by a model of distributed polarizabilities (light lines), with and without de novo optimization of the participating Lennard-Jones parameters. The electrostatic contribution to the binding energy is depicted as dashed lines. (b) Components of the interaction energy, determined from an SAPT expansion using the 6-311++G(d, p) basis set (dark solid line) and the polarizable potential energy function (light solid line). The undamped induction contribution corresponds to the pure induction term of equations (1) and (2). At the quantum-mechanical level, the damped induction contribution consists of a sum of induction and exchange-induction terms. At the molecular-mechanical level, it stands for the pure induction component corrected by a damping function. The van der Waals contribution encompasses at the quantum-mechanical level the dispersion, the exchange and the exchange-dispersion terms of the SAPT expansion. In the classical description, it coincides with the Lennard-Jones potential.

Figure 3 :Figure 4 :

 34 Figure 3: Interaction of a divalent calcium ion with formaldehyde. (a) Comparison of the binding energies determined from BSSE-corrected MP2/6-311++G(d, p) calculations (dark solid line), the classical, nonpolarizable CHARMM force field (black lines), and the latter supplemented by a model of distributed polarizabilities (light lines), with and without de novo optimization of the participating Lennard-Jones parameters. The electrostatic contribution to the binding energy is depicted as dashed lines. (b) Components of the interaction energy, determined from an SAPT expansion using the 6-311++G(d, p) basis set (dark solid line) and the polarizable potential energy function (light solid line). The undamped induction contribution corresponds to the pure induction term of equations (1) and(2). At the quantum-mechanical level, the damped induction contribution consists of a sum of induction and exchange-induction terms. At the molecular-mechanical level, it stands for the pure induction component corrected by a damping function. The van der Waals contribution encompasses at the quantum-mechanical level the dispersion, the exchange and the exchange-dispersion terms of the SAPT expansion. In the classical description, it coincides with the Lennard-Jones potential.

  term of the macromolecular force field the challenging role of compensating the flaws of all other individual contributions has necessarily stringent limitations. The ability of the proposed strategy to hierarchize correctly at the minimum of the interaction energy the binding affinity of a divalent calcium ion towards three distinct electron donors, nevertheless, provides very encouraging results, underlining the strength of ab initio polarizable energy functions for handling accurately induction phenomena. The shortcomings brought to light in the course of this investigation provide new guidelines for the improvement of the interaction models, paving the way for the development of increasingly accurate polarizable force fields and opening new vistas for the realistic description of the catalytic site of metalloproteins with the guarantee to preserve its structure over long time scales.

F

Figure 1 :Figure 2 :

 12 Figure 1: Interaction of a divalent calcium cation with a formate anion. Approach of the former towards the latter is considered (a) following a monodentate chelation motif, and (b) along the axis perpendicular to the plane of the anion.

Table I :

 I Models of distributed multipoles and regenerated molecular multipole moments of water, formaldehyde and formate at the MP2/6-311++G(d, p)//MP2/6-311++G(2d, 2p) level of approximation.

		distributed	molecular
		multipoles	multipoles	MP2/6-311++G(d, p)
	Water	Q O 00 Q H 00 o F RMSD	-0.782 0.391 1.029	Q 10 Q 20	-0.869 -0.219	-0.857 -0.274
		r ∆ε	51.119
	Formaldehyde	Q C 00 P Q H 00 e 0.429 -0.005 Q O 00 -0.419 e RMSD 0.432 r ∆ε 7.568	Q 10 Q 20 Q 30	-0.949 -0.284 -1.332	-0.946 -0.488 0.086
	Formate	Q C 00 Q H 00 Q O 00 RMSD ∆ε	0.864 R -0.210 e -0.827 v Q 20 Q 40 -17.507 4.456 0.468 i w e 0.283	4.563 -13.313
					O n
					l y

) and its classical counterpart coincide, and (ii) the damped induction contribution, i.e. ∆U ind + ∆U damp , reproduces appropriately the sum of the two SAPT terms, ∆U ind + ∆U exch-ind , the van der Waals component of ∆U MM tot was determined in the form of a "6-

Table II :

 II Models of distributed polarizabilities and regenerated molecular polarizabilities of water, formaldehyde and formate at the MP2/Sadlej//MP2/6-311++G(2d, 2p) level of approximation.

			Molecular Physics Molecular Physics
		distributed	molecular
		polarizabilities	polarizabilities	MP2/Sadlej
	Water	F o α OH 00,00 α OO 1κ,1κ F o	-0.808 8.180	α 10,10 α 11c,11c 11.483 10.177 α 11s,11s 8.180	9.751 10.063 9.542
		r RMSD r	0.127
	Formaldehyde	P e e r ∆ε 7.000 α CH 00,00 -1.984 P α CO 00,00 -2.070 α CC 1κ,1κ 8.133 α OO 1κ,1κ 4.087 e e r	α 10,10 α 11c,11c 12.220 27.934 α 11s,11s 24.680	22.419 12.844 17.966
	Formate	RMSD ∆ε α CH 00,00 α CO 00,00 α CC 1κ,1κ α OO 1κ,1κ RMSD	R e v i e w 0.136 14.228 -2.201 α 10,10 -2.104 α 1c,1c 10.459 α 1c,1c 7.134 0.265 R e v i e w	38.900 24.732 44.321	36.551 24.684 42.080
		∆ε	2.700	O n O n
					l y l y

Page 8 of 35

URL: http://mc.manuscriptcentral.com/tandf/tmph RMSD (in 10 -3 a.u.) between the reference quantum-mechanical induction energy mapped on a grid of points and that regenerated from the model of distributed polarizabilities; ∆ε (in percents) is the corresponding mean error additive character of the metal-ligand interactions involved. Polarization and charge transfer are known, indeed, to represent significant contributions to the overall interaction energy

[34,35]

. Standard pairwise additive force fields used to investigate biological systems cannot, however, capture in any effective way the underlying physical complexity of these interactions. So far, common solutions imagined to address this problem consisted either in considering a bonded model of the active site -i.e. all metal-ligand interactions are treated as non-dissociable bonds, or to derive ad hoc electrostatic atomic charge models of the metal-ligand moieties

[36]

. Recently Dal Peraro et al. have introduced an elegant method to model metalloproteins

[5]

. Deriving point charges from the Atoms in Molecules theory of Bader

[START_REF] Bader | Atoms in Molecules -A Quantum Theory[END_REF]

for differ-Page 9 of 35 URL: http://mc.manuscriptcentral.com/tandf/tmph

Table III :

 III Comparison of the contributions to the binding energies a of the complexes formed by Ca 2+ with water, formaldehyde and formate, determined from an MP2/6-311++G(d, p) SAPT expansion and a polarizable potential energy function.

				Molecular Physics
		∆U ele		∆U b ind		∆U c vdW			∆U tot
		SAPT	MM	SAPT	MM	SAPT MM	δHF d	MP2 e	SAPT	MM f
	Ca 2+ • • •H 2 O Ca 2+ • • •H 2 CO	F -51.7 o -51.8 -48.4 -45.1 r	-22.1 -23.0 -37.5 -31.3	19.8 22.7 (3.6) 24.2 15.3 (4.8)	-0.8 -1.6	-53.4 -61.2	-54.7 -63.3	-52.2 (-64.6) -61.1 (-71.6)
	Ca 2+ • • •HCOO --308.7 -290.0 P -62.2 -58.3 r e e	65.8 62.8 (1.4)	-0.5	-297.8 -305.6	-285.6 (-301.6)
				R e v i e w
						O n
						. On the contrary, the classical, nonpolarizable force l y
	field favored the binding of the water to the calcium cation by about 2-3 kcal/mol. As a consequence, the
	CHARMM model allows carbonyl groups to exchange freely with water, when, in reality, quantum mechan-
	ics indicates that this process is highly improbable. It is, therefore, becoming quite clear that CHARMM

would hardly preserve the integrity of a catalytic site involving a calcium ion interacting concomitantly with an aspartate or a glutamate residue and water.

The ab initio polarizable potential energy function that utilizes the standard Lennard-Jones parameters Page 12 of 35 URL: http://mc.manuscriptcentral.com/tandf/tmph a All quantities are expressed in kcal/mol, with respect to the quantum-mechanical energy minima of Table IV. b The SAPT value includes the pure induction and the exchange-induction contributions. The MM model consists of the pure induction term, supplemented by a damping correction. c The SAPT value corresponds to the sum of the exchange, the dispersion and the exchange-dispersion contributions. The MM value is simply the Lennard-Jones component of the CHARMM force field. d This contribution encompasses the third and higher order induction and exchange-induction terms of the SAPT expansion. It is clearly absent in the MM description. e BSSE-corrected MP2/6-311++G(d, p) interaction energies. f The values in bold were determined after de novo optimization of the Lennard-Jones parameters. The values in parentheses correspond to estimates obtained with the standard Lennard-Jones parameters of the CHARMM force field.

Table IV :

 IV viz. -53.4 kcal/mol. An equally remarkable agreement is obtained when the calcium ion is bound to formaldehyde. In the case of the strongest interaction featuring formate, the observed deviation of less than 5% remains satisfactory. It is also noteworthy that the SAPT binding energies overshoot the MP2/6-311++G(d, p) estimates by up to 7.8 kcal/mol in the particular instance of the formate-calcium ion complex. This deviation might be ascribed to the frozen core approximation used in MP2 calculations. The repulsive part of the total binding energy curves computed at the molecular-mechanical level is, however, largely overestimated in view of the quantum-mechanical target. A closer look at the components of the interaction energy in Figures2, 3 and 4and in TableIIIilluminates the origin of this discrepancy. Intermolecular separations a in the complexes formed by Ca 2+ with water, formaldehyde and formate, determined from MP2/6-311++G(d, p) potential energy surfaces, the classical CHARMM force field and the ab initio polarizable force field.

			Molecular Physics
	F o
	r
	P
	e e r	d X•••Y
	Ca 2+ • • •H 2 O Ca 2+ • • •H 2 CO Ca 2+ • • •HCOO -	QM 2.3 2.2 2.5	CHARMM e R 2.3 v 2.1 i w e 2.3	CHARMM + polz. 2.4 2.2 2.6
				O n
					l y
					2, 3 and 4 computed at the
	MP2/ 6-311++G(d, p) level of approximation and using the CHARMM force field supplemented by atomic
	polarizabilities.

Design of a new set of Lennard-Jones parameters based on the SAPT expansion yields binding energy profiles in qualitative agreement with the quantum-mechanical references. As can be seen in Table

III

, the binding energy estimated for the water-calcium ion complex from the ab initio polarizable force field viz. -52.2 kcal/mol, compares nicely with the quantum-mechanical energy minima of Table

IV

-URL: http://mc.manuscriptcentral.com/tandf/tmph a All quantities are expressed in Å.

  sities overlap. The nuclei on one molecule are no longer shielded by the electron density, and, hence, is

	Molecular Physics	Page 14 of 35
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

URL: http://mc.manuscriptcentral.com/tandf/tmph

URL: http://mc.manuscriptcentral.com/tandf/tmph Molecular Physics

Acknowledgments

János G. Ángyán is gratefully acknowledged for stimulating discussions. The authors are indebted to the Centro Nacional de Supercomputación, Barcelona, Spain, for provision of generous amounts of computational time. This work was supported by a joint CNRS/UIUC collaborative grant.

Supplementary Material:

An ab initio strategy for handling induction phenomena in metal ion complexes


El Mendoub

El Bahloul

Verlaine Metz

Lpmd Charpentier

E B El Mendoub

J.-F Wax

I Charpentier

N Jakse

O n l y Integral Equation Study of the Square-Well Fluid for Varying Attraction Range

Keywords: Integral equations, square-well fluid, liquid-vapour coexistence, colloids

In the present work, the properties of the square-well ‡uid (SW) for varying attraction range are investigated in the framework of the integral equation method. The pair-correlation functions are calculated for a wide range of the potential attraction including the colloidal regime, and the liquid-vapour phase diagram as well as the critical parameters are determined for attraction ranges spanning from the van der Waals limit down to the metastable region. An overall good agreement is found compared to available simulation data as well as several theoretical approaches.

Complex systems such as colloidal suspensions or protein solutions have attracted much attention from the scienti…c community. An important issue in this …eld is to understand the role played by the range of attraction in controlling of the thermodynamics and the phase behavior in such systems driven by complex e¤ective interactions 1 . It has been argued that the details of the interactions might have little in ‡uence, and a generalized law of corresponding states seems to govern the system. The size of the particles mapped by an equivalent hard-sphere (HS) diameter and the second virial coe¢ cient, containing a simple energetic parameter, act as the scaling variables [2][3][4] . This reinforces the possibility to describe such kind of systems with very simpli…ed models of interactions.

Interestingly, the square-well (SW) potential whose the potential energy as a function of the distance of separation r reads

u(r) = 8 > > > < > > > : 1; r ; "; < r ; 0 r > ; (1)
which is the simplest model that exhibits a liquid-vapour equilibrium and obeys a law of corresponding states. Its three parameters are the size of the hard spheres, , the strength ", and the width of the potential-well. Obviously, the HS potential is recovered when " = 0 , otherwise the range of attraction can be controlled by varying such that various kind of ‡uids can be treated.

The in…nitely short-range limit ( ! 1, and " ! 1) corresponds to the sticky HS model 5 and short ranges of the potential ( < 1:25) are taken as models of colloidal suspensions 6 and protein solutions 7 . Intermediate ranges (1:5 < 2) intend to mimic simple ‡uids, and longer ranges correspond to very shallow potentials, for instance with ! 1, corresponding to the van der Waals mean …eld theory, a limit which is numerically reached 8 for = 3. Therefore, due to the simplicity of its mathematical expression and its applicability to a wide class of ‡uids, the SW model represents a useful benchmark for theoretical developments. Most of the works reported in the literature were devoted to the intermediate-ranged potentials for which the model behaves in a qualitative manner as a normal ‡uid. Some interest has regained recently since it has been shown that short ranges are particularly well-suited to the description of colloids 6,9,10 .

A large amount of work was devoted to the understanding of the structural and thermodynamic properties including the liquid-vapour phase coexistence of the SW ‡uids. Numerical simulations, based on molecular dynamics and Monte-Carlo (MC) methods 8,[11][12][13][14][15] , were mostly used to study the SW potential for the intermediate-ranges of the potential. Only a few of them were devoted to the short-ranged domain 2, 15,16 . Analytical theories for the structure 17,18 , thermodynamics 19 , and liquid-vapour coexistence 20 were developed on the basis of either the van der Waals theory or HS perturbation method or the renormalization group (RG) approach 21 for various ranges of the potential. Besides, semi-analytical methods, based on the statistical physics of ‡uids such as the integral equations (IE), which consists in solving numerically the Ornstein-Zernike (OZ) equation together with an approximate closure relation have been thoroughly used in the seventies and eighties 22 . In this framework, the hierarchical reference theory (HRT) 23 and the self-consistent OZ approximation (SCOZA) 24 are successful to deal with the critical behaviour of the SW ‡uids for various ranges of attraction.

Of particular interest is the integral equation proposed by Sarkisov 25 to treat the structure and the thermodynamics of ‡uids, which has proven to achieve a good degree of thermodynamic consistency 26 . It has been previously applied to a wide class of Lennard-Jones type potentials 27 and also to discontinuous potentials such as the HS 28 and the attractive Yukawa potentials 29 .

The aim of the present work is to investigate the properties of the SW ‡uids covering attraction part of the potential from the short-range to the van der Waals limit on the same footing. For this purpose, the have carried out calculations using the Sarkisov 25,29 integral equation method coupled with e¢ cient numerical tools, namely, the tangent linear method to get exact derivatives and an adaptative technique 27 . We show that within the present IE approach: (i) a correct description of the pair-correlation functions can be obtained for a wide range of the potential attraction including the colloidal regime and (ii) the liquid-vapour phase diagram as well as the critical parameters can be calculated for the whole domain of for which coexistence exists.

The layout of the paper is the following. Section 2 brie ‡y describes the integral equation method, the thermodynamic properties as applied to the SW model, as well as the numerical techniques to construct the liquid-vapour phase diagrams. The results for the structural properties and the liquid vapour coexistence are presented and discussed in Section 3. In Section 4 we summarize and conclude.

F

h(r) = c(r) + Z h(r 0 )c( r r 0 )dr 0 = c(r) + (r): (2)
The integral equation method consists in solving Eq. (2) together with a closure relation whose formal expression for the pair-correlation function g(r) (= h(r) + 1) is

g(r) = exp [ u(r)] exp[ (r) + B(r)] ; (3)
= 1=k B T being the inverse temperature and k B Boltzmann's constant. The bottleneck of the method is that Eq. (3) depends on the bridge function, B(r), which represents an in…nite sum of elementary bridge diagrams [START_REF] Martynov | Fundamental Theory of Fluids. Methods of Distribution Functions[END_REF] . Following Sarkisov's scheme 25 the bridge function is written as

B(r) = [1 + 2 (r)] 1=2 1 (r): (4)
The function B(r) is a functional in terms of the indirect correlation function and also of the interaction potential. In this closure relation, a renormalization concept of the indirect correlation function has been introduced by considering (r) rather than (r) in B(r). Such an idea was implemented by Lee et al. 32 for the Lennard-Jones with (r) = (r) u 2 (r), where u 2 (r) is the attractive part of the pair potential according to the Weeks, Chandler and Andersen (WCA) splitting [START_REF] Weeks | [END_REF] . For the hard-sphere (HS) potential, an e¢ cient renormalization was found, namely (r) = (r) + f (r)=2, f (r) (= exp( u(r)) 1) being the Mayer function 33 . In a recent work 29 , we have extended this renormalization to the hard-core Yukawa potential as

(r) = 8 < : (r) + f ( + )=2; r < ; (r) + f (r)=2; r ; (5)
that we adopt also here for the SW model of interaction. The renormalized indirect correlation function is then a discontinuous function. However this poses no conceptual problem since the numerical solution is obtained through the LVM method in the Fourier space where the correlation functions are continuous. Speci…c forms of the thermodynamic properties have to be used for the SW potential. The excess energy E is expressed as

E N = 2 " Z g(r)r 2 dr (6)
The compressibility factor Z is determined from the virial equation of states such as :

Z = P = 1 + 2 3 3 g( + ) 3 g ( ) + g ( ) + ; ( 7
)
where P is the pressure. The inverse isothermal compressibility, 1= T , is obtained from a density derivative of the pressure

@P @ = 1 + 4 3 3 g( + ) + 2 3 2 3 @g( + ) @ 4 ( ) 3 3 g ( ) + g ( ) + 2 @g(( ) + ) @ @g(( ) ) @ (8)
Finally, the excess chemical potential is calculated using the approximate form proposed by Kyseliov and Martynov 34

ex = 4 Z 1 0 (r) + B(r) h(r) + 1 2 h(r) (r) + 4 3 B(r) r 2 dr: (9)
The construction of the liquid-vapour coexistence curves is based on conditions of thermal, mechanical and chemical equilibria. The thermodynamic points of the liquid (l) and vapour (v)

phases belonging to the binodal line satisfy the relations:

T l = T v ; P l = P v ; l = v : (10)
The chemical potentials l and v are obtained from Eq.( 9) by

= ex + k B T ln( 3 ).
The domain of mechanical stability is limited by the spinodal line, which can be obtained from the state points corresponding to the extrema of the virial pressure along each isotherm. These are the points in the metastable region where the inverse isothermal compressibility given by Eq. ( 8) becomes zero. As argued by Sarkisov 26 , when approaching the instable region, correlation functions including g(r) itself display a diverging correlation radius that a¤ects the thermodynamic calculations due to long-range truncation e¤ects. However, the inverse compressibility calculated from the virial route and given by Eq. ( 8) does not depend explicitely on the asymptotic form of the correlation functions, as only the values of g(r) at the discontinuities are involved, and therefore it is not a¤ected by these long-range e¤ects. Consequently, the right-hand side of Eq. ( 8) can be calculated as soon as a valid solution of the IE exists, as it is not subject to long range truncation e¤ects when approaching the unstable region. We thus de…ne the spinodal line 35 as a zero value of Eq. ( 8). Our experience is that the quantity 1= T can be negative when a solution of IE is found in the instable region being therefore beyond the spinodal line.

C. Numerical techniques

The accurate calculation of the structure, thermodynamics and the liquid-vapour coexistence curve relies on e¢ cient numerical tools that are presented brie ‡y in this paragraph.

First, a numerical solution of the pair-correlation function within the integral equation is obtained through the algorithm of Labik et al. 36 . It combines the Newton-Raphson and the successive substitution (standard Picard scheme) methods guaranteeing a good stability and a fast solution.

In a …rst stage, successive approximations of the indirect correlation function (r; ; T ) for a given state point ( ; T ) are generated through the following iterative process:

i+1 (r; ; T ) = M[ i (r; ; T )]; (11)
where M is a nonlinear functional derived from the OZ equation in the reciprocal space as

~ (q) = c2 (q) 1 c(q) : (12)
The direct correlation function c(q) is expressed as a function of (r) in the direct space by means of the closure relation, which in the present work is given by Eqs. ( 4) and (5). The second stage corresponds to the Newton-Raphson part, in which the functional M is linearized with respect to as:

M[ ] ' M[ 0 ] + : @M @ [ 0 ]; (13)
about initial guess 0 (r) assumed to be su¢ ciently close to the solution (r), i.e. (r) = 0 (r) + (r) with (r) small (the interested reader is referred to 36 for the speci…c expression of M as well as technical details). The NR method, which consists in estimating the unknown di¤erence at each step of the iterative process given by Eq. ( 11) converges quadratically. Finally, a solution at state point ( ; T ) is found when vanishes. From a technical point of view, the correlation functions are discretized in the r space on a grid of 4096 points with a mesh of 0:0125 , which implies a spatial extension R max = 51:2 , su¢ cient for an accurate numerical integration and a fast Fourier transform. Secondly, the determination of the thermodynamic quantities of interest have expressions containing not only the pair-correlation function g(r) but also its thermodynamic derivatives either with respect to density or temperature that have to be evaluated. This is the case for the isothermal compressibility given by Eq. (8). For this purpose the tangent linear method [START_REF]Automatic Di¤ erentiation of Algorithms: Theory, Implementation and Applications[END_REF] gives us the possibility to compute locally all the thermodynamic functions in an "exact" manner for each state point (T; ), i.e. without having recourse to a …nite di¤erence method. The tangent linear technique is applied directly to the code of Labik et al., using an automatic di¤erentiation tool, namely, the software Odyssee 39 (we refer the interested reader to previous papers for a more detailed description [START_REF] Charpentier | [END_REF] ).

Finally, as the thermodynamic properties are calculated locally at any state point (T; ), an adaptative technique enables us to perform an e¢ cient and accurate systematic mapping of the thermodynamic space, in order to construct easily the binodal and spinodal lines as close to the critical point as possible. The detailed description can be found in a reference paper 27 , and will not be repeated here. It is worth emphasizing that the adaptative technique used here drives the whole process and makes calls to the di¤erentiated code of Labik et al. at each state point needed for the mapping. The latter is started with an initial guess 0 (r; ; T ) and the solution is computed at the nearby state point ( ; T + T ) where T is small enough to get a small di¤erence , which is one of the keys for a fast resolution of the Labik et al. algorithm for successive state points.

Technically, a good compromise between accuracy and speed for the two-dimensional grid is to use an initial temperature step T = 0:01k B =" and density step = 0:01 3 .

III. RESULTS AND DISCUSSION

The central quantities entering the calculation of the thermodynamic functions needed to construct the liquid-vapour phase diagram are the pair-correlation function g(r) and its temperature and density derivatives. A comparison with simulation data o¤ers the possibility of probing the IE directly at the microscopic level. In Figures 1 and2 2 show a reasonable agreement at low density. At high densities, the IE method underestimates substantially the MC simulation data regarding the height of the …rst peak of g(r) and overestimates the oscillations for distances just beyond r = . It is worth mentioning however that the concordance is still reasonable for = 1:2.

As can be seen in Eqs. ( 6) to (8), the values of g( + ), g( ) and g( + ) drive the thermodynamics of the SW ‡uid. The disagreement with the simulation data limits to = 1:05 and = 1:1 at high density, as already seen for the pair-correlation function. Finally, the results found in Table I demonstrate the reliability of the present IE scheme, i.e. its ability to describe the structural properties of the SW ‡uid over a wide domain of thermodynamic states as well as of interaction ranges , except for the two lowest values of .

A stringent test for the reliability of the present integral equation approach is to determine the liquid-vapour phase diagram and to compare the predicted results to available simulation data. We recall that the accuracy of the phase diagram lies on the quality of the calculated thermodynamic functions which in turn depend essentially on the values of g( + ), g( ) and g( + ) studied above, as well as their thermodynamic derivatives calculated in an exact manner through tangent linear technique. On the one hand, the pressure and the chemical potential given respectively by Eqs.

(7) and ( 9) are needed for the binodal line and, on the other hand, the inverse compressibility (Eq. ( 8)), that controls the mechanical stability of the system, yields the spinodal line. For the determination of the critical parameters, we follow the technique described in our preceeding work 27 and we refer the reader to it for further details. In this method, the error bars are controled by …xing the minimum value of the temperature step and that of the density step of the adaptative grid. In the vicinity of the critical point these are 0:001k B =" and 0:001 3 , respectively.

Let us …rst examine the results for the long attraction range (2 < 3). To our best knowledge, no simulation data exist for > 3, above the van der Waals limit 8 . Figure 3 shows the IE curves which are compared to those of the SCOZA method as well as to the Gibbs ensemble Monte-Carlo (GEMC) simulation data of Patel et al. 13 . For = 3 (upper curve) the GEMC simulation data of Orkoulas and Panagiotopoulos 12 as well as the molecular dynamics results of Benavides et al. 8 have also been included. For the four values of covering this range, the IE method yields a binodal line in very good agreement with the simulation data. It appears that the IE and SCOZA 24 predictions are of comparable accuracy for the high density branch, while the SCOZA seem to be slightly better for the low density branch for = 2:8 and 3. The critical parameters are gathered in Table II along with the simulation data when available and with those of theoretical approaches, namely SCOZA, HRT 23 and RG 21 . These have been inferred within the IE scheme from the position of the maximum of the spinodal line. The obtained values are close to SCOZA, HRT and RG ones, especially, the critical temperature.

Figure 4 shows the liquid-vapour phase diagram for the intermediate attraction range (1:5

2). An overall good agreement is found with respect to GEMC data of Patel et al. 13 for the three values of considered. Interesting enough, for the high density (liquid) branch, it appears that the IE gives better results than the SCOZA for = 1:5 and = 1:75. For the coexistence curves, calculated for the shorter range of attraction and shown in Figure 5, a di¤erent situation emerges. In the case of = 1:25, both the IE and SCOZA, which give similar results, fail to reproduce the high density branch when compared to the simulation data of Patel et al. 13 and del Rio et al. 14 . For = 1:375, this is also the case, nevertheless the IE approach gives better predictions. It is worth mentioning however that, for 1:25 2 the critical parameters inferred from the IE (see Table II) are close to the SCOZA ones, and are also in reasonable agreement with the MC data. It appears that, away from the critical region, the IE is more reliable than the SCOZA for the liquid branch of the coexistence curve.

As can be seen on Figures 3,4, and 5, the temperature range of the liquid-vapour coexistence is reduced with decreasing values of . The question of the disappearance of the liquid-vapour coexistence with the narrowing of the potential well has been studied for various type of potentials 41 .

For the SW potential, it is worth mentioning that Tejero and co-workers 42 shown that such a "metastable" liquid-vapour coexistence could be determined by simulation 16 down to = 1:15. In the same spirit, we have tested the ability of the IE method to determine the liquid-vapour coexistence curves for values of below 1:25, and Figure 5 shows the binodal lines obtained for = 1:20 and = 1:175. The critical temperatures and densities predicted by our IE approach are lower that those obtained from HRT. We mention that below = 1:175, no binodal line could be calculated and the value of 1:175 can by considered as the limit for the present integral equation scheme, which is very close to that found from the computer simulation 16 .

IV. CONCLUSION

To summarize the integral equation method is used to study the properties of the squarewell ‡uid for varying attraction range. The calculations are carried out with the self-consistent closure relation proposed by Sarkisov 25 . We have shown that the pair-correlation functions can be accurately described for a wide range of the potential attraction including the colloidal regime by comparison with simulation data. For the smallest values of (1:05 and 1:1) correct results are obtained only at low densities and high temperatures.

The liquid-vapour phase diagram as well as the critical parameters are obtained for 1:175 3, using the adaptative technique associated to the tangent linear method to get exact thermodynamic derivatives of the pair-correlation function. In particular, for long attraction range

(2 < 3) the integral equations predicts accurately the binodal line with respect to the simulation data. Our results are close to those obtained from the SCOZA method 24 , including the critical parameters. In the intermediate range of attraction (1:5 2), a reliable description of the liquid-vapour coexistence curves and critical parameters is found, while for shorter values

(1:25 < 1:5), namely down to the stability limit 16 the results are only in qualitative agreement with the simulation data, especially for the liquid branch of the binodal. Nevertheless, it appears that for the intermediate and short ranges, the latter is generally better described by the IE scheme than the SCOZA. Results for the liquid-vapour phase diagram is also obtained in the region with < 1:25, and a limit of = 1:175 is found, close to that obtained by state-of the-art simulation technique 16 . The dash-dotted lines represent the curve m 3 = ( l 3 + v 3 )=2, where l 3 and v 3 are the densities of the liquid and vapour branches of binodal line. The thin lines are the SCOZA results 24 , while the open circles, up triangles and squares correspond respectively to the simulation data of Patel et al. 13 , Orkoulas and Panagiotopoulos 12 , and Benavides et al. 8 . Table I: Comparison of the contact values g( + ), g( + ) and g( ) of the pair-correlation function obtained from our IE method with Monte Carlo data of 18 for several values of the attraction range , density 3 and temperatures k B T =".

Table II: Coordinates of the critical point inferred from the integral equations method for di¤erent values of , along with the MC simulation data of Bolhuis et al. 11 , Reyes and Orea 15 (a)

and Orkoulas and Panagiotopoulos 12 (b), as well as theoretical results from the SCOZA 24 , HRT 23 and renormalization group (RG) 21 approaches. The error bars are 0:001k B =" for the temperature and 0:001 3 for the density (see text).

B T c =" c 3 k B T c =" c 3 k B T c =" c 3 k B T c =" c 3 k B T c =" c 3

  Let us consider a ‡uid at a temperature T , and number density in which the particles interact through a pair-potential u(r). According to the OZ equation, the total correlation function, h(r), between two particles is de…ned as the sum of the direct correlation function, c(r), and the indirect correlation function, (r):

Page 5 of 28 URL

 28

  , the calculated g(r) are respectively plotted for values of in the domain [1:3; 2:0], belonging to the intermediate attraction range, and [1:05; 1:2] corresponding to the colloidal regime. Our results are compared with existing Monte-Carlo simulations in the canonical ensemble (N V T ) by Largo et al.18 . In the intermediate attraction range, the curves reveal an excellent agreement with the MC simulation data at both low and high densities, whatever the values of . The height and shape of the …rst peak is well reproduced, as well as the subsequent oscillations for r > . In the short range of attraction range, the curves of Figure

  have shown that the simple van der Waals theory predicts a vanishing of the coexistence around = 1:25. More recently, Liu et al.16 using GEMC simulation in conjunction with an aggregation volume biased technique have pointed out that the triple point disappears for = 1:25 and a sublimation line takes place for all temperatures. The liquid-vapour coexistence, which can nevertheless be constructed, become metastable and hidden under it. Such …ndings could explain the di¢ culties encountered by both the IE and SCOZA methods to reproduce the binodal line for this value of .Liu et al. have

  1 M. H. J. Hagen, E. J. Meijer, G. C. A. M. Mooij, D. Frenkel, and H. N. W. Lekkerkerker, Nature (London) 365, 425 (1993).

10 Page 10 of 28 URLFigure 1 :

 10281 Figure 1: Pair-correlation function g(r) in the intermediate attraction range for density 3 = 0:2 (a), and 3 = 0:8 (b). In both panels, the solid, dashed and dotted lines correspond respectively to the IE results for = 1:30 and k B T =" = 1:0, = 1:50 and k B T =" = 1:0, and = 2:0, and k B T =" = 3:0. The symbols (circles, up triangles and diamonds respectively) are the corresponding MC simulation data 18 .

Figure 2 :

 2 Figure 2: Pair-correlation function g(r) in the short attraction range for density 3 = 0:2 (a), and 3 = 0:8 (b). In both panels, the solid, dashed and dotted lines correspond respectively to the IE results for = 1:05 and k B T =" = 0:5, = 1:10 and k B T =" = 0:5, and = 1:20 and k B T =" = 0:7. The symbols (circles, up triangles and diamonds respectively) are the corresponding MC simulation data 18 .

Figure 3 :

 3 Figure 3: Liquid-vapour phase diagram of a SW ‡uid for long attraction range, namely = 2:3, = 2:5, = 2:8, and = 3:0 from the bottom to the top. The thick lines and dashed lines correspond to the results of the IE method respectively for the binodal and spinodal lines.

Figure 4 :

 4 Figure 4: Liquid-vapour phase diagram of a SW ‡uid for intermediate attraction range, namely = 1:5, = 1:75, and = 2:0 from the bottom to the top. The thick lines and dashed lines correspond to the results of our IE method respectively for the binodal and spinodal lines. The dash-dotted lines represent the curve m 3= ( l 3 + v 3 )=2, where l 3 and v 3 are the densities of the liquid and vapour branches of binodal line. The thin lines are the SCOZA results24 , while the open circles correspond to the simulation data of Patel et al.13

Figure 5 :

 5 Figure 5: Liquid-vapour phase diagram of a SW ‡uid for intermediate and short attraction range, namely = 1:175, = 1:20, = 1:25 and = 1:375 from the bottom to the top. The thick lines and dashed lines correspond to the results of our IE method respectively for the binodal and spinodal lines. The dash-dotted lines represent the curve m 3 = ( l 3 + v 3 )=2, where l 3 and v 3 are the densities of the liquid and vapour branches of binodal line. The thin lines are the SCOZA results 24 , while the open circles and open squares correspond respectively to the simulation data of Patel et al. 13 and del Rio et al. 14 .

  Figure 1

  TableIgathers these values as well as the corresponding MC data of Largo et al.18 to provide a general and quantitative view of the performance of the present IE scheme. Di¤erent values of are considered in the short and intermediate attraction range, as a function of density for several isotherms. The general trends that emerge from the analysis of Table I con…rms the one inferred from the curves of g(r). The results are especially in good accordance for the intermediate range of . They remain good with the simulation data in the domain of corresponding to the short range of attraction at low and moderate densities, provided that the temperature is not too low.

Page 14 of 28 URL: http://mc.manuscriptcentral.com/tandf/tmph

	Molecular Physics
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

TABLE I

 I

		k B T ="	3	g( + )	g( )	g( + )
				IE MC	IE MC	IE MC
	1.2	0.7	0.1 4.093 4.126 4.002 4.030 0.965 0.968
			0.5 3.151 3.449 2.844 3.122 0.693 0.748
			0.7 3.308 3.281 2.579 2.810 0.708 0.675
		1.0	0.1 2.679 2.671 2.624 2.628 0.965 0.968
			0.5 2.773 2.799 2.346 2.415 0.890 0.888
	F o r 1.5	0.7 3.284 3.173 2.269 2.372 0.901 0.871 0.1 2.009 2.016 1.953 1.956 1.003 1.004 0.5 2.539 2.541 2.009 2.042 1.062 1.050 0.7 3.253 3.168 2.010 2.063 1.087 1.059
	P 1.5 1.5	0.1 1.955 1.952 1.836 1.832 0.941 0.941
		0.5 2.174 1.989 1.378 1.382 0.648 0.709 e 0.7 3.069 2.783 1.191 1.115 0.454 0.590 e r 2.0 0.1 1.661 1.661 1.563 1.563 0.944 0.945
			0.5 2.167 2.006 1.261 1.274 0.732 0.771 R 0.7 3.150 2.880 1.063 1.063 0.563 0.646
		3.0	e 0.1 1.448 1.44 1.350 1.351 0.966 0.968
			w (Continued) 0.5 2.154 2.046 1.160 1.171 0.813 0.839 v i e 0.7 3.169 2.998 0.974 0.982 0.651 0.704
					O n l
						y

TABLE I (

 I Continued)

	Page 17 of 28		Molecular Physics Molecular Physics		Page 16 of 28
	IE	k B T ="	3 MC	g( + )	g( ) SCOZA	g( + ) HRT	RG
	k			IE MC	IE MC	IE MC
	1.8	2.0	0.1 1.844 1.839 1.606 1.599 0.975 0.974
			0.5 2.293 2.199 1.201 1.175 0.691 0.714
			0.7 3.348 3.487 1.184 1.141 0.671 0.693
		3.0	0.1 1.458 1.449 1.318 1.320 0.944 0.942
			0.5 2.211 2.162 1.086 1.079 0.760 0.773
	F o r 5.0 F o r	0.7 3.271 3.392 1.067 1.043 0.741 0.747 0.1 1.299 1.297 1.172 1.169 0.958 0.956 0.5 2.171 2.152 1.005 1.005 0.815 0.822 0.7 3.216 3.318 0.982 0.970 0.794 0.794
	P 3.0 3.0 P	0.1 1.563 1.564 1.352 1.346 0.969 0.968
		e e r 0.5 2.362 2.424 1.213 1.199 0.865 0.861 e 0.7 3.398 3.669 1.331 1.322 0.964 0.949 5.0 0.1 1.320 1.325 1.178 1.178 0.960 0.960 e r
			R 0.5 2.151 2.300 1.114 1.109 0.912 0.910 0.7 3.290 3.471 1.228 1.225 1.015 1.005 R
				e e
				v i e v i e
						w w
						O n l O n l
							y y

URL: http://mc.manuscriptcentral.com/tandf/tmph URL: http://mc.manuscriptcentral.com/tandf/tmph