Pressure Profile Calculation with Mesh Ewald Methods
PPressure Profile Calculation with Mesh Ewald Methods
Marcello Sega, ∗ , † Bal´azs F´abi´an, ‡ , ¶ and P´al Jedlovszky § , (cid:107) , ⊥ † Computational Physics Group, University of Vienna, Sensengasse 8/9, 1170 Vienna, Austria ‡ Institut UTINAM (CNRS UMR 6213), Universit´e de Franche-Comt´e, 16 route de Gray, F-25030Besan¸con, France ¶ Department of Inorganic and Analytical Chemistry, Budapest University of Technology and Economics,Szt. Gell´ert t´er 4, H-1111 Budapest, Hungary § EKF Department of Chemistry, Le´anyka u. 6, H-3300 Eger, Hungary (cid:107)
Laboratory of Interfaces and Nanosize Systems, Institute of Chemistry, E¨otv¨os Lor´and University,P´azm´any P. Stny 1/A, H-1117 Budapest, Hungary ⊥ MTA-BME Research Group of Technical Analytical Chemistry, Szt. Gell´ert t´er 4, H-1111 Budapest,Hungary
E-mail: [email protected]
Abstract
The importance of calculating pressure profiles across liquid interfaces is increasingly gaining recognition,and efficient methods for the calculation of long-range contributions are fundamental in addressing systemswith a large number of charges. Here, we show how to compute the local pressure contribution for mesh-basedEwald methods, retaining the typical N log N scaling as a function of the lattice nodes N . This is a considerableimprovement on existing methods, which include approximating the electrostatic contribution using a largecut-off and the, much slower, Ewald calculation. As an application, we calculate the contribution to thepressure profile across the water/vapour interface, coming from different molecular layers, both including andremoving the effect of thermal capillary waves. We compare the total pressure profile with the one obtainedusing the cutoff approximation for the calculation of the stresses, showing that the stress distribution obtainedby the Harasima and Irving-Kirkwood are quite similar and shifted with respect to each other at most 0.05 nm. Calculating the profiles of various physical quantitiessuch as mass, charge, or electrostatic potential, is a stan-dard approach used in computer simulations in orderto characterize inhomogeneous systems including fluid(soft) interfaces, like the water/vapour one depictedin Fig.1, membranes, micelles or other self-associates.While some of these profiles have been routinely cal-culated in the past decades, the importance of deter-mining pressure profiles is only recently being recog-nized, for example, for the calculation of the mechanicalproperties of macromolecules.
Conjectures related topressure profiles also play a key role in the possible ex-planations of several phenomena. Thus, for instance,it was proposed by Cantor that the molecular mecha-nism of anesthesia is related to the alteration of the lat-eral pressure profile inside the cell membrane; Imre etal. claimed that the spinodal pressure of a liquid phasecan be extracted from the lateral pressure profile ob-tained at the liquid-vapour interface of the same systemat the same temperature. Experimental verification orfalsification of these conjectures is an extremely chal-lenging task as it requires the measurement of the meanlateral pressure across interfaces, or macromolecules at atomic resolution, a task that is rarely accomplished(see, for example Ref. ).The calculation of pressure profiles in computer sim-ulations, on the other hand, is complicated by the factthat it requires the local determination of a quantitythat is inherently non-local. More precisely, a localpressure tensor cannot, in general, be uniquely defined.Instead, it can only be given up to a divergence-freesecond rank tensor in a path dependent form. Severalintegration paths, such as the Irving-Kirkwood andHarasima paths, were shown to provide comparableprofiles. Besides its conceptual simplicity, the Irving-Kirkwood path has the advantage of leading to a for-mula for the local pressure that gives access to both thenormal and lateral components of the virial, distributedhomogeneously along the paths connecting pairs of par-ticles. On the other hand, while the Harasima pathallows one to calculate only the lateral components ofthe pressure tensor, these are concentrated at the par-ticles’ positions. In other words, a certain amount ofthe total lateral pressure is associated with each parti-cle, making the calculation of the lateral pressure profilestraightforward and computationally efficient.
Another serious technical problem in calculating pres-sure profiles in charged systems is how to take into ac-1 a r X i v : . [ c ond - m a t . s o f t ] J un igure 1: Simulation snapshot of a water/vapour interface.The first five molecular layers in the upper half of the simulationbox (lower half not shown) are highlighted with different colors.One molecule identified to be in the vapour phase by a cluster-search algorithm can be also seen above the first layer of thesurface count the long range correction term of the electrostaticinteraction. It is now well-established that neglectingthis correction to the energy and forces can lead to im-portant systematic errors. The method of reaction fieldcorrection, which uses the dielectric constant of thecontinuum beyond a suitably chosen cutoff sphere, isunsuitable for strongly inhomogeneous systems. TheEwald summation can, in principle, be used in suchsystems; the contribution of its correction term to thepressure profile was derived by Sonne and coworkers forthe Harasima path. They also showed that potentialsthat are not strictly pairwise additive, such as the cor-rection term of the Ewald summation, cannot be usedin combination with the Irving-Kirkwood profile. How-ever, because of its poor scaling properties, the use ofthe Ewald summation is prohibitive for systems with alarge number of charges.In theory, this problem can be overcome using amethod based on the Fast Fourier Transform (FFT),such as mesh Ewald methods.
However, the localpressure from the reciprocal space contribution of meshEwald methods has not, to the best of our knowledge,been derived yet.For approaches based on the Irving-Kirkwood path,as it is not possible to take into account non-pair addi-tive potentials, several authors have opted for runningthe simulations using their method of choice for the eval-uation of long range contributions, but considering cut-off electrostatics for the purpose of evaluating the localpressure (see, for example, Refs. ). This is usu-ally done by re-running the simulation and using cutoffvalues that are larger than usual, thus reducing the ar-tifacts. However, in order to reach a precision of theorder of 1 bar, cutoff values larger than 2.5 nm mustbe used (see Fig. 7 in Ref. ). Sonne and collabora-tors showed that the pressure profile calculated usingthe Irving–Kirkwood approach and a cutoff radius of 2nm is qualitatively similar to the one calculated using the Harasima path with long-range electrostatic contri-butions, although with an indetermination of roughly100 to 200 bar, which makes it only a qualitative test.In this paper we show how the correction term of thesmooth Particle Mesh Ewald (sPME) method can betaken into account in calculating lateral pressure pro-files using the Harasima path in an efficient way thatdoes not alter the scaling of the long-range correctioncalculation. We describe the method itself in Sec. 2and compare the pressure contributions associated toparticles with those obtained with a reference, simpleEwald calculation. In Sec. 3 we apply the method tothe case of the water/vapour interface, and in Sec. 4we describe the scaling properties of the algorithm. Fi-nally, in Sec. 5, we summarize the results presented inthis article. We start by decomposing the instantaneous pressuretensor P , as customary, in the ideal gas and virial ( Ξ )contributions, P V = (cid:88) i m i v i ⊗ v i − Ξ (1)The difference between lateral and normal componentsof the pressure tensor plays an important role for inter-faces, as it describes the surface tension, and its localvariant can be written as γ ( z ) = P N − P T ( z ) , (2)where we have chosen P N = P zz and P T = ( P xx + P yy ) /
2. This formula applies only to planar interfacesand we will assume, as it is customary, that it appliesalso in the presence of capillary waves, as long as thesurface remains macroscopically flat. The presence ofcurved interfaces introduces additional difficulties in thedefinition and calculation of the surface tension.
Here we note that while P N is constant through thesystem to ensure mechanical stability, Ξ zz is not.In the Ewald sum, the reciprocal space contributionsto energy and virial are expressed as sums over the re-ciprocal lattice vectors m U rec = (cid:88) m (cid:54) =0 πV f ( m ) S ( − m ) S ( m )Ξ rec µν = 1 V (cid:88) m (cid:54) =0 πV f ( m ) g µν ( m ) S ( − m ) S ( m ) (3)where the functions f and g µν are defined as f ( m ) = e − π m /β m (4) g µν ( m ) = (cid:18) δ µν − π m /β m m µ m ν (cid:19) , (5)2nd the structure factor is S ( m ) = (cid:88) i q i e πi m · r i . (6)In mesh-based algorithms, the reciprocal space term inthe Ewald sum is calculated with the aid of the FFT by replacing the point charges, located in the contin-uum at positions r i , with a discretized distribution ona lattice with nodes at positions r p ˜ ρ ( r p ) = 1 h (cid:88) i W ( r p − r i ) , (7)where h is the lattice spacing and W ( r p − r i ) is a suit-able charge assignment function. The structure fac-tor S ( m ) is then replaced in Eqs.(3) by its estimate˜ S ( m ) = B ( m )FFT[˜ ρ ], which is now evaluated on thepoints of the reciprocal lattice, and where the function B ( m ) depends on the specific interpolation scheme. The choice of the assignment function for the discretecharge distribution is not unique, a fact that has ledto the appearance of families of methods such as theParticle-Particle-Particle-Mesh (P3M) of Hockney andEastwood, and the Particle Mesh Ewald (PME) in itsoriginal variant that uses Lagrange interpolation orin its “smooth” version (sPME) based on cardinal B-splines. Once a physical quantity pertaining to singleparticles has been calculated on the lattice, it has to beinterpolated back from the lattice to the real positionof the particles. This is usually done by the same as-signment function used to generate the distribution onthe lattice. If a function A ( r p ) is known at the latticenodes, it can be distributed back to the atomic positionsusing A ( r i ) = (cid:88) r p A ( r p ) W ( r i − r p ) . (8)The reciprocal space term of the virial in the Ewaldsum can be easily written in a form that is suitablefor the Harasima path formulation, which associates toparticle i the contribution Ξ rec , i µν = q i V (cid:88) m (cid:54) =0 Re (cid:26) e − i m · r i πV S ( m ) (cid:27) f ( m ) g µν ( m ) . (9)As we are interested in the diagonal elements of thetensor only, µ = ν , the real part operator in Eq.(9) issuperfluous because the f and g µν are even functions of m . Another approach has been previously used in liter-ature to take into account the reciprocal space contribu-tion to the local pressure for the plain Ewald method, but consisted of distributing the reciprocal space virialcontribution equally to each particle, which is not a jus-tified assumption, as the reciprocal space force contri-bution is different for each particle.Note that, because the contribution in reciprocalspace is not pairwise additive, it is not possible to use anIrving-Kirkwood formulation either for the Ewald sum,or for its mesh-based approximations. We now proceed to extending the reciprocal spacecontribution to mesh-based algorithms, by noting thatEq.(9) is, for µ = ν , the discrete inverse Fourier trans-form of the complex function S ( m ) f ( m ) g µν ( m ), evalu-ated at the particles’ positions, and scaled by the factor q i . On the lattice, this expression is replaced by theinverse FFT of the function ˜ S ( m ) f ( m ) g µν ( m ). Thisquantity is eventually interpolated back to the real po-sitions using Eq.(8), leading to the approximate expres-sion for the reciprocal space contribution˜Ξ rec , i µµ = q i V (cid:88) r p ˜Ξ rec µµ (r p ) W ( r i − r p ) , (10)with˜Ξ rec µµ (r p ) = FFT − (cid:2) B ( m )FFT[˜ ρ ] f ( m ) g µµ ( m ) (cid:3) . (11)Notice that in order to be consistent with the sumEq.(3), the argument of the inverse FFT in Eq.(11) for m = 0 must be set explicitly to zero. The correction andexclusion terms in real space are expressed as sum ofpair contributions and, thus, must be taken into accountsimilarly to other pairwise, short range interactions.We would like to stress once more that, thanks to theuse of the Harasima formulation, it is possible to as-sociate a virial contribution to each particle, and thelocal stresses are thus distributed in the continuum.In the present case, a grid is used only for the mesh-Ewald charge spreading procedure and, once the back-interpolation is performed, the reciprocal space contri-bution is associated to each of the particles.We implemented this approach in our de-velopment version (available free of charge at github.com/Marcello-Sega/gromacs/tree/virial )of the popular GROMACS molecular simulation pack-age, version 5.0, which extends the trajectoryfiles (usually containing positions, velocities and forces)to include also the local virial components Ξ xx , Ξ yy andΞ zz .In order to check for the correctness and proper im-plementation of the algorithm, we tested it on a sys-tem of point charges (in absence of any other type ofinteraction but the electrostatic one). As a first step,we checked that the sum of the pressure terms associ-ated to the atoms indeed yields the global virial up toroundoff errors. Subsequently, we computed the rootmean square error ∆ P per particle, following Holm andDeserno, as∆ P = (cid:115) N (cid:88) i (cid:88) µ (cid:16) P iµµ − ˆ P iµµ (cid:17) , (12)where N is the number of charges and ˆ P iµµ is a refer-ence value, namely the pressure contribution of particle i calculated on the same configurations using the plainEwald method, with a large number of reciprocal spacecontributions (12 along x and y direction, 36 along z di-rection), and a relative accuracy at the cutoff (1.5 nm)3 − − − − ∆ P / q h P i mesh spacing / nm Figure 2: Root mean square error per particle for the pressure,as a function of the mesh spacing (charge interpolation orderset to 4). − − − − ∆ P / q h P i charge interpolation order Figure 3: Root mean square error per particle for the pressure,as a function of the charge interpolation order (mesh spacingset to 0.15 nm). of 10 − . In Fig. 2 and Fig. 3 we report, respectively, therelative root mean square error ∆ P/ (cid:112) (cid:104) P (cid:105) per parti-cle, as a function of the mesh spacing and of the chargeinterpolation order, where (cid:104) P (cid:105) = N (cid:80) i (cid:80) µ ˆ P i µµ .The average error that is made in estimating the pres-sure contribution is less than 0 .
5% for all choices of pa-rameters and decreases monotonically with decreasingmesh spacing and increasing interpolation order, downto 0 . We proceeded to run a simulation of the water/vapourinterface, modeling water molecules with the SPC/E potential. While other models compare better with sev-eral experimentally known quantities, the SPC/Emodel represents an optimal choice for testing newmethods, thanks to its modest computational require-ments.The initial configuration was generated simply by in-creasing the z box edge of a 5 × × equili-brated water simulation box to 15 nm. To integratethe equations of motion we used an integration step of1 fs, constraining all bonds by means of the SHAKEalgorithm (the local pressure contribution of the con-straints were being also taken into account), and keep-ing the temperature at the constant value of 300 K bymeans of the Nos´e–Hoover thermostat. The short-range part of the potentials were cut off either at 1.5 orat 2.0 nm and, at that distance, the relative accuracyof the direct space contribution to the electrostatic po-tential was set to 10 − . The long-range correction tothe electrostatic potential was calculated using sPMEfor all trajectories, with a target lattice spacing of 0.15nm. No long-range corrections were used for the vander Waals potential. Configurations, including atomicpositions and virial contributions, were saved to diskevery 0.1 ps. The contributions to the lateral virialstemming from the Lennard-Jones interaction and fromconstraint forces were taken into account as usualwith the expression of the virial for the Harasima pathΞ iµµ = − / (cid:80) j (cid:54) = i f ijµ r ijµ , where f ijµ and r ijµ are the com-ponents of the pair forces and of the vectors connectingthe two interacting atoms, respectively.We calculated the surface tension profile with the Ha-rasima path, including the full sPME contributions andLennard-Jones interactions cut off at r c = 1 . using the Irving-Kirkwoodpath and cut-off values of 1.5 and 2.0 nm, respectively.The resulting profiles are reported in Fig.4. The twoprofiles calculated using the Irving-Kirkwood path arequalitatively similar, differing only in the peak height.On the other hand, the profile calculated with the Ha-rasima path is shifted by about 0.05 nm towards thecenter of the fluid slab. From the surface tension profiles γ ( z ) it is possible to calculate both the surface tension γ = (1 / (cid:82) γ ( z ) dz as well as the position of the surfaceof tension z s = (1 / γ ) (cid:82) γ ( z ) | z | dz . The values of thesurface tension γ as computed from the integral of theprofile (thus neglecting the long-range contributions inthe Irving-Kirkwood case) and from the global pressuretensor are reported in Tab. 1, together with the valuesof z s and of its distance δ = z e − z s (the Tolman length)from the location of the equimolar surface z e . Theequimolar surface is defined as the one that divides avolume V in two regions (liquid and vapour) of volumes4 [ P N − P T ( z ) ]/ k ba r ( z − z cm ) / nmIrving-Kirkwood, 2.0 nmIrving-Kirkwood, 1.5 nmHarasima, 1.5 nm + sPME Figure 4: Comparison between the surface tension profile γ ( z ) obtained with two different paths (Irving-Kirkwood and Ha-rasima) and two different cut-off values for the short-rangeforces and (in case of the Irving-Kirkwood path) for the cal-culation of the electrostatic contribution to the pressure. Theprofile is the average of the two halves of the simulation box,and the origin is in the middle of the water slab. V l and V v , respectively, such that ρV = ρ v V v + ρ l V l , with ρ v and ρ l the densities of the two phases as measuredfar from the interface. Since we have two interfaces,we apply this criterion to half of the simulation box (theliquid phase being centerd in the box), and obtain that ρ l z e + ρ v ( L/ − z e ) = ρL/
2, or z e = L ( ρ − ρ v ) / ρ l − ρ v ).Since in our case ρ v /ρ (cid:39) ρ v /ρ l (cid:39) − , we can safelyuse the approximation z e (cid:39) ρL/ (2 ρ l ).The equimolar surface was found to be located at z e =2 . ± .
01 and 2 . ± .
01 nm for the 1.5 and 2.0 nmcutoff cases, respectively.For the sake of comparison, in Tab. 1 we report somedata from the literature, obtained with different meshparameters and cutoff treatment, for the surface ten-sion of the SPC/E model at 300 K. The values reportedhere are those which do not include analytical tail cor-rections, as in the present work.It is interesting to note that analytical results esti-mate the value of δ to be in the range of the molecularsize r m (at least, for pair potentials). More precisely, δ (cid:39) r m / r m / This means, in the case ofwater, δ (cid:39) . extending its analysis capabilities to the Table 1: Surface tension (bar nm) as derived from the integralof the surface tension profile ( γ † ) and from the global pressuretensor components ( γ ), position of the surface of tension z s (nm), and Tolman length δ (nm). The quantities have beencalculated using the Irving-Kirkwood (I-K) and Harasima (H)paths, for different values of the cutoff radius r c . r c γ † γ z s δ I-K 1.5 565 588 ± ± . ± ± .
05H 1.5 588 588 ± ± . − − − − − − [ P N − P T ( z ) ]/ k ba r ( z − z cm ) / nm Figure 5: Lateral pressure profile ( P N − P T ( z ) ) of the wa-ter/vapour interface. The contribution of each of the first 5layers is reported (shifted along the vertical axis by 0.5 kbareach), using the same color scheme as for the simulation snap-shot, Fig.1. The total profile is the curve at the top. extended GROMACS trajectory format that includesvirial contributions (this software is also freely available,at marcello-sega.github.io/gitim/ ). The interfa-cial analysis was performed on a molecular basis (i.e. ifat least one atom in the molecule is at the interface, thecomplete molecule is considered to be an interfacial one)using a probe sphere radius of 0.2 nm. Vapour phasemolecules were identified using a neighbor list cluster-search algorithm based on a simple oxygen-oxygen dis-tance cutoff of 0.35 nm (corresponding to the first min-imum of the oxygen-oxygen radial distribution functionof water). A simulation snapshot, with molecules fromdifferent molecular layers highlighted in different colors,as well as molecules in the vapour phase as identified bythe cluster-search, is reported in Fig.1.The presence of successive layers below the surfaceone is not just a feature encountered in solids, but isalso present in liquids, and can be clearly seen by cal-culating the density profile in a local reference framelocated at the interface, effectively removing the smear-ing introduced by the presence of capillary waves (whichis more important the larger the allowed wavelengths inthe simulation box are). If the interface position5 − − − − − − − − − [ ( P N − P T ( z )) ]/ k ba r ( z − ξ ) / nm Figure 6: Intrinsic lateral pressure profile ( P N − P T ( z ) ) of thewater/vapour interface. The contribution of each of the first 5layers is reported (shifted along the vertical axis by 3 kbar each),using the same color scheme as for the simulation snapshot,Fig.1. The total profile is the curve in the inset. Delta-likecontributions at the interface ( z = ξ ) are replaced by a grayband in the corresponding bin. is z = ξ ( x, y ), the intrinsic density profile can then bedefined in terms of the position ( x i , y i , z i ) of the i -thparticle as ρ I ( z ) = 1 A (cid:42)(cid:88) i δ ( z − z i + ξ ( x i , y i )) (cid:43) , (13)where A is the simulation box cross-sectional area alongthe macroscopic interface normal vector. Similarly, wedefine the intrinsic pressure profile as γ I ( z ) = P N − VA (cid:42)(cid:88) i P iT δ ( z − z i + ξ ( x i , y i )) (cid:43) . (14)Having access to the virial contribution of all parti-cles, it is straightforward to compute the pressure pro-file, both with respect to the distance from the center ofmass, (non-intrinsic profile) and with respect to the lo-cal position of the interface (intrinsic profile). In Fig.5,the total (top curve) value of the profile of γ ( z ) is re-ported, together with its decomposition into the suc-cessive layers contributions, up to the fifth layer. It isimportant to notice that while the total normal pressureprofile must be a constant to ensure mechanical stabil-ity, this is not true for the separate layer contributions.As we do not have direct access to the normal compo-nent using the Harasima path, the layer decompositionhas to be understood, strictly speaking, as that of thelateral pressure, which we offset by the value P N forconvenience. Nevertheless, this decomposition is stillinstructive, as one can find that the first layer profilecontributes the most to the total profile, even thoughthe latter is more narrowly distributed than the former. The net contribution of each of the layers far from theinterface (i.e., the integral of the surface tension distri-bution) is close to zero, consistently with the fact thatin this region the fluid is homogeneous and isotropic.The distribution itself, however, is far from being iden-tically zero, and reflects the fact that in the oscillatinglayer, atoms are undergoing different stresses dependingon their distance from the average position.The width of the contribution of the different layersis, on the other hand, quite constant and rather closeto the size of the first water coordination shell (of ra-dius 0.35 nm). Note that the distribution is normal-ized such that its integral, from the center of the liquidslab ( z = 0) up to the box edge length, yields the sur-face tension. The intrinsic profile γ I ( z ), Fig. 6, givesanother point of view on the pressure profile distribu-tion. The intrinsic profile is calculated as a functionof the distance from the local surface position, ξ ( x, y ),according to Eq.(14). In this way, the effect of thermalcapillary waves, which are corrugating the surface, isremoved. Note that the delta-like contribution at thesurface has been removed, and the corresponding binblanked. From a comparison between the two plots,one can appreciate the considerable narrowing of thepressure distribution, which decreases from 0.5 nm inthe non-intrinsic case, to 0.3 nm in the intrinsic case,clearly showing the smearing effect of capillary waves.A striking difference between the pressure profile (bothintrinsic and not) of a molecular liquid such as water,and that of an atomic liquid like argon, is that thelayers’ contributions are much less separated. Since thekissing condition between the probe sphere and a sur-face atom takes into account the atomic radius, some(point-like) hydrogen atoms in the second layer will belocated further than an oxygen atom of the first layer,contributing to the enhanced layer interpenetration. Inaddition, since we are considering molecular layers, thecontribution to the first layer comes also from atomsthat are not right at the surface and, therefore, the firstlayer contribution is not just a delta function, as it istypical for simple liquids. The second, qualitative dif-ference with respect to simple liquids, is that the majorpart of the lateral tension arises from atoms located in anarrow interval of 0.3 nm around the interfacial centres,a range which is comparable with one molecular diame-ter. The first molecular layer, in particular, contributesto more than 90% of the total surface tension. In simpleliquids, the distribution of the lateral pressure extendsto almost two atomic diameters, and the contributionof the first layer is slightly less than 80%.
The present algorithm for the calculation of the recipro-cal space virial contributions follows closely that of thecalculation of forces in mesh-Ewald methods. There-fore, the scaling is of the type N log( N ). Regarding theprefactor, it is easy to make an estimate of the rela-6 N / N N log( N ) e x e c u t i on t i m e / s number of chargesPME + particle pressurePMEEwald Figure 7: Time needed to integrate 100 steps for systems of dif-ferent sizes, from about 500 to more than water molecules.Circles and squares refer to the sPME execution time with-out and with local pressure calculation, respectively. Trianglesrepresent the execution time when a plain Ewald method isused. The solid thin lines are the result of a fit to the function N log( N ) . The typical scaling behaviors for the plain Ewaldcalculation N / and N are also shown. Inset: the Ewaldscaling in linear scale. tive computational load with respect to the case whereonly forces are being computed. The functions f ( m )and g µν ( m ) are usually being pre-calculated in orderto provide energy and forces, and this fact means thatfor the pressure profile computation, the additional cal-culations that need to be performed are one additionalinverse FFT per tensor element, and the correspondingback-interpolation. In terms of computational cost, thecalculation of each of the tensor elements has the sameload as the force calculation in the sPME method (whichrequires only one FFT, at a difference with the PMEmethod, which requires three). Since we calculate onlythe diagonal elements of the virial tensor, the theoreticalprefactor is precisely 4 times larger than for a standardsPME calculation. This is, however, the scaling onewould expect for the reciprocal part of the mesh-basedalgorithm alone. As a matter of fact, van der Waalsinteractions, as well as the short-range part of the elec-trostatic potential, are always present in a typical sim-ulation, with the result that in a real simulation setup,this factor of four represents only an upper limit of themeasured scaling. To test the actual performance of thecode, we ran a series of simulations of bulk water at adensity of about 1 g / cm with different box sizes, witha water content ranging from about 500 to more than10 water molecules, with the simulation protocol andparameters already used for producing the pressure pro-files, but a cut-off radius of 0.9 nm. In Fig.7, we reportthe time required to perform 10 integration timesteps,with and without calculating the local pressure contri-bution. For comparison, the timings of runs using theplain Ewald method are also reported. The result of afit to a N log( N ) scaling gives, in the large N limit, the prefactors are 5 . . µ s/timestep, with and withoutpressure calculation, respectively. The simulations in-cluding the calculation of the local pressure with sPMEare in this case 2.7 times slower than without includingthe pressure calculation. It should be noted, however,that in most practical cases it is not necessary to com-pute the local pressure at every timestep and, therefore,the method would introduce only a little overhead withrespect to a regular simulation. We presented a method to calculate the long-range con-tribution of the electrostatic interaction to the localpressure tensor, based on mesh-based algorithms. Weimplemented our approach in the GROMACS molecu-lar simulation package and have tested it on a model ofwater/vapour interface. We have calculated the lateralpressure profile with our approach, and compared it tothe one computed using the Irving-Kirkwood path, re-porting also the location of the surface of tension for thedifferent models. We have also shown how it is possibleto take advantage of the Harasima path decompositionto calculate the layer-by-layer contribution to the in-trinsic and non-intrinsic lateral pressure profile. Thepresent method is particularly efficient, as it displaysthe same scaling (in this case, N log( N )) of the algo-rithm used to compute long-range correction to forcesand energy.This project is supported by the Hungarian OTKAFoundation under project No. 109732, by the Ac-tion Austria Hungary Foundation under project No.93¨ou3104 and through ETN COLLDENSE, GrantH2020-MSCA-ITN-2014 No. 642774. The computa-tional results presented have been achieved in using theVienna Scientific Cluster (VSC). We thank M. McCaf-frey for proofreading the manuscript. References (1) Gullingsrud, J.; Schulten, K.
Biophys. J. , ,3496–3509.(2) Ollila, O. S.; Risselada, H. J.; Louhivuori, M.; Lin-dahl, E.; Vattulainen, I.; Marrink, S. J. Phys. Rev.Lett. , , 078101.(3) Hatch, H. W.; Debenedetti, P. G. J. Chem. Phys. , , 035103.(4) Torres-S´anchez, A.; Vanegas, J. M.; Arroyo, M. Phys. Rev. Lett. , , 258102.(5) Cantor, R. S. Biochemistry , , 2339–2344.(6) Imre, A.; Mayer, G.; H´azi, G.; Rozas, R.;Kraska, T. J. Chem. Phys. , , 114708.(7) H. Templer, R.; J. Castle, S.; Rachael Curran, A.;Rumbles, G.; R. Klug, D. Faraday Discuss. , , 41–53.(8) Irving, J.; Kirkwood, J. G. J. Chem. Phys. , , 817–829.(9) Harasima, A. Adv. Chem. Phys , , 203–237.710) Schofield, P.; Henderson, J. Proc. Royal Soc. Lon-don. A , , 231–246.(11) Sonne, J.; Hansen, F. Y.; Peters, G. H. J. Chem.Phys. , , 124903.(12) Sega, M.; F´abi´an, B.; Jedlovszky, P. J. Chem.Phys. , , 114709.(13) Onsager, L. J. Am. Chem. Soc. USA , ,1486–1493.(14) Ewald, P. Ann. d. Phys. , , 253.(15) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. , , 10089–10092.(16) Essmann, U.; Perera, L.; Berkowitz, M. L.; Dar-den, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. , , 8577–8593.(17) Hockney, R. W.; Eastwood, J. W. Computer sim-ulation using particles ; Taylor & Francis: NewYork, 1988.(18) Kasson, P. M.; Hess, B.; Lindahl, E.
Chem. Phys.Lipids , , 106–112.(19) Ollila, S.; Hyv¨onen, M. T.; Vattulainen, I. J. Phys.Chem. B , , 3139–3150.(20) Vanegas, J. M.; Torres-S´anchez, A.; Arroyo, M. J.Chem. Theory Comput. , , 691–702.(21) Rowlinson, J.; Widom, B. Molecular Theory ofCapillarityOxford Univ. Press ; Oxford UniversityPress: Oxford, 1982.(22) Sodt, A. J.; Pastor, R. W.
J. Chem. Phys. , , 234101.(23) Cooley, J. W.; Tukey, J. W. Math. Comput. , , 297–301.(24) Deserno, M.; Holm, C. J. Chem. Phys. , ,7678–7693.(25) Alejandre, J,; Tildesley, D J,; Chapela, G A, J.Chem. Phys. , , 4574–4583.(26) Berendsen, H. J.; van der Spoel, D.; vanDrunen, R. Comput. Phys. Comm. , , 43–56.(27) P´all, S.; Abraham, M.; Kutzner, C.; Hess, B.; Lin-dahl, E. In Solving Software Challenges for Exas-cale ; Markidis, S., Laure, E., Eds.; Lecture Notesin Computer Science; Springer International Pub-lishing: Cham (ZG), Switzerland, 2015; Vol. 8759.(28) Berendsen, H.; Grigera, J.; Straatsma, T.
J. Phys.Chem. , , 6269–6271.(29) Vega, C.; Abascal, J.; Nezbeda, I. J. Chem. Phys. , , 034503.(30) Alejandre, J.; Chapela, G. A. J. Chem. Phys. , , 014701.(31) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. J.Comput. Phys. , , 327–341.(32) Nos´e, S. Mol. Phys. , , 255–268.(33) Hoover, W. Phys. Rev. A , , 1695–1697.(34) Vega, C.; de Miguel, E. J. Chem. Phys. , ,154707.(35) Chen, F.; Smith, P. E. J. Chem. Phys. , ,221101.(36) P´artay, L. B.; Hantal, G.; Jedlovszky, P.;Vincze, ´A.; Horvai, G. J. Comput. Chem. , , 945–956. (37) Sega, M.; Kantorovich, S. S.; Jedlovszky, P.;Jorge, M. J. Chem. Phys. , , 044110.(38) Chac´on, E.; Tarazona, P. Phys. Rev. Lett. ,91