Ewald Summation Research Articles

Determination of Ionic Hydration Free Energies with Grand Canonical Monte Carlo/Molecular Dynamics Simulations in Explicit Water.

Grand canonical Monte Carlo (GCMC) simulations of ionic solutions with explicit solvent models are known to be challenging. One challenge arises from the treatment of long-range electrostatics and finite-box size in Monte Carlo simulations when periodic boundary condition and Ewald summation methods are used. Another challenge is that constant excess chemical potential GCMC simulations for charged solutes suffer from inadequate insertion and deletion acceptance ratios. In this work, we address those problems by implementing an oscillating excess chemical potential GCMC algorithm with smooth particle mesh Ewald and finite-box-size corrections to treat the long-range electrostatics. The developed GCMC simulation program was combined with GROMACS to perform GCMC/MD simulations of ionic solutions individually containing Li+, Na+, K+, Rb+, Cs+, F-, Cl-, Br-, I-, Ca2+, and Mg2+, respectively. Our simulation results show that the combined GCMC/MD approach can approximate the ionic hydration free energies with proper treatment of long-range electrostatics. Our developed simulation approach can open up new avenues for simulating complex chemical and biomolecular systems and for drug discovery.

A physically grounded damped dispersion model with particle mesh Ewald summation.

Accurate modeling of dispersion is critical to the goal of predictive biomolecular simulations. To achieve this accuracy, a model must be able to correctly capture both the short-range and asymptotic behavior of dispersion interactions. We present here a damped dispersion model based on the overlap of charge densities that correctly captures both regimes. The overlap damped dispersion model represents a classical physical interpretation of dispersion: the interaction between the instantaneous induced dipoles of two distinct charge distributions. This model is shown to be an excellent fit with symmetry adapted perturbation theory dispersion energy calculations, yielding an RMS error on the S101x7 database of 0.5 kcal/mol. Moreover, the damping function used in this model is wholly derived and parameterized from the electrostatic dipole-dipole interaction, making it not only physically grounded but transferable as well.

Fast Ewald summation for Green\u2019s functions of Stokes flow in a half-space

Recently, Gimbutas et al. (J Fluid Mech, 2015. https://doi.org/10.1017/jfm.2015.302) derived an elegant representation for the Green’s functions of Stokes flow in a half-space. We present a fast summation method for sums involving these half-space Green’s functions (stokeslets, stresslets and rotlets) that consolidates and builds on the work by Klinteberg et al. (Res Math Sci 4(1):1, 2017. https://doi.org/10.1186/s40687-016-0092-7) for the corresponding free-space Green’s functions. The fast method is based on two main ingredients: The Ewald decomposition and subsequent use of FFTs. The Ewald decomposition recasts the sum into a sum of two exponentially decaying series: one in real space (short-range interactions) and one in Fourier space (long-range interactions) with the convergence of each series controlled by a common parameter. The evaluation of short-range interactions is accelerated by restricting computations to neighbours within a specified distance, while the use of FFTs accelerates the computations in Fourier space thus accelerating the overall sum. We demonstrate that while the method incurs extra costs for the half-space in comparison with the free-space evaluation, greater computational savings is also achieved when compared to their respective direct sums.

A fast and accurate computational method for the linear-combination-based isotropic periodic sum

An isotropic periodic sum (IPS) is a powerful technique to reasonably calculate intermolecular interactions for wide range of molecular systems under periodic boundary conditions. A linear-combination-based IPS (LIPS) has been developed to attain computational accuracy close to an exact lattice sum, such as the Ewald sum. The algorithm of the original LIPS method has a high computational cost because it needs long-range interaction calculations in real space. This becomes a performance bottleneck for long-time molecular simulations. In this work, the combination of an LIPS and fast Fourier transform (FFT) was developed, and evaluated on homogeneous and heterogeneous molecular systems. This combinational approach of LIPS/FFT attained computational efficiency close to that of a smooth particle mesh Ewald while maintaining the same high accuracy as the original LIPS. We concluded that LIPS/FFT has great potential to extend the capability of IPS techniques for the fast and accurate computation of many types of molecular systems.

Communication: Computing the Helmholtz capacitance of charged insulator-electrolyte interfaces from the supercell polarization.

Supercell modeling of an electrical double layer (EDL) at electrified solid-electrolyte interfaces is a challenge. The net polarization of EDLs arising from the fixed chemical composition setup leads to uncompensated EDLs under periodic boundary condition and convolutes the calculation of the Helmholtz capacitance [C. Zhang and M. Sprik, Phys. Rev. B 94, 245309 (2016)]. Here we provide a new formula based on the supercell polarization at zero electric field Ē = 0 (i.e., standard Ewald boundary condition) to calculate the Helmholtz capacitance of charged insulator-electrolyte interfaces and validate it using atomistic simulations. Results are shown to be independent of the supercell size. This formula gives a shortcut to compute the Helmholtz capacitance without locating the zero net charge state of EDL and applies directly to any standard molecular dynamics code where the electrostatic interactions are treated by the Ewald summation or its variants.

Analysis of interleukin 23 and 7G10 interactions for computational design of lead antibodies against immune-mediated inflammatory diseases

Wealth of structural data on theurapeutic targets in complex with monoclonal antibodies (mAbs) and advances in molecular modeling algorithms present exciting opportunities in the field of novel biologic design. Interleukin 23 (IL23), a well-known drug target for autoimmune diseases, in complex with mAb 7G10 offers prospect to design potent lead antibodies by traversing the complete epitope–paratope interface. Herein, key interactions aiding antibody-based neutralization in IL23–7G10 complex are resolute through PyMOL, LigPlot+, Antibody i-Patch, DiscoTope and FoldX. Six amino acids Ser31, Val33, Asn55, Lys59 in heavy chain and His34, Ser93 in light chain are subjected to in silico mutagenesis with residues Met, Trp, Ile, Leu and Arg. A set of 431 mutant macromolecules are outlined. Binding affinities of these molecules with IL23 are estimated through protein–protein docking by employing ZDOCK, ClusPro and RosettaDock. Subsequently, the macromolecules revealed comparable result with 7G10 are cross validated through binding free-energy calculations by applying Molecular Mechanics/Poisson Boltzman Surface Area method in CHARMM. Thirty nine designed theoretical antibodies showed improved outcome in all evaluations; from these, top 10 molecules showed at least nine unit better binding affinity compared to the known mAb. These molecules have the potential to act as lead antibodies. Subsequent molecular dynamics simulations too favored prospective of best ranked molecule to have therapeutic implications in autoimmune and inflammatory diseases.Abbreviations: IL23: interleukin 23; IL17: interleukin17; Ab: antibody; Ag: antigen; mAbs: monoclonal antibodies; STAT3: signal transducer and activator of transcription 3; STAT4: signal transducer and activator of transcription 4; PDB: protein databank; MM/PBSA: molecular mechanics Poisson-Boltzmann surface area; Ag–Ab: antigen– antibody complex; SPC/E: extended simple point charge; SD: steepest descents; PME: particle mesh ewald; dG: binding free energies; Fv: variable fragment

Replicate or alternate? Isotropic boundary conditions and implications for Ewald summation

Periodic boundary conditions (PBC) are well suited to describe repetitive structures in space, yet flawed for isotropic setups in which it is commonly used. We present an approach based on exploiting the congruency of a unit cell which combines the periodic nature of PBC with isotropy of an alternating cell. By applying these isotropic periodic boundary conditions (IPBC) to the standard Ewald summation for electrostatic interactions, we show a marked reduced artificial ordering and a need for significantly less wave vectors to an otherwise equivalent PBC system. The methodology can be trivially implemented in existing molecular dynamics or Metropolis Monte Carlo Ewald summation codes, and is applicable also for particle mesh Ewald summation.

Magnetic response of Fe nanoparticles embedded in artificial SiO2 opals

Applying the Stöber method, silica spheres were obtained and arranged in an fcc lattice. Metallic Fe was introduced into the opal voids. Magnetization curves of the synthesized metallic nanoparticles and of the artificial opals with infiltrated nanoparticles were measured by using the vibrating sample magnetometry (VSM) technique with an applied magnetic field oriented along the [111] direction of the opal matrix. Nonlinear Landau-Lifshitz-Gilbert (LLG) equation was numerically solved, in order to reproduce the experimental results. In solving LLG equation we have used Ewald summations to take into account the dipolar interaction between the spatially arranged magnetic nanoparticles. The theory reproduces the general features of the measured magnetization curves. In particular, the S-like shape of the magnetization curves produced by the ferromagnetic response of Fe nanoparticles and the diamagnetic one of the SiO2 opal matrix is well described. The observed small areas of the magnetization hysteresis loops are explained as a consequence of the random distribution of the magnetic anisotropy axes of the Fe nanoparticles. The effect of the dipolar interaction between the magnetic nanoparticles is predicted to be important at magnitudes of the external applied field, larger than the coercive magnetic field.

Bicanonical ensemble Monte Carlo simulation of water condensation in the field of crystal lattice defects

The free energy, formation work and entropy dependences for water condensate formed from the vapor over the defect-free and containing surface defects basal face of β-AgI crystal at initial stage of condensation at temperatures of 260 K and 400 K are calculated using the bicanonical statistical ensemble method, with Ewald summation for long-range electrostatic and polarization interactions with the substrate. The effect of surface defects in form of rectangular “towers” as a part of regular structure on the stability of condensed phase embryos is investigated. In contrast to small-scale structures, relatively larger coarse-grained nanostructure of crystal surface demonstrates an unconditional advantage in the ability to stimulate condensation compared to defect-free surface. The formation of condensed phase embryos on the surface with multiple defects begins at vapor pressures 5–6 orders of magnitude lower than that on corresponding defect-free surface, and this effect is resistant to temperature variations. The condensate on the surface of the crystal is thermodynamically stable, both on defect-free and nanostructured surfaces, with the exception of short initial stage of the monomolecular film on the defect-free surface.

Molecular Dynamics Simulation of Mahkota Dewa (Phaleria Macrocarpa) Extract in Subcritical Water Extraction Process

Mahkota Dewa (Phaleria Macrocarpa), a good source of saponin, flavanoid, polyphenol, alkaloid, and mangiferin has an extensive range of medicinal effects. The intermolecular interactions between solute and solvents such as hydrogen bonding considered as an important factor that affect the extraction of bioactive compounds. In this work, molecular dynamics simulation was performed to elucidate the hydrogen bonding exists between Mahkota Dewa extracts and water during subcritical extraction process. A bioactive compound in the Mahkota Dewa extract, namely mangiferin was selected as a model compound. The simulation was performed at 373 K and 4.0 MPa using COMPASS force field and Ewald summation method available in Material Studio 7.0 simulation package. The radial distribution functions (RDF) between mangiferin and water signify the presence of hydrogen bonding in the extraction process. The simulation of the binary mixture of mangiferin:water shows that strong hydrogen bonding was formed. It is suggested that, the intermolecular interaction between OH2O••HMR4(OH1) has been identified to be responsible for the mangiferin extraction process.

Towards Electronic Detection of DNA Conformational Transition

The conformation of DNA plays an important role in its chemical and biological function as well as for its potential applications in nanotechnology1. Different conformations of DNA exhibit distinct values for the conductance2. Such a distinct conductance can also be used to detect the conformation. DNA mainly adopts A and B conformations at lower and higher hydration levels, respectively1,2. Experimental and theoretical studies suggest that A to B transformation is smooth and continuous1, whereas B to A transition is sequence specific, cooperative, and reversible3,4. To the best of our knowledge, there is no experimental method that can precisely evaluate the electronic structure of DNA during its conformation transition. A recent theoretical study discusses the change in electronic properties during the transformation1. A recent work found an unexpected result that A-DNA conducts better than B-DNA. To provide insight into the time evolution of the electronic properties during change in conformation, we present a computational investigation of the electrical conductance of 3’-CCCGCGCCC-5’ DNA during various time step of the MD trajectory as the molecule undergoes A to B transition. The initial structure of A-DNA is built using nucleic acid builder (NAB) in AMBER15. For DNA, the parmbsc1 force field is used and SPC/E force field is used for water and ions. The negative charge on DNA backbone is neutralized with sodium (Na) ions. Additionally, Na+ and Cl- ions are added to achieve 150mM salt concentration. For A to B DNA transformation, A-DNA is placed in water. Further minimization and equilibration is performed according to ABC Minimization5. 50 ns molecular dynamics (MD) trajectory is generated with particle mesh Ewald (PME) algorithm and periodic boundary condition (PBC) under constant pressure (isobaric-isothermal ensemble). The transformed B-DNA at the end of 50ns simulation is then placed in 85%ethanol+15%water with 150mM of salt concentration. Packmol is used to model the 85%ethanol+15%water solvent. The temperature is increased to 325K and the simulation is continued for another 48ns. The snapshots during A to B transformations (at 1ps, 50 ps, 100ps, 125ps, 150ps, 200ps and 50000ps, respectively) are considered for electronic structure calculations, where the DNA encounters major structural changes (also done for B to A transformation). Density Functional Theory (DFT) calculations are performed on these structures using GAUSSIAN 09, with the B3LYP functional and 6-31G(d,p) basis set. The electrical conductance of DNA is determined using Landauer formalism with Non-equilibrium Green's function approach2. Figure 1 (a) shows the sugar pucker in the range of 0o-36o (representing A-DNA) in the initial phase of simulation that increases step by step to 140o-180o (representing B-DNA) and remains in that range till 50ns. We also find that the rise per base pairs increases from ~2.6Å to ~3.4 Å and twist increases from 30-34 Å to 35-38 Å, which confirm the transformation. We could reverse this transformation when the B-DNA is placed in 85%ethanol+15%water at 325K (results not presented here). A-DNA at 1ps has HOMO at -5.04eV and LUMO at -1.04eV. During the transformation, mean value of HOMO energy is -5.06eV and that of LUMO is -1.21eV, indicating significant change in the LUMO energy (~20%). We next look at the charge transport during the transformation, i.e. for the frames at 1 ps, 50 ps, 100 ps, 125 ps, 150ps, 200ps and 50000ps. Overall conductance is seen to decrease (increase) with the A(B) to B(A) transformation (Figure 1 (b) for A to B transformation). It is to be noted that this is the same DNA sequence investigated experimentally2 for which the A-conformation exhibits higher conductance than B-conformation. Our simulations clearly verify these experimental findings and provide further information on the evolution of conductance with time. The fluctuations in between could be related to the spatial configurations of counterions, which can alter the conductance of DNA6. References 1 H. Wang, T.E. Cheatham, P.M. Gannett, and J.P. Lewis, Soft Matter 5, 685 (2009). 2 J.M. Artés, Y. Li, J. Qi, M.P. Anantram, and J. Hihath, Nat. Commun. 6, (2015). 3 J.T. Waters, X.-J. Lu, R. Galindo-Murillo, J.C. Gumbart, H.D. Kim, T.E. Cheatham, and S.C. Harvey, J. Phys. Chem. B 120, 8449 (2016). 4 M. Kulkarni and A. Mukherjee, Prog. Biophys. Mol. Biol. 128, 63 (2017). 5 M. Pasi, J.H. Maddocks, D. Beveridge, T.C. Bishop, D.A. Case, I.I.I.T. Cheatham, P.D. Dans, B. Jayaram, F. Lankas, C. Laughton, J. Mitchell, R. Osman, M. Orozco, A. Pérez, D. Petkevičiūtė, N. Spackova, J. Sponer, K. Zakrzewska, and R. Lavery, Nucleic Acids Res. 42, 12272 (2014). 6 R.N. Barnett, C.L. Cleveland, A. Joy, U. Landman, and G.B. Schuster, Science (80-. ). 294, 567 (2001). Figure 1

Critical test of isotropic periodic sum techniques with group-based cut-off schemes

Truncation is still chosen for many long-range intermolecular interaction calculations to efficiently compute free-boundary systems, macromolecular systems and net-charge molecular systems, for example. Advanced truncation methods have been developed for long-range intermolecular interactions. Every truncation method can be implemented as one of two basic cut-off schemes, namely either an atom-based or a group-based cut-off scheme. The former computes interactions of “atoms” inside the cut-off radius, whereas the latter computes interactions of “molecules” inside the cut-off radius. In this work, the effect of group-based cut-off is investigated for isotropic periodic sum (IPS) techniques, which are promising cut-off treatments to attain advanced accuracy for many types of molecular system. The effect of group-based cut-off is clearly different from that of atom-based cut-off, and severe artefacts are observed in some cases. However, no severe discrepancy from the Ewald sum is observed with the extended IPS techniques.

Fast evaluation of the fundamental singularities of two-dimensional doubly periodic Stokes flow

Analytical representations of the flows associated with doubly periodic arrangements of point singularities of two-dimensional Stokes flow are derived. The analysis makes use of a conformal mapping from a concentric annulus to a rectangle. A natural mathematical object known as the Schottky–Klein prime function is used to derive the solutions avoiding the theory of elliptic functions. The new expressions have the advantage of being immediately amenable to fast evaluation without the need for Ewald or other fast summation methods or acceleration techniques for lattice sums.

Modification of the Wolf Method and Evaluation for Molecular Simulation of Vapor–Liquid Equilibria

A modification of the Wolf method [ Wolf et al. J. Chem. Phys. 1999 , 110 , 8254 - 8282 ], a spherically truncated pairwise summation to evaluate electrostatic interactions efficiently, is proposed. This method achieves better results for the energy, and the approach is used for determining phase equilibria in the grand canonical ensemble. To assess optimal parameters for the Wolf summation we propose a simple iterative approach using the Ewald summation as a reference. We show that phase equilibrium properties for pure components as well as mixtures can be calculated accurately using the Wolf summation. Our study considers molecular fluids with partial charges but without net charge.

Rapid Convergence of Energy and Free Energy Profiles with Quantum Mechanical Size in Quantum Mechanical-Molecular Mechanical Simulations of Proton Transfer in DNA.

In recent years, a number of quantum mechanical-molecular mechanical (QM/MM) enzyme studies have investigated the dependence of reaction energetics on the size of the QM region using energy and free energy calculations. In this study, we revisit the question of QM region size dependence in QM/MM simulations within the context of energy and free energy calculations using a proton transfer in a DNA base pair as a test case. In the simulations, the QM region was treated with a dispersion-corrected AM1/d-PhoT Hamiltonian, which was developed to accurately describe phosphoryl and proton transfer reactions, in conjunction with an electrostatic embedding scheme using the particle-mesh Ewald summation method. With this rigorous QM/MM potential, we performed rather extensive QM/MM sampling, and found that the free energy reaction profiles converge rapidly with respect to the QM region size within ca. ±1 kcal/mol. This finding suggests that the strategy of QM/MM simulations with reasonably sized and selected QM regions, which has been employed for over four decades, is a valid approach for modeling complex biomolecular systems. We point to possible causes for the sensitivity of the energy and free energy calculations to the size of the QM region, and potential implications.

Application of Screening Functions as Cutoff-Based Alternatives to Ewald Summation in Molecular Dynamics Simulations Using Polarizable Force Fields.

The range-dependent screening of the charge-charge, charge-induced dipole, and induced dipole-induced dipole interactions was examined for a variety of liquids modeled using polarizable force fields. A cutoff-based method for calculation of the electrostatic interactions in molecular dynamics (MD) is presented as an alternative to Ewald-type summation for simulations of the disordered materials modeled using many-body polarizable force fields with permanent charges and induced point dipoles. The proposed approach was tested on bulk water, room-temperature ionic liquids, and solutions of ions in polar solvents. The smooth, short-range, and atom-type independent screening functions for interactions between the charges and induced dipoles were obtained using the force matching approach. An excellent agreement for both the magnitude and directionality of forces, structural and dynamic properties, was found in MD simulations utilizing the developed screening functions, compared to those with Ewald summation. While similar in shape and range, the charge-charge screening functions were somewhat dependent on the material chemistry. On the other hand, the charge-induced dipole and induced dipole-induced dipole screening functions were found to be nearly universal for the tested materials.

Dissipative particle dynamics: Effects of thermostating schemes on nano-colloid electrophoresis

A novel fully explicit approach using dissipative particle dynamics (DPD) method is introduced in the present study to model the electrophoretic transport of nano-colloids in an electrolyte solution. Slater type charge smearing function included in 3D Ewald summation method is employed to treat electrostatic interaction. Performance of various thermostats are challenged to control the system temperature and study the dynamic response of colloidal electrophoretic mobility under practical ranges of external electric field (0.072<E<0.361 v/nm) covering linear to non-linear response regime, and ionic salt concentration (0.049 <SC<0.69 [M]) covering weak to strong Debye screening of the colloid. System temperature and electrophoretic mobility both show a direct and inverse relationships respectively with electric field and colloidal repulsion; although they each respectively behave direct and inverse trends with salt concentration under various thermostats. Nosé–Hoover–Lowe–Andersen and Lowe–Andersen thermostats are found to function more effectively under high electric fields (E>0.145v/nm) while thermal equilibrium is maintained. Reasonable agreements are achieved by benchmarking the system radial distribution function with available EW3D modellings, as well as comparing reduced mobility against conventional Smoluchowski and Hückel theories, and numerical solution of Poisson–Boltzmann equation.

Finite Maxwell field and electric displacement Hamiltonians derived from a current dependent Lagrangian

ABSTRACTIn the common Ewald summation technique for the evaluation of electrostatic forces, the average electric field E is strictly zero. Finite uniform E can be accounted for by adding it as a new degree of freedom in an extended Lagrangian. Representing the uniform polarization P as the time integral of the internal current and E as the time derivative of a uniform vector field A, we define such an extended Lagrangian coupling A to the total current jt (internal plus external) and hence derive a Hamiltonian resembling the minimal coupling Hamiltonian of electrodynamics. Next, applying a procedure borrowed from nonrelativistic molecular electrodynamics the jt · A coupling is transformed to P · D form where D is the electric displacement acting as an electrostatic boundary condition. The resulting Hamiltonian is identical to the constant-D Hamiltonian obtained by Stengel, Spaldin and Vanderbilt (SSV) using thermodynamic arguments. The corresponding SSV constant-E Hamiltonian is derived from an alternative extended Lagrangian.

Rapid sampling of stochastic displacements in Brownian dynamics simulations with stresslet constraints.

Brownian Dynamics simulations are an important tool for modeling the dynamics of soft matter. However, accurate and rapid computations of the hydrodynamic interactions between suspended, microscopic components in a soft material are a significant computational challenge. Here, we present a new method for Brownian dynamics simulations of suspended colloidal scale particles such as colloids, polymers, surfactants, and proteins subject to a particular and important class of hydrodynamic constraints. The total computational cost of the algorithm is practically linear with the number of particles modeled and can be further optimized when the characteristic mass fractal dimension of the suspended particles is known. Specifically, we consider the so-called "stresslet" constraint for which suspended particles resist local deformation. This acts to produce a symmetric force dipole in the fluid and imparts rigidity to the particles. The presented method is an extension of the recently reported positively split formulation for Ewald summation of the Rotne-Prager-Yamakawa mobility tensor to higher order terms in the hydrodynamic scattering series accounting for force dipoles [A. M. Fiore et al., J. Chem. Phys. 146(12), 124116 (2017)]. The hydrodynamic mobility tensor, which is proportional to the covariance of particle Brownian displacements, is constructed as an Ewald sum in a novel way which guarantees that the real-space and wave-space contributions to the sum are independently symmetric and positive-definite for all possible particle configurations. This property of the Ewald sum is leveraged to rapidly sample the Brownian displacements from a superposition of statistically independent processes with the wave-space and real-space contributions as respective covariances. The cost of computing the Brownian displacements in this way is comparable to the cost of computing the deterministic displacements. The addition of a stresslet constraint to the over-damped particle equations of motion leads to a stochastic differential algebraic equation (SDAE) of index 1, which is integrated forward in time using a mid-point integration scheme that implicitly produces stochastic displacements consistent with the fluctuation-dissipation theorem for the constrained system. Calculations for hard sphere dispersions are illustrated and used to explore the performance of the algorithm. An open source, high-performance implementation on graphics processing units capable of dynamic simulations of millions of particles and integrated with the software package HOOMD-blue is used for benchmarking and made freely available in the supplementary material.

Comparison of Additive and Polarizable Models with Explicit Treatment of Long-Range Lennard-Jones Interactions Using Alkane Simulations.

Long-range Lennard-Jones (LJ) interactions have a significant impact on the structural and thermodynamic properties of nonpolar systems. While several methods have been introduced for the treatment of long-range LJ interactions in molecular dynamics (MD) simulations, increased accuracy and extended applicability is required for anisotropic systems such as lipid bilayers. The recently refined Lennard-Jones particle-mesh Ewald (LJ-PME) method extends the particle-mesh Ewald (PME) method to long-range LJ interactions and is suitable for use with anisotropic systems. Implementation of LJ-PME with the CHARMM36 (C36) additive and CHARMM Drude polarizable force fields improves agreement with experiment for density, isothermal compressibility, surface tension, viscosity, translational diffusion, and 13C T1 relaxation times of pure alkanes. Trends in the temperature dependence of the density and isothermal compressibility of hexadecane are also improved. While the C36 additive force field with LJ-PME remains a useful model for liquid alkanes, the Drude polarizable force field with LJ-PME is more accurate for nearly all quantities considered. LJ-PME is also preferable to the isotropic long-range correction for hexadecane because the molecular order extends to nearly 20 Å, well beyond the usual 10-12 Å cutoffs used in most simulations.