Long-range Lennard-Jones and electrostatic interactions in interfaces: application of the isotropic periodic sum method

Abstract

Molecular dynamics (MD) simulations of heptane/vapor, hexadecane/vapor, water/vapor, hexadecane/water, and dipalmitoylphosphatidylcholine (DPPC) bilayers and monolayers are analyzed to determine the accuracy of treating long-range interactions in interfaces with the isotropic periodic sum (IPS) method. The method and cutoff (rc) dependences of surface tensions, density profiles, water dipole orientation, and electrostatic potential profiles are used as metrics. The water/vapor, heptane/vapor, and hexadecane/vapor interfaces are accurately and efficiently calculated with 2D IPS (rc=10 A). It is demonstrated that 3D IPS is not practical for any of the interfacial systems studied. However, the hybrid method PME/IPS [Particle Mesh Ewald for electrostatics and 3D IPS for Lennard-Jones (LJ) interactions] provides an efficient way to include both types of long-range forces in simulations of large liquid/vacuum and all liquid/liquid interfaces, including lipid monolayers and bilayers. A previously published pressure-based long-range LJ correction yields results similar to those of PME/IPS for liquid/liquid interfaces. The contributions to surface tension of LJ terms arising from interactions beyond 10 A range from 13 dyn/cm for the hexadecane/vapor interface to approximately 3 dyn/cm for hexadecane/water and DPPC bilayers and monolayers. Surface tensions of alkane/vapor, hexadecane/water, and DPPC monolayers based on the CHARMM lipid force fields agree very well with experiment, whereas surface tensions of the TIP3P and TIP4P-Ew water models underestimate experiment by 16 and 11 dyn/cm, respectively. Dipole potential drops (DeltaPsi) are less sensitive to long-range LJ interactions than surface tensions. However, DeltaPsi for the DPPC bilayer (845+/-3 mV proceeding from water to lipid) and water (547+/-2 mV for TIP4P-Ew and 521+/-3 mV for TIP3P) overestimate experiment by factors of 3 and 5, respectively, and represent expected deficiencies in nonpolarizable force fields.

Figure 1.

Schematic of the two types interfaces examined in this study.

Figure 2.

Density of the heptane/vapor interface at 298.15 K for 2D IPS, 3D IPS, and PME for selected values of the cutoff r_c in angstroms.

Figure 3.

Orientation of the water dipole with respect to the positive z axis scaled by the electron density for PME (top), 2D IPS (middle), and 3D IPS and PME/IPS (bottom) for TIP4P-Ew at 323.15 K. Results for TIP3P are included for 2D IPS (middle panel).

Figure 4.

Electrostatic potential profiles of the water/vapor interface for TIP4P-Ew at 323.15 K (top) and the hexadecane/water interface at 323.15 K (bottom; the dotted line is for TIP4P-Ew, and the remainder is TIP3P). In both panels, simulations for hexadecane/water were run with an interface of water/hexadecane/water, and the values shown in this figure were shifted to place water at the center.

Figure 5.

Orientation of the water dipole with respect to the positive z axis, scaled by the electron density, for hexadecane/water at 323.15 K. The two PME/IPS results are with r_c = 10 Å, and the PME result is with r_c = 25 Å.

Figure 6.

Electron density of the DPPC bilayer at 323.15 K for PME/ IPS with r_c = 10 Å (red), for PME/LRC (green), and for PME(r_c=10 Å) (black) compared to the experimental H2 model density (blue). The individual component densities are also shown: CH₃ = methyl, CH₂ = methylene, CG = carbonyl-glycerol, P = phosphate, and W + Chol = water and choline.

Figure 7.

Simulated DPPC lipid monolayer surface tension at three surface areas for PME/IPS and PME with r_c = 10 Å. The experimental data are also shown for 323.15 and 321.15 K from refs 31 and 32, respectively.

Figure 8.

(Top) Orientation of the water dipole with respect to the positive z axis, which is scaled by the electron density for DPPC at 323.15 K. (Bottom) Electrostatic potential profile of the fully hydrated DPPC bilayer. A value of r_c = 10 Å is used for PME/IPS and PME. The bilayer center is located at z = 0.

Figure 9.

Electrostatic potential profile of DPPC monolayers for PME/ IPS and PME with r_c = 10 Å. For clarity, a comparison is made between PME/IPS (solid) and PME (dashed) only at A = 64 Å ²/lipid. The center of the water layer is located at z = 0.