Treatment Of Electrostatic Interactions Research Articles

An Iterative Fragment Scheme for the ACKS2 Electronic Polarization Model: Application to Molecular Dimers and Chains.

The treatment of electrostatic interactions is a key ingredient in the force field-based simulation of condensed phase systems. Most approaches used fixed, site-specific point charges. Yet, it is now clear that many applications of force fields (FFs) demand more sophisticated treatments, prompting the implementation of charge equilibration methods in polarizable FFs to allow the redistribution of charge within the system. One approach allowing both, charge redistribution and site-specific polarization, while at the same time solving methodological shortcomings of earlier methods, is the first-principles-derived atom-condensed Kohn-Sham density functional theory method approximated to the second order (ACKS2). In this work, we present two fragment approaches to ACKS2, termed f-ACKS2 and a self-consistent version, scf-ACKS2, that treat condensed phase systems as a collection of electronically polarizable molecular fragments. The fragmentation approach to ACKS2 not only leads to a more transferable and less system-specific collection of electronic response parameters but also opens up the method to large condensed phase systems. We validate the accuracies of f-ACKS2 and scf-ACKS2 by comparing polarization energies and induced dipole moments for a number of charged hydrocarbon dimers against DFT reference calculations. Finally, we also apply both fragmented ACKS2 variants to calculate the polarization energy for electron-hole pair separation along a chain of anthracene molecules and find excellent agreement with reference DFT calculations.

Rational Design of Methodology-Independent Metal Parameters Using a Nonbonded Dummy Model.

A nonbonded dummy model for metal ions is highly imperative for the computation of complex biological systems with for instance multiple metal centers. Here we present nonbonded dummy parameters of 11 divalent metallic cations, namely, Mg(2+), V(2+), Cr(2+), Mn(2+), Fe(2+), Co(2+), Ni(2+), Zn(2+), Cd(2+), Sn(2+), and Hg(2+), that are optimized to be compatible with three widely used water models (TIP3P, SPC/E, and TIP4P-EW). The three sets of metal parameters reproduce simultaneously the solvation free energies (ΔGsol), the ion-oxygen distance in the first solvation shell (IOD), and coordination numbers (CN) in explicit water with a relative error less than 1%. The main sources of errors to ΔGsol that arise from the boundary conditions and treatment of electrostatic interactions are corrected rationally, which ensures the independence of the proposed parameters on the methodology used in the calculation. This work will be of great value for the computational study of metal-containing biological systems.

QM/MM simulation of the amide-I band in the Raman spectrum of insulin

ABSTRACTRaman spectroscopy is an effective tool to detect conformational changes and secondary structures of biological molecules. The amide-I band representing the amide carbonyl (C=O) stretching, with smaller contributions of C–N stretching and N–H bending is a signature band for protein secondary structure conformation. We have simulated the Raman spectra of insulin by a hybrid quantum-mechanics and molecular-mechanics (QM/MM) method with an aim to provide an accurate description of the amide-I band. To fulfil this aim we have considered three different QM/MM models with increasingly accurate description of the electrostatic environment for tyrosine (TYR), phenylalanine (PHE) and cystine (CYS) residues of insulin. All three models successfully describe the experimental Raman spectral features associated with the vibrational modes of the amino acid side chains. However, an accurate simulation of the amide-I band is achieved only in one of the three models, where the peptide backbone atoms together with its hydrogen bonding partners are treated with QM method. This work indicates that the accurate treatment of electrostatic interactions of the peptide backbone is crucial for correct simulation of the amide-I region, which acts as a spectral signature of proteins.

Evaluating the Strength of Salt Bridges: A Comparison of Current Biomolecular Force Fields

Recent advances in computer hardware and software have made rigorous evaluation of current biomolecular force fields using microsecond-scale simulations possible. Force fields differ in their treatment of electrostatic interactions, including the formation of salt bridges in proteins. Here we conducted an extensive evaluation of salt bridge interactions in the latest AMBER, CHARMM, and OPLS force fields, using microsecond-scale molecular dynamics simulations of amino acid analogues in explicit solvent. We focused on salt bridges between three different pairs of oppositely charged amino acids: Arg/Asp, Lys/Asp, and His(+)/Asp. Our results reveal considerable variability in the predicted KA values of the salt bridges for these force fields, as well as differences from experimental data: almost all of the force fields overestimate the strengths of the salt bridges. When amino acids are represented by side-chain analogues, the AMBER ff03 force field overestimates the KA values the least, while for complete amino acids, the AMBER ff13α force field yields the lowest KA value, most likely caused by an altered balance of side-chain/side-chain and side-chain/backbone contacts. These findings confirm the notion that the implicit incorporation of solvent polarization improves the accuracy of modeling salt bridge interactions.

Numerical simulation of photocurrent generation in bilayer organic solar cells: Comparison of master equation and kinetic Monte Carlo approaches

Numerical approaches can provide useful information about the microscopic processes underlying photocurrent generation in organic solar cells (OSCs). Among them, the Kinetic Monte Carlo (KMC) method is conceptually the simplest, but computationally the most intensive. A less demanding alternative is potentially represented by so-called Master Equation (ME) approaches, where the equations describing particle dynamics rely on the mean-field approximation and their solution is attained numerically, rather than stochastically. The description of charge separation dynamics, the treatment of electrostatic interactions and numerical stability are some of the key issues which have prevented the application of these methods to OSC modelling, despite of their successes in the study of charge transport in disordered system. Here we describe a three-dimensional ME approach to photocurrent generation in OSCs which attempts to deal with these issues. The reliability of the proposed method is tested against reference KMC simulations on bilayer heterojunction solar cells. Comparison of the current-voltage curves shows that the model well approximates the exact result for most devices. The largest deviations in current densities are mainly due to the adoption of the mean-field approximation for electrostatic interactions. The presence of deep traps, in devices characterized by strong energy disorder, may also affect result quality. Comparison of the simulation times reveals that the ME algorithm runs, on the average, one order of magnitude faster than KMC.

Atomistic Simulations of the ZnO(12̅10)/Water Interface: A Comparison between First-Principles, Tight-Binding, and Empirical Methods.

We investigate the adsorption behavior of water over the zinc oxide (12̅10) surface starting from single molecules up to bulk liquid by means of atomistic molecular dynamics simulations. We compare results obtained with density-functional theory, density-functional tight binding, and a recently developed reactive force field. The methods perform comparably up to the level of a single monolayer of adsorbed water, predicting only small differences in adsorption energies and, as a consequence, adsorption geometries. These lie within the error bars of typical quantum mechanical calculations performed with different exchange-correlation functionals. However, the discrepancies among the methods have a dramatic effect on the dissociation equilibria and the structuring of liquid water layers in contact with the surface. Especially the different treatment of electrostatic interactions via self-consistent atomic point charges appears to heavily influence the simulation outcomes. Critical comparisons with experimental studies and possibly ad hoc reparametrizations of the semiempirical functionals may thus be necessary to study phenomena such as dissolution or biomolecular adsorption at ZnO surfaces within statistically relevant time and size scales.

Calculation of Derivative Thermodynamic Hydration and Aqueous Partial Molar Properties of Ions Based on Atomistic Simulations.

The raw ionic solvation free energies calculated on the basis of atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions and treatment of electrostatic interactions used during these simulations. However, as shown recently [Kastenholz, M. A.; Hünenberger, P. H. J. Chem. Phys.2006, 124, 224501 and Reif, M. M.; Hünenberger, P. H. J. Chem. Phys.2011, 134, 144104], the application of an appropriate correction scheme allows for a conversion of the methodology-dependent raw data into methodology-independent results. In this work, methodology-independent derivative thermodynamic hydration and aqueous partial molar properties are calculated for the Na(+) and Cl(-) ions at P° = 1 bar and T(-) = 298.15 K, based on the SPC water model and on ion-solvent Lennard-Jones interaction coefficients previously reoptimized against experimental hydration free energies. The hydration parameters considered are the hydration free energy and enthalpy. The aqueous partial molar parameters considered are the partial molar entropy, volume, heat capacity, volume-compressibility, and volume-expansivity. Two alternative calculation methods are employed to access these properties. Method I relies on the difference in average volume and energy between two aqueous systems involving the same number of water molecules, either in the absence or in the presence of the ion, along with variations of these differences corresponding to finite pressure or/and temperature changes. Method II relies on the calculation of the hydration free energy of the ion, along with variations of this free energy corresponding to finite pressure or/and temperature changes. Both methods are used considering two distinct variants in the application of the correction scheme. In variant A, the raw values from the simulations are corrected after the application of finite difference in pressure or/and temperature, based on correction terms specifically designed for derivative parameters at P° and T(-). In variant B, these raw values are corrected prior to differentiation, based on corresponding correction terms appropriate for the different simulation pressures P and temperatures T. The results corresponding to the different calculation schemes show that, except for the hydration free energy itself, accurate methodological independence and quantitative agreement with even the most reliable experimental parameters (ion-pair properties) are not yet reached. Nevertheless, approximate internal consistency and qualitative agreement with experimental results can be achieved, but only when an appropriate correction scheme is applied, along with a careful consideration of standard-state issues. In this sense, the main merit of the present study is to set a clear framework for these types of calculations and to point toward directions for future improvements, with the ultimate goal of reaching a consistent and quantitative description of single-ion hydration thermodynamics in molecular dynamics simulations.

On simulations of complex interfaces: Molecular dynamics simulations of stationary phases

Methodological considerations for molecular dynamics simulations of complex interfaces are presented in this article. A slab geometry is examined in the context of stationary phases where selectivity occurs predominantly in pores within silica beads. Specifically, we examine the Whelk-O1 interface with n-hexane/2-propanol, the TMA-(Pro)(2)-N(CH(3))-tether interface with n-hexane/2-propanol, and the C(18)H(37)Si interface with water/methanol. The following methodological issues are considered in detail: The assessment of solvent density within the confined region and excluded volume of the interface; the structural equilibration of surface-bound moieties; solvent equilibration for binary mixtures; surface size effects, and periodic boundary conditions; the treatment of electrostatic interactions; and the impact of pore size.

On the Treatment of Electrostatic Interactions of Non-spherical Molecules in Equation of State Models

For predictive applications, equation of state (EOS) models have to describe all relevant physical interactions accurately. In this contribution, the vapor–liquid equilibria of various dipolar two-center Lennard-Jones model molecules are determined by molecular simulation, as function of molecular elongation and deflection angle of the dipole. It is shown that present PC-SAFT-based EOS models require additional adjustable parameters in order to describe the orientational effects of the dipole-moment. We present extensions of the model to avoid the additional parameters and apply the extended equations to model systems and real molecules.

Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant.

The chemical composition of small organic molecules is often very similar to amino acid side chains or the bases in nucleic acids, and hence there is no a priori reason why a molecular mechanics force field could not describe both organic liquids and biomolecules with a single parameter set. Here, we devise a benchmark for force fields in order to test the ability of existing force fields to reproduce some key properties of organic liquids, namely, the density, enthalpy of vaporization, the surface tension, the heat capacity at constant volume and pressure, the isothermal compressibility, the volumetric expansion coefficient, and the static dielectric constant. Well over 1200 experimental measurements were used for comparison to the simulations of 146 organic liquids. Novel polynomial interpolations of the dielectric constant (32 molecules), heat capacity at constant pressure (three molecules), and the isothermal compressibility (53 molecules) as a function of the temperature have been made, based on experimental data, in order to be able to compare simulation results to them. To compute the heat capacities, we applied the two phase thermodynamics method (Lin et al. J. Chem. Phys.2003, 119, 11792), which allows one to compute thermodynamic properties on the basis of the density of states as derived from the velocity autocorrelation function. The method is implemented in a new utility within the GROMACS molecular simulation package, named g_dos, and a detailed exposé of the underlying equations is presented. The purpose of this work is to establish the state of the art of two popular force fields, OPLS/AA (all-atom optimized potential for liquid simulation) and GAFF (generalized Amber force field), to find common bottlenecks, i.e., particularly difficult molecules, and to serve as a reference point for future force field development. To make for a fair playing field, all molecules were evaluated with the same parameter settings, such as thermostats and barostats, treatment of electrostatic interactions, and system size (1000 molecules). The densities and enthalpy of vaporization from an independent data set based on simulations using the CHARMM General Force Field (CGenFF) presented by Vanommeslaeghe et al. (J. Comput. Chem.2010, 31, 671) are included for comparison. We find that, overall, the OPLS/AA force field performs somewhat better than GAFF, but there are significant issues with reproduction of the surface tension and dielectric constants for both force fields.

Computation of methodology-independent single-ion solvation properties from molecular simulations. IV. Optimized Lennard-Jones interaction parameter sets for the alkali and halide ions in water

The raw single-ion solvation free energies computed from atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions and treatment of electrostatic interactions used during these simulations. However, as shown recently [M. A. Kastenholz and P. H. Hünenberger, J. Chem. Phys. 124, 224501 (2006); M. M. Reif and P. H. Hünenberger, J. Chem. Phys. 134, 144103 (2010)], the application of appropriate correction terms permits to obtain methodology-independent results. The corrected values are then exclusively characteristic of the underlying molecular model including in particular the ion-solvent van der Waals interaction parameters, determining the effective ion size and the magnitude of its dispersion interactions. In the present study, the comparison of calculated (corrected) hydration free energies with experimental data (along with the consideration of ionic polarizabilities) is used to calibrate new sets of ion-solvent van der Waals (Lennard-Jones) interaction parameters for the alkali (Li(+), Na(+), K(+), Rb(+), Cs(+)) and halide (F(-), Cl(-), Br(-), I(-)) ions along with either the SPC or the SPC/E water models. The experimental dataset is defined by conventional single-ion hydration free energies [Tissandier et al., J. Phys. Chem. A 102, 7787 (1998); Fawcett, J. Phys. Chem. B 103, 11181] along with three plausible choices for the (experimentally elusive) value of the absolute (intrinsic) hydration free energy of the proton, namely, ΔG(hyd)(⊖)[H(+)] = -1100, -1075 or -1050 kJ mol(-1), resulting in three sets L, M, and H for the SPC water model and three sets L(E), M(E), and H(E) for the SPC/E water model (alternative sets can easily be interpolated to intermediate ΔG(hyd)(⊖)[H(+)] values). The residual sensitivity of the calculated (corrected) hydration free energies on the volume-pressure boundary conditions and on the effective ionic radius entering into the calculation of the correction terms is also evaluated and found to be very limited. Ultimately, it is expected that comparison with other experimental ionic properties (e.g., derivative single-ion solvation properties, as well as data concerning ionic crystals, melts, solutions at finite concentrations, or nonaqueous solutions) will permit to validate one specific set and thus, the associated ΔG(hyd)(⊖)[H(+)] value (atomistic consistency assumption). Preliminary results (first-peak positions in the ion-water radial distribution functions, partial molar volumes of ionic salts in water, and structural properties of ionic crystals) support a value of ΔG(hyd)(⊖)[H(+)] close to -1100 kJ·mol(-1).

Computation of methodology-independent single-ion solvation properties from molecular simulations. III. Correction terms for the solvation free energies, enthalpies, entropies, heat capacities, volumes, compressibilities, and expansivities of solvated ions

The raw single-ion solvation free energies computed from atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions (finite or periodic system, system or box size) and treatment of electrostatic interactions (Coulombic, lattice-sum, or cutoff-based) used during these simulations. However, as shown by Kastenholz and Hünenberger [J. Chem. Phys. 124, 224501 (2006)], correction terms can be derived for the effects of: (A) an incorrect solvent polarization around the ion and an incomplete or/and inexact interaction of the ion with the polarized solvent due to the use of an approximate (not strictly Coulombic) electrostatic scheme; (B) the finite-size or artificial periodicity of the simulated system; (C) an improper summation scheme to evaluate the potential at the ion site, and the possible presence of a polarized air-liquid interface or of a constraint of vanishing average electrostatic potential in the simulated system; and (D) an inaccurate dielectric permittivity of the employed solvent model. Comparison with standard experimental data also requires the inclusion of appropriate cavity-formation and standard-state correction terms. In the present study, this correction scheme is extended by: (i) providing simple approximate analytical expressions (empirically-fitted) for the correction terms that were evaluated numerically in the above scheme (continuum-electrostatics calculations); (ii) providing correction terms for derivative thermodynamic single-ion solvation properties (and corresponding partial molar variables in solution), namely, the enthalpy, entropy, isobaric heat capacity, volume, isothermal compressibility, and isobaric expansivity (including appropriate standard-state correction terms). The ability of the correction scheme to produce methodology-independent single-ion solvation free energies based on atomistic simulations is tested in the case of Na(+) hydration, and the nature and magnitude of the correction terms for derivative thermodynamic properties is assessed numerically.

Methodological assessment of kinetic Monte Carlo simulations of organic photovoltaic devices: The treatment of electrostatic interactions

The kinetic Monte Carlo (KMC) method provides a versatile tool to investigate the mechanisms underlying photocurrent generation in nanostructured organic solar cells. Currently available algorithms can already support the development of more cost-efficient photovoltaic devices, but so far no attempt has been made to test the validity of some fundamental model assumptions and their impact on the simulation result. A meaningful example is given by the treatment of the electrostatic interactions. In most KMC models, electrostatic interactions are approximated by means of cutoff based potentials, irrespective of the long-range nature of the Coulomb interaction. In this paper, the reliability of such approximation is tested against the exact Ewald sum. The results under short-circuit and flat-band conditions show that use of cutoff-based potentials tends to underestimate real device performance, in terms of internal quantum efficiency and current density. Together with this important finding, we formalize other methodological aspects which have been scarcely discussed in the literature.

Treatment of electrostatic interactions in simulations of the interface between transition metal surfaces and organic matter: the discrete classical model and the reaction-field dielectric method

A more efficient molecular dynamics model is devised to study the interface between platinum(111) and liquid isopropanol. The faster reaction-field dielectric method (RFD) is investigated as the alternative to the very accurate discrete classical model (DCM) (M. W. Finnis, Surf. Sci., 1991, 241, 61) to account for the electrostatic interactions at the interface. Differences in static and dynamic properties of the interface are examined. The mass-density and number-density distributions showed only a little sensitivity towards the choice of the method. The orientational structure revealed notable but still small quantitative changes for highly-polarized parts of the isopropanol molecules if the RFD is used instead of the DCM. The general qualitative picture is still preserved. The simulations of the platinum(111)/isopropanol interface with the DCM or the Ewald summation were carried out for longer times than in our previous contribution (K. B. Tarmyshov and F. Müller-Plathe, J. Chem. Phys., 2007, 126, 74702). The analysis disclosed that the influence of chemical or electrostatic interactions at the interface on the relative orientation of the neighbouring hydroxyl groups, which are adsorbed on the surface, is notably smaller than reported before. Overall, the reaction-field dielectric method can be used to study the properties of the interface between transition metals and organic matter, if there appears a chemical bond of 10 kT or higher between the surface and adsorbate molecules or their parts.

Comparison of different potential models in the calculations of the infrared and polarized Raman spectra and the transient infrared absorption anisotropy of the OH stretching mode of liquid water

The infrared (IR) and polarized Raman spectra of the OH stretching mode of neat liquid H 2O, as well as the time dependence of the transient IR (TRIR) absorption anisotropy of the same vibrational mode of both neat liquid H 2O and a dilute liquid mixture of H 2O/HDO/D 2O, are calculated by using four widely-used potential models (SPC/E, TIP3P, TIP4P, and TIP5P) to see the effects of the difference in the treatment of electrostatic interactions on the calculated spectral profiles. In the calculations, three important factors for simulating vibrational spectra of coupled-oscillators systems in the liquid state (frequency shifts and fluctuations of individual molecules, resonant vibrational couplings between molecules, and liquid dynamics) are taken into account with a time-domain computational method. It is shown that, although the calculations with all the four potential models commonly reproduce the main observed spectral features, there are also some significant differences. Those differences are partly a matter of the calculated magnitudes of the electric fields operating on OH bonds from hydrogen-bonded molecules, and are discussed in relation to the differences in the calculated radial distribution function. It is also shown that the time constant of the long-time decay (in the region of t ≥ 1.5 ps) of the TRIR anisotropy in a dilute liquid mixture of H 2O/HDO/D 2O is calculated as 2.7–2.8 ps by using the SPC/E and TIP5P potentials, in reasonable agreement with the experimental results in the literature. An improved and self-consistent treatment of the electric field-induced enhancement of dipole derivatives, which is another important factor for simulating the OH stretching band of liquid water, is also proposed.

Dynamical motions of lipids and a finite size effect in simulations of bilayers.

Molecular dynamics (MD) simulations of dipalmitoylphosphatidylcholine bilayers composed of 72 and 288 lipids are used to examine system size dependence on dynamical properties associated with the particle mesh Ewald (PME) treatment of electrostatic interactions. The lateral diffusion constant Dl is 2.92 x 10(-7) and 0.95 x 10(-7) cm2/s for 72 and 288 lipids, respectively. This dramatic finite size effect originates from the correlation length of lipid diffusion, which extends to next-nearest neighbors in the 288 lipid system. Consequently, diffusional events in smaller systems can propagate across the boundaries of the periodic box. The internal dynamics of lipids calculated from the PME simulations are independent of the system size. Specifically, reorientational correlation functions for the slowly relaxing phosphorus-glycerol hydrogen, phosphorus-nitrogen vectors, and more rapidly relaxing CH vectors in the aliphatic chains are equivalent for the 72 and 288 lipid simulations. A third MD simulation of a bilayer with 72 lipids using spherical force-shift electrostatic cutoffs resulted in interdigitated chains, thereby rendering this cutoff method inappropriate.

Computation of methodology-independent ionic solvation free energies from molecular simulations. II. The hydration free energy of the sodium cation

The raw ionic solvation free energies computed from atomistic (explicit-solvent) simulations are extremely sensitive to the boundary conditions (finite or periodic system, system shape, and size) and treatment of electrostatic interactions (Coulombic, lattice sum, or cutoff based) used during these simulations. In the present article, it is shown that correction terms can be derived for the effect of (A) an incorrect solvent polarization around the ion due to the use of an approximate (not strictly Coulombic) electrostatic scheme; (B) the finite size or artificial periodicity of the simulated system; (C) an improper summation scheme to evaluate the potential at the ion site and the possible presence of a liquid-vacuum interface in the simulated system. Taking the hydration free energy of the sodium cation as a test case, it is shown that the raw solvation free energies obtained using seven different types of boundary conditions and electrostatic schemes commonly used in explicit-solvent simulations (for a total of 72 simulations differing in the corresponding simulation parameters) can be corrected so as to obtain a consistent value for this quantity.

On the Ewald Artifacts in Computer Simulations. The Test-Case of the Octaalanine Peptide With Charged Termini

The treatment of electrostatic interactions in molecular simulations is of fundamental importance. Ewald and related methods are being increasingly used to the detriment of cutoff schemes, which are known to produce several artifacts. A potential drawback of the Ewald method is the spatial periodicity that is imposed to the system, which could produce artifacts when applied in the simulation of liquids. In this work we analyze the octaalanine peptide with charged termini in explicit solvent, for which severe effects due to the use of Ewald sums were predicted using continuum electrostatics. Molecular Dynamics simulations for a total of 158 nanoseconds were performed in cells of different sizes. From the comparison of the results of different system sizes, no significant periodicity-induced artifacts were observed. It is argued that in current biomolecular simulations, the incomplete sampling is likely to affect the results to a larger extent than the artifacts induced by the use of Ewald sums.

The Influence of Different Treatments of Electrostatic Interactions on the Thermodynamics of Folding of Peptides

Replica exchange molecular dynamics simulations were performed to investigate the effects of different electrostatic treatments on the structure and thermodynamics of a small beta-hairpin forming peptide. Three different electrostatic schemes were considered: regular cutoffs, generalized reaction field (GRF), and particle mesh Ewald (PME), with the peptide modeled using OPLS/AA all-atom force field with explicit TIP3P water. Both the GRF and PME methods yielded results consistent with experiment, with free energy surfaces displaying a single minimum corresponding to the native beta-hairpin structure. In contrast, use of straight cutoffs led to the population of an additional local minimum corresponding to nonhairpin conformations that compete with the formation of the native beta-hairpin at low temperatures. This extra minimum would not be apparent in conventional constant-temperature molecular dynamics simulations run for a few nanoseconds. This result points to the critical need of careful sampling of conformational space to assess the quality of different numerical treatments of long-range forces. While differences emerged in the nature of the unfolded states populated using PME and GRF approaches, simulations on the beta-hairpin forming peptide and on two additional control peptides indicate that the GRF treatment of electrostatics offers a satisfactory compromise between accuracy and computational speed for the identification of low-energy conformations. A GRF-based approach emerges as a viable means for treating larger biological systems that would be prohibitively costly to simulate using PME methods.

Lipid Bilayers Driven to a Wrong Lane in Molecular Dynamics Simulations by Subtle Changes in Long-Range Electrostatic Interactions

We provide compelling evidence that different treatments of electrostatic interactions in molecular dynamics simulations may dramatically affect dynamic properties of lipid bilayers. To this end, we consider a fully hydrated pure dipalmitoylphosphatidylcholine bilayer through 50-ns molecular dynamics simulations and study various dynamic properties of individual lipids in a membrane, including the velocity autocorrelation function, the lateral and rotational diffusion coefficients, and the autocorrelation function for the area per molecule. We compare the results based on the Particle-Mesh Ewald (PME) and reaction field (RF) techniques with those obtained by an approach where the electrostatic interactions are truncated at rcut = 1.8, 2.0, and 2.5 nm. We find that the overall performance of PME is very good; its results are consistent with the expected behavior. The RF method performs rather well, too, despite certain inherent problems and the fact that its results differ from those obtained by PME. Never...