Long-range dispersion effects on the water/vapor interface simulated using the most common models
LLong-range dispersion effects on the water/vaporinterface simulated using the most common models
Marcello Sega and Christoph Dellago ∗ Faculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090, Vienna, Austria
E-mail: [email protected]
Abstract
The long-range contribution to dispersion forcesis known to have a major impact on the proper-ties of inhomogeneous fluids, and its correct treat-ment is increasingly recognized as being a neces-sary requirement to avoid cutoff-related artefacts.Although analytical corrections for quantities likethe surface tension are known, these can not takeinto account the structural changes induced by thelong-range contributions. Here, we analyze the in-terfacial properties of seven popular water models,comparing the results with the cut-off version ofthe dispersion potential. The differences in sur-face tension estimates are in all cases found to beless than 2 mN/m.
Introduction
Efficiency considerations have often been at theroot of the practice, customary in computer sim-ulations of molecular systems, to truncate disper-sion forces between pairs of particles separated bymore than a given cut-off. In contrast to the elec-trostatic force, whose contribution from far par-ticles can be comparable to that of close ones, dispersion forces are decaying much faster, pro-viding often a justification for the truncation. Theimportance of taking into account the full disper-sion forces, however, has been known since longin the field of crystallography, where it is essen-tial for the precise estimation of binding energies(it is no wonder that perhaps the most celebrated ∗ To whom correspondence should be addressed method for the calculation of the full electrostaticenergy, developed by Ewald, was devised pre-cisely to compute the energy of ionic crystals). Forhomogeneous, isotropic liquids, the analytical tailcorrections to the energy and (more importantly)pressure of dispersion forces are usually sufficientto remove artifacts due to truncation. Analyticalcorrections to the surface tension are also known, that apply to inhomogeneous systems in slab con-figuration.Explicit simulations of the liquid/vapor coexis-tence are usually performed at fixed volumes, andthe correction to the surface tension cannot beused as a feedback to change the density of thefluid (as in the case of constant pressure simula-tions). In’t Veld, Ismail and Grest showed thattaking into account the full dispersion forces hasa dramatic impact on the density profile of theLennard-Jones liquid/vapor interface. Similarly,Wennberg and coworkers demonstrated the impor-tance of the long-range part of dispersion forcesin membranes. These effects, although less pro-nounced in the case of pure water because thedominant contribution to the cohesion forces is ofelectrostatic nature, follow a clear trend, observ-able in the binodal line, of increasing the densityof the liquid phase. As we will show, the analyt-ical corrections tend to overestimate the surfacetension for temperatures T lower than (cid:39) K ,and to underestimate it for larger ones, althoughalways within about 2 mN/m.In the following sections we will first discussour methodological approach, including the sim-ulation details and the calculation of the analyticaltail contribution, and then we will present the re-1 a r X i v : . [ c ond - m a t . s o f t ] D ec igure 1: Snapshots of two equilibrium config-urations of TIP4P water (550 K, no long-rangedispersion forces). Upper panel: the system isclearly composed of one liquid phase and one va-por phase. Lower panel: the system separates intotwo separate liquid slabs.sults of our simulations on the binodal line and onthe surface tension of several, commonly used wa-ter models. Methods
We simulated seven popular rigid watermodels, including three-, four-, and five-point ones, namely, SPC, SPC/E, TIP3P, TIP4P, TIP4P/2005, TIP5P and, TIP5P-E. All simulations were performed using the GRO-MACS molecular dynamics simulation pack-age, version 5.1, and consisted of 1000 watermolecules in a slab configuration in a rectangularsimulation box of 5 × ×
15 nm , with periodicboundary conditions applied in all directions. Thesimulations were performed in the canonical en-semble by integrating the equation of motions us-ing the leapfrog algorithm (1 fs timestep), impos-ing the temperature equilibrium value by meansof a Nosé–Hoover thermostat (2 ps relax-ation time), and keeping the molecules rigid usingthe SETTLE algorithm. The electrostatic en-ergy, force, and pressure were calculated using thesmooth Particle Mesh Ewald method (sPME)with a real space mesh spacing of 0.15 nm, a 10 − relative accuracy of the potential energy at thereal space cutoff of 1.3 nm, and an interpolationscheme of order 4. In all simulated models the dis- persion forces are taken into account by using theLennard-Jones pair interaction. For each modeland each chosen temperature T (300, 350, 400,450, 500 and 550 K ) we performed one simulationby truncating dispersion forces at r c = . , with a 10 − relative accuracy at the real space cut-off of 1.3 nm(the same as in Ref. ). After an equilibration of2 ns, an equilibrium trajectory of 50 ns was gener-ated, saving configurations to disk at 1 ps intervals,and pressure every 0.1 ps. The density profileswere calculated as ρ ( z ) = (cid:104) ∑ i δ ( z − z i + z cm ) (cid:105) ,where the sum is extended over all molecules, z cm is the position of the center of mass of thesystem, and the angular brackets represent the en-semble average. Prior calculation of the profile,the liquid phase was determined in every frameusing a cluster search based on a distance crite-rion: two molecules that dist less than 0.35 nmare considered to be in the same cluster, and thelargest cluster in the system identifies the liquidphase. The liquid phase slab was then shifted bya suitable amount along the macroscopic interfacenormal, in order not to cross the box boundaries,thus removing the ambiguity in the definition ofthe center of mass of a periodic system.The analytical correction to the surface ten-sion can be obtained by calculating the integral γ tail = πεσ ( ρ L − ρ V ) (cid:90) ds (cid:90) ∞ r c dr ×× coth ( rs / d )( s − s ) r − . (1)The densities of the liquid ( ρ L ) and of the va-por ( ρ V ) phases as well as the interface thickness d can be estimated by performing a Marquardt-Levenberg least square fit of the sampled densityprofiles (see Fig. 2 for an example) to the function ρ ( z ) = ( ρ L + ρ V ) − ( ρ L − ρ V ) tanh [( z − z ) / d ] , (2)where z is the position along the interface normal,and z identifies the location of the middle of theinterface. In order to perform the integral, Eq. (1),we used the general purpose adaptive quadraturefor infinite intervals (QAGI) of the QUADPACKlibrary, as simpler integration schemes showed2 T = 300K T = 450K T = 550K ρ n m z / nm Figure 2: Density profile in the interfacial regionfor the TIP4P-2005 model, with (full lines) andwithout (dashed lines) long range Lennard-Jonescontribution. The origin of the z -axis is located atthe center of mass of the system.a marked dependence on the choice of the upperlimit of integration.At the highest investigated temperature of550 K, some of the systems, especially after longsimulation times, separate occasionally into two ormore liquid slabs, preventing to perform a mean-ingful fit of the density profiles with Eq. (2) and, atthe same time, preventing to compute the surfacetension with the help of the expression for planarinterfaces γ = L ( p N − p T ) , (3)where L is the length of the box edge parallel tothe surface normal, and p N and p T are the normaland lateral pressure components, respectively. Inthese cases, we did not determine the data for sur-face tension and densities of the liquid and vaporphases. In Fig.1 we report two snapshots of theTIP4P water/vapor interface at T =550 K, show-ing a case where the system is partitioned intotwo distinct liquid and vapor phases (top) and onecase where the system separated into two liquiddroplets surrounded by vapor (bottom). Results and Discussion
In Tab. 1 we report the average densities of the liq-uid ( ρ L ) and of the vapor ( ρ V ) phases for differ-ent temperatures, both for the case when disper-sion forces are treated with sPME, and when theyare truncated at r c = . T / K ρ / kg/m SPCwithout LJPMEwith LJPME
Figure 3: Density of the coexisting liquid and va-por phases for the SPC water model, as a functionof the temperature. Full symbols: results fromthe simulations with truncated dispersion forces;open symbols: results from the simulations withfull long-range dispersion forces. The error barsare much smaller than the symbols.from the simulation with truncated forces are re-ported as differences from the sPME results ( ∆ ρ L and ∆ ρ V ). The same convention is used for theother reported quantities (interfacial width d , sur-face tension γ ). With the exception of the TIP4P-2005 model, the inclusion of long-range disper-sion forces with sPME shifted the equilibrium den-sity of the liquid phase systematically to highervalues, and the density of the vapor phase to lowervalues.This effect is the more pronounced, the higherthe temperature, and can be easily appreciated byplotting the densities along the coexistence curve,as a function of the temperature (see Fig. 3 for theSPC model and the supplementary material for theother systems).The fact that the density of the liquid phase in-creases is a consequence of the stronger cohesiveforces arising by inclusion of long-range contribu-tions. As a consequence, the difference ρ L − ρ V increases, and this implies a shift of the criticaltemperature to higher values and, therefore, of ahigher surface tension, as we will see later. In ret-rospect, such a change in density could have beenexpected, as the analytical tail corrections for thesurface tension are known to be positive. Another effect of the increased cohesive forceis the change in interfacial width, which also de-creases systematically when the long-range dis-3 γ / m N / m SPCSPCETIP3Pexp. interp. γ / m N / m T /K TIP4PTIP4P2005TIP5PTIP5PEexp. interp. Figure 4: Values of the surface tension as a func-tion of the temperature for the seven simulated wa-ter models, when long-range contribution of thedispersion forces are taken into account.persion forces are taken into account. A generaltrend appears also in this case, with the interfacialwidth at higher temperature being characterized bylarger differences, although in absolute terms thechanges are rather small, in most cases being be-low one Angstrom, and at most half of a nanometer(see Tab. 1).The case of TIP4P-2005, as already mentioned,is qualitatively different from the other models, inthat not only the density of the liquid increasesupon introduction of the long-range dispersionforces, but also that of the vapor. This fact doesnot contradict the increase in surface tension, asthe difference between liquid and vapor densitiesis still larger due to the truncated dispersion forcecases, but we have found no explanation for thisopposite trend.The surface tension for the seven different watermodels, simulated including the long-range dis-persion forces, is reported in Fig. 4, together withthe continuous line that interpolates experimentalresults using the formula, γ = B τ µ ( + b τ ) , (4)where τ = − T / T c , T c = .
096 K, B = . b = − .
625 and µ = . − − − − − − − δ γ / m N / m TIP4PTIP4P2005TIP5P T /K TIP5PE Figure 5: Surface tension difference δ γ betweenthe results of the simulations with truncated dis-persion forces, corrected with Eq. (1) and thosefrom the simulations with full dispersion forcecomputed using sPME. Note that it holds also δ γ = γ tail − ∆ γ .one, passing through an inflection point locatedexperimentally at T / T c = − b ( − µ ) / ( + µ ) = .
929 and, for several rigid water molecules, at T / T c = . ± . In Tab. 1 we report thecomplete set of measured surface tensions, includ-ing those simulated with the truncated dispersionforce. The long sampling times of 50 ns allowedto reach a remarkably high accuracy, with standarddeviations in all cases below 0.2 mN/m.In Tab. 1 we report also the value of the ana-lytical correction estimated from the simulationswith truncated dispersion forces. In order to com-pare the two set of results, it is more convenient toplot the difference δ γ between the corrected sur-face tension in the truncated case, and the surfacetension obtained with the full-long range contri-butions. We have reported these values for the4 able 1: Various properties for the different models simulated with inclusion of long-range disper-sion forces and the difference with the same properties calculated with cutoff (analytical correctionsnot included). Temperatures are expressed in K, densities in kg/m , distances in nm, and surfacetensions in mN/m. model T ρ L ∆ ρ L ρ V ∆ ρ V d ∆ d γ ∆ γ γ tail300 975.2(2) 3.5 0.024(1) -0.007 1.644(3) -0.016 54.8(1) 3.1 4.3350 932.6(2) 4.1 0.297(4) -0.016 2.062(3) -0.023 46.5(1) 3.1 3.9SPC 400 877.9(2) 5.3 1.721(10) -0.042 2.592(4) -0.025 36.8(2) 2.5 3.3450 808.4(2) 5.9 6.328(19) -0.624 3.338(5) -0.047 26.8(2) 2.3 2.6500 717.4(2) 8.6 20.10(4) -1.520 4.666(5) -0.090 16.0(2) 2.1 1.8550 576.9(2) 23.3 70.22(27) -11.9 7.789(14) -0.578 5.6(1) 1.3 0.6300 996.8(2) 3.2 0.003(0) -0.004 1.473(4) -0.003 61.8(1) 2.8 4.5350 964.0(2) 3.5 0.103(2) -0.006 1.837(3) -0.007 54.3(2) 3.3 4.2SPC/E 400 919.1(2) 4.1 0.686(7) 0.043 2.269(3) -0.015 45.4(2) 2.7 3.7450 863.0(2) 5.1 2.700(13) -0.130 2.814(4) -0.024 35.5(2) 2.4 3.1500 792.9(2) 7.4 8.784(23) -0.480 3.620(4) -0.064 25.1(2) 1.9 2.4550 698.8(2) 10.2 25.61(4) -1.958 5.132(5) -0.154 14.5(2) 1.8 1.6300 983.1(2) 3.8 0.033(1) 0.003 1.733(3) -0.012 51.4(1) 3.4 4.6350 933.1(2) 4.2 0.359(4) -0.040 2.183(4) -0.012 43.0(1) 2.7 4.0TIP3P 400 871.5(2) 5.4 1.875(9) -0.090 2.751(4) -0.010 33.8(1) 2.5 3.4450 795.0(2) 7.5 7.116(19) -0.197 3.560(5) -0.060 23.7(1) 2.1 2.6500 694.2(2) 11.3 21.99(4) -2.064 5.158(6) -0.038 13.7(1) 2.0 1.6300 991.7(2) 3.7 0.031(1) -0.003 1.609(3) -0.010 56.5(1) 3.1 4.6350 953.3(2) 3.9 0.348(4) -0.046 2.024(3) -0.013 48.0(1) 3.3 4.2TIP4P 400 899.2(2) 4.6 2.049(11) -0.297 2.553(3) -0.021 38.3(2) 3.3 3.6450 828.8(2) 6.9 7.577(22) -0.483 3.326(4) -0.056 27.2(2) 2.3 2.8500 734.7(3) 10.1 24.289(36) -1.033 4.724(5) -0.136 16.0(2) 2.1 1.9550 582.4(3) - 82.47(28) - 8.316(17) - 4.9(1) - -300 995.3(3) 2.2 0.005(1) -0.080 1.407(4) -0.011 68.4(2) 3.8 4.6350 971.2(2) 3.8 0.072(2) 0.072 1.734(2) -0.013 60.9(2) 3.8 4.3TIP4P/ 400 932.7(2) 4.7 0.570(6) 0.570 2.130(3) -0.012 51.8(2) 3.4 3.82005 450 882.9(2) 6.0 2.320(12) 0.686 2.624(4) -0.015 41.6(2) 3.1 3.3500 820.0(2) 6.8 7.490(21) 0.446 3.325(4) -0.030 30.8(2) 2.7 2.7550 739.7(2) 9.1 20.88(3) -0.202 4.484(4) -0.121 19.6(2) 2.5 1.9300 981.8(2) 3.9 0.078(2) -0.006 1.750(3) -0.016 53.0(2) 2.9 4.6350 936.7(2) 4.5 0.730(6) 0.004 2.289(4) -0.022 41.5(2) 2.9 4.1TIP5P 400 863.7(2) 6.3 3.819(15) -0.094 3.030(5) -0.025 30.0(1) 2.4 3.3450 763.0(2) 8.9 13.978(28) -0.897 4.241(5) -0.082 18.0(1) 2.0 2.2500 614.0(3) 19.8 52.46(18) -7.128 7.165(10) -0.456 6.7(1) 1.2 0.9300 1000.1(2) 4.0 0.051(2) -0.012 1.712(3) -0.006 56.5(2) 3.9 4.5350 956.8(2) 4.3 0.597(6) -0.028 2.215(4) -0.023 45.3(2) 3.3 4.0TIP5PE 400 887.6(2) 6.0 3.075(13) -0.280 2.889(5) -0.022 33.4(1) 2.7 3.3450 794.8(2) 8.5 11.510(26) -0.359 3.957(5) -0.089 21.4(1) 2.1 2.4500 664.4(3) 15.6 38.033(61) -5.736 6.269(7) -0.212 9.8(1) 1.8 1.2 δ γ indicates that the analytical correctionis overestimating the real surface tension. In allcases, the analytical tail correction tends to over-estimate the surface tension at temperatures lowerthan (cid:39)
500 K, where a crossover seems to occur.At first sight this seems to contradict the data ondensity and interfacial width, which differ moreat high temperatures, but one should not forgetthat the tail correction has an explicit dependenceon the squared difference ( ρ L − ρ V ) , which di-minishes greatly the magnitude of the correctionat high temperature, when the density of the twophases tend to converge to the same critical value. Conclusions
We have calculated some interfacial properties ofseven popular water models with and without in-clusion of the long-range part of the dispersionforces. As the cohesion of water molecules in theliquid state is mainly determined by electrostaticforces, the binodal line is affected in a less pro-nounced way than in uncharged liquids, but shows,nevertheless, the clear tendency of the liquid phaseto be more dense, especially at high temperature,with respect to the simulations with truncated dis-persion forces. The analytical tail corrections tothe surface tension in simulations with truncatedforces (with a cutoff of 1.3 nm) yield results whichcan differ up to 1.5 mN/m from from those ob-tained using the full dispersion interaction, andthis can represent about 20-30% of the analyticalcorrection term itself at 300 K.The use of mesh Ewald methods like sPME forthe calculation of long-range dispersion forces ap-pears therefore to be advantageous even in caseof liquids like water, which are dominated by theelectrostatic interaction: it allows to remove thedependence of the liquid and vapor densities onthe choice of the cutoff, which represents certainlyan important step towards more transferable mod-els, and it also allows to incorporate autmaticallythe long-range contributions to the surface tensionalso for those geometries for which analytical cor-rections are not available.A final note on the computational cost of the in-troduction of long-range dispersion forces is due: compared to the simulations with cutoff, the per-formance of the simulations with long-range dis-persion forces decreased by 24% for the simplest,3 point charges models, and by 13% for the 5 pointcharges models, where the overhead associated tothe calculation of all other interactions is larger.The choice of the sPME parameters was, how-ever, not aimed at maximizing the performances,but rather, to achieve accurate results.
Acknowledgement
We acknowledge finan-cial support by the Austrian Science Foundation(FWF) within the SFB Vienna Computational Ma-terials Laboratory (Vicom, Grant F41) and byETN-COLLDENSE (H2020-MCSA-ITN-2014,Grant No. 642774). Calculations were carried outon the Vienna Scientific Cluster.
References (1) Hansen, J. P.; McDonald, I. R.
Theory of Sim-ple Liquids , 2nd ed.; Academic Press: Lon-don, UK, 1990.(2) de Leeuw, S. W.; Perram, J. W.; Smith, E. R.Simulation of Electrostatic Systems in Peri-odic Boundary Conditions. I. Lattice Sumsand Dielectric Constants.
Proceedings of theRoyal Society A: Mathematical, Physical andEngineering Sciences , , 27–56.(3) Ewald, P. P. Die Berechnung optischer undelektrostatischer Gitterpotentiale. Ann. Phys. , , 253–287.(4) Blokhuis, E.; Bedeaux, D.; Holcomb, C.;Zollweg, J. Tail corrections to the surfacetension of a Lennard-Jones liquid-vapour in-terface. Molecular Physics , , 665–669.(5) in ’t Veld, P. J.; Ismail, A. E.; Grest, G. S.Application of Ewald summations to long-range dispersion forces. Journal of ChemicalPhysics , , 144711.(6) Wennberg, C. L.; Murtola, T.; Hess, B.; Lin-dahl, E. Lennard-Jones Lattice Summationin Bilayer Simulations Has Critical Effectson Surface Tension and Lipid Properties. J.Chem. Theory Comput. , , 3527–3537.67) Hermans, J.; Berendsen, H. J. C.; van Gun-steren, W. F.; Postma, J. P. M. A consistentEmpirical Potential for Water-Protein inter-actions. Biopolymers , , 1513–1518.(8) Berendsen, H. J. C.; Grigera, J. R.;Straatsma, T. P. The missing term in effec-tive pair potentials. J. Phys. Chem. , ,6269–6271.(9) Jorgensen, W. L.; Chandrasekhar, J.;Madura, J. D.; Impey, R. W.; Klein, M. L.Comparison of simple potential functionsfor simulating liquid water. The Journal ofChemical Physics , , 926.(10) Abascal, J. L.; Vega, C. A general pur-pose model for the condensed phases of wa-ter: TIP4P/2005. The Journal of chemicalphysics , , 234505.(11) Mahoney, M. W.; Jorgensen, W. L. A five-site model for liquid water and the reproduc-tion of the density anomaly by rigid, nonpo-larizable potential functions. J. Chem. Phys. , , 8910–8922.(12) Rick, S. W. A reoptimization of the five-sitewater potential (TIP5P) for use with Ewaldsums. J. Chem. Phys. , , 6085.(13) Abraham, M. J.; Murtola, T.; Schulz, R.;Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E.GROMACS: High performance molecularsimulations through multi-level parallelismfrom laptops to supercomputers. SoftwareX , , 19–25.(14) Nosé, S. A molecular dynamics method forsimulations in the canonical ensemble. Mol.Phys. , , 255–268.(15) Hoover, W. G. Canonical dynamics: Equilib-rium phase-space distributions. Phys. Rev. A , , 1695–1697.(16) Miyamoto, S.; Kollman, P. A. SETTLE: AnAnalytical Version of the SHAKE and RAT-TLE Algorithms for Rigid Water Models. J.Comp. Chem. , , 952–962. (17) Essmann, U.; Perera, L.; Berkowitz, M. L.;Darden, T.; Lee, H.; Pedersen, L. G. Asmooth particle mesh Ewald method. J ChemPhys , , 8577–8593.(18) Wennberg, C. L.; Murtola, T.; Pall, S.; Abra-ham, M. J.; Hess, B.; Lindahl, E. Direct-Space Corrections Enable Fast and Ac-curate Lorentz-Berthelot Combination RuleLennard-Jones Lattice Summation. Journalof Chemical Theory and Computation , , 5737–5746.(19) Vega, C.; de Miguel, E. Surface tension of themost popular models of water by using thetest-area simulation method. J. Chem. Phys. , , 154707.(20) R. Piessens, E. de Doncker-Kapenga, C. W.Uberhuber, D. K. K. Quadpack: a SubroutinePackage for Automatic Integration ; SpringerScience & Business Media: Berlin, 1983.(21) Vargaftik, N. B.; Volkov, B. N.; Voljak, L. D.International Tables of the Surface Tensionof Water.
Journal of Physical and ChemicalReference Data , , 817–820.(22) Sega, M.; Horvai, G.; Jedlovszky, P. Micro-scopic origin of the surface tension anomalyof water. Langmuir , , 2969–2972.7 raphical TOC Entryraphical TOC Entry