Simplistic Coulomb Forces in Molecular Dynamics: Comparing the Wolf and Shifted-Force Approximations

Modified and generalized single-element Maxwell viscoelastic model

In this Letter, the single-element Maxwell model is generalized with respect to the wave vector and extended with a correction function that measures the reduced viscous response. This model has only two free parameters and avoids the attenuation-frequency locking present in the original model. Through molecular simulations it is shown that the model satisfactory predicts the transverse dynamics of the binary Lennard-Jones system at different temperatures, as well as water and toluene at ambient conditions. The correction function shows that the viscous response is significantly reduced compared to the predictions of the original Maxwell model and that there exists a characteristic length scale of minimum dissipation.

Unification of Ewald and shifted force methods to calculate Coulomb interactions in molecular simulations

Three new Ewald series are derived using a new strategy that does not start with a proposed charge spreading function. Of these, the Ewald series produced using shifted potential interactions for the Ewald real space series converges relatively slowly, while the corresponding expression using a shifted force (SF) interaction does not converge. A comparison is made between several approximations of the Ewald method and the SF route to include Coulomb interactions in molecular dynamics (MD) computer simulations. MD simulations of a model bulk molten salt and water were carried out. The recently derived α' variant of Ewald, by K. D. Hammonds and D. M. Heyes [J. Chem. Phys. 157, 074108 (2022)], has been developed analytically and found to be more accurate and computationally efficient than SF in part due to the smaller real space truncation distance that can be used. In addition, with α', the number of reciprocal lattice vectors required is reduced considerably compared with the standard Ewald implementations to give the same accuracy. The invention of the α' method shifts the computational balance back toward using an Ewald construction. The SF method shows greater errors in the Coulomb pressure and time dependent fluctuation properties compared to α'. It does not conserve the shadow Hamiltonian in a microcanonical MD simulation, whereas the α' method does, which facilitates long time stability and insignificant drift of properties over time. The speed of the Ewald computer code is improved by using a new lookup table method.

Quantifying dynamical and structural invariance in a simple molten salt model

Recent experimental results for the structure in the ionic liquid PYR₁₄ ⁺TFSI^- have shown invariance in the main structure factor peak along curves of equal electrical conductivity [Hansen et al., Phys. Chem. Chem. Phys. 22, 14169 (2020)]. The charge peak decreases slightly with increasing temperature at fixed conductivity, however. For simple liquids, curves with invariant dynamics and structure, known as isomorphs, can be identified as configurational adiabats. While liquids with strong-Coulomb interactions do not have good isomorphs, ionic liquids could be an intermediate case with approximate isomorphs along which some aspects of structure and dynamics are invariant. We study a simple molten salt model using molecular dynamics simulations to test this hypothesis. Simple measures of structure and dynamics are investigated along with one transport property, the shear viscosity. We find that there is a substantial degree of invariance of the self-intermediate scattering function, the mean square displacement, and the viscosity along configurational adiabats over a wide range of densities for the three adiabats simulated. The density range studied is more than a factor of two and extends from the strong-Coulomb regime at low densities to the weak-Coulomb regime at high densities. The structure is not invariant over the full range of density, but in the weak-Coulomb regime, we see behavior similar to that seen experimentally over density changes of order 15%. In view of the limited structural invariance but substantial dynamical invariance, we designate the configurational adiabats as isodynes.

Shadow Hamiltonian in classical NVE molecular dynamics simulations involving Coulomb interactions

Microcanonical ensemble (NVE) Molecular Dynamics (MD) computer simulations are performed with negligible energy drift for systems incorporating Coulomb interactions and complex constraint schemes. In principle, such systems can now be simulated in the NVE ensemble for millisecond time scales, with no requirement for system thermostatting. Numerical tools for assessing drift in MD simulations are outlined, and drift rates of 10^-6 K/μs are demonstrated for molten salts, polar liquids, and room temperature ionic liquids. Such drift rates are six orders of magnitude smaller than those typically quoted in the literature. To achieve this, the standard Ewald method is slightly modified so the first four derivatives of the real space terms go smoothly to zero at the truncation distance, r_c. New methods for determining standard Ewald errors and the new perturbation errors introduced by the smoothing procedure are developed and applied, these taking charge correlation effects explicitly into account. The shadow Hamiltonian, E_s, is shown to be the strictly conserved quantity in these systems, and standard errors in the mean of one part in 10¹⁰ are routinely calculated. Expressions for the shadow Hamiltonian are improved over previous work by accounting for O(h⁴) terms, where h is the MD time step. These improvements are demonstrated by means of extreme out-of-equilibrium simulations. Using the new methodology, the very low diffusion coefficients of room temperature 1-hexyl-3-methyl-imidazolium chloride are determined from long NVE trajectories in which the equations of motion are known to be integrated correctly, with negligible drift.

Dielectric continuum model examination of real-space electrostatic treatments

Electrostatic interaction is long ranged; thus, the accurate calculation is not an easy task in molecular dynamics or Monte Carlo simulations. Though the rigorous Ewald method based on the reciprocal space has been established, real-space treatments have recently become an attractive alternative because of the efficient calculation. However, the construction is not yet completed and is now a challenging subject. In an earlier theoretical study, Neumann and Steinhauser employed the Onsager dielectric continuum model to explain how simple real-space cutoff produces artificial dipolar orientation. In the present study, we employ this continuum model to explore the fundamental properties of the recently developed real-space treatments of three shifting schemes. The result of the distance-dependent Kirkwood function G_K(R) showed that the simple bare cutoff produces a well-known hole-shaped artifact, whereas the shift treatments do not. Two-dimensional mapping of electric field well explained how these shift treatments remove the hole-shaped artifact. Still, the shift treatments are not sufficient because they do not produce a flat G_K(R) profile unlike ideal no-cutoff treatment. To test the continuum model results, we also performed Monte Carlo simulations of dipolar particles. The results found that the continuum model could predict the qualitative tendency as to whether each electrostatic treatment produces the hole-shaped artifact of G_K(R) or not. We expect that the present study using the continuum model offers a stringent criterion to judge whether the primitive electrostatic behavior is correctly described or not, which will be useful for future construction of electrostatic treatments.

Assessment of the Wolf method using the Stillinger-Lovett sum rules: From strong electrolytes to weakly charged colloidal dispersions

The Ewald method has been the cornerstone in molecular simulations for modeling electrostatic interactions of charge-stabilized many-body systems. In the late 1990s, Wolf and collaborators developed an alternative route to describe the long-range nature of electrostatic interactions; from a computational perspective, this method provides a more efficient and straightforward way to implement long-range electrostatic interactions than the Ewald method. Despite these advantages, the validity of the Wolf potential to account for the electrostatic contribution in charged fluids remains controversial. To alleviate this situation, in this contribution, we implement the Wolf summation method to both electrolyte solutions and charged colloids with moderate size and charge asymmetries in order to assess the accuracy and validity of the method. To this end, we verify that the proper selection of parameters within the Wolf method leads to results that are in good agreement with those obtained through the standard Ewald method and the theory of integral equations of simple liquids within the so-called hypernetted chain approximation. Furthermore, we show that the results obtained with the original Wolf method do satisfy the moment conditions described by the Stillinger-Lovett sum rules, which are directly related to the local electroneutrality condition and the electrostatic screening in the Debye-Hückel regime. Hence, the fact that the solution provided by the Wolf method satisfies the first and second moments of Stillinger-Lovett proves, for the first time, the reliability of the method to correctly incorporate the electrostatic contribution in charge-stabilized fluids. This makes the Wolf method a powerful alternative compared to more demanding computational approaches.

Electrospray molecular dynamics simulations using an octree-based Coulomb interaction method

A new octree-based Coulomb interaction model is developed to model the electrospray of ionic liquids (ILs) in molecular dynamics. Using an octree-based method, Coulomb interactions are categorized as intra- and interleaf Coulomb interactions based on a criterion related to the Bjerrum length of the IL. The octree-based method is found capable of reproducing Coulomb energy in agreement with established and computationally more expensive models, such as the direct Coulomb and the damped shifted force (DSF) method in the absence of an external electric field. In the presence of an external electric field, the octree-based method produces distinctly different results compared to that obtained by the direct Coulomb method. The time required to form Taylor's cone was shorter for the octree method compared to the direct Coulomb approach. While no emission larger than monomers was observed from the direct Coulomb simulation, emission of larger species such as dimers and trimers was observed when the octree-based Coulomb interaction model was used. Furthermore, the octree-based model forms a smaller ion emission cone compared to that from the direct Coulomb method.

A Modified Shifted Force Approach to the Wolf Summation

The Wolf method for calculation of electrostatic interactions in molecular simulations is known to describe the energy well, whereas the forces have discontinuities. For a more reliable description of the forces this method can be extended with a shifted force approach. This leads to a good description of the forces and precise molecular dynamics simulation, but the description of the energy becomes poorer. In this study we propose a modification of a shifted force extension to describe the energy as well as the forces in better agreement to reference data as determined from the Ewald summation. We show that vapor-liquid phase equilibria (VLE) calculated with Monte Carlo simulations in the grand canonical ensemble and dynamic properties calculated with molecular dynamics simulations can be calculated reliably using this modification to describe the electrostatic interactions.

Single-pass incremental force updates for adaptively restrained molecular dynamics

Adaptively restrained molecular dynamics (ARMD) allows users to perform more integration steps in wall-clock time by switching on and off positional degrees of freedoms. This article presents new, single-pass incremental force updates algorithms to efficiently simulate a system using ARMD. We assessed different algorithms for speedup measurements and implemented them in the LAMMPS MD package. We validated the single-pass incremental force update algorithm on four different benchmarks using diverse pair potentials. The proposed algorithm allows us to perform simulation of a system faster than traditional MD in both NVE and NVT ensembles. Moreover, ARMD using the new single-pass algorithm speeds up the convergence of observables in wall-clock time. © 2017 Wiley Periodicals, Inc.

Electric fields can control the transport of water in carbon nanotubes

The properties of water confined inside nanotubes are of considerable scientific and technological interest. We use molecular dynamics to investigate the structure and average orientation of water flowing within a carbon nanotube. We find that water exhibits biaxial paranematic liquid crystal ordering both within the nanotube and close to its ends. This preferred molecular ordering is enhanced when an axial electric field is applied, affecting the water flow rate through the nanotube. A spatially patterned electric field can minimize nanotube entrance effects and significantly increase the flow rate.

Hybrid Cassie-Wenzel model for droplets on surfaces with nanoscale roughness

Several models have been developed to predict the contact angle of a droplet sitting on a roughened surface; however, no such model has been developed for substrates with nanoscale surface structures. In this paper we propose a hybrid Cassie-Wenzel model, which considers two factors attributed to the breakdown of macroscopic predictions, including the width of the wall-fluid depletion region and the coexistence of Cassie and Wenzel states in cases where the wall-fluid interface presents nanoscale structures. At the molecular scale, the parameter of surface roughness can be corrected by treating the wall-fluid interface as a hybrid Cassie-Wenzel state in which the fraction in the Wenzel state depends on fluid density within the cavities. A more general model developed using data fitted to fluid density is able to account for deviating tendencies induced by nanoscale surface features. A comparison of predicted results obtained in this study with those from previous works demonstrates that the proposed hybrid Cassie-Wenzel model is applicable to the evaluation of wettability in a wide range of substrates with nanoscale surface structures, corresponding to a Cassie state, a Wenzel state, and a mixed state. More importantly, the present work provides a quantitative approach to the estimation of wettability even amidst nanoscale effects, which can have a significant influence in cases with surface features at the molecular scale.

Using molecular dynamics to quantify the electrical double layer and examine the potential for its direct observation in the in-situ TEM

Understanding the fundamental processes taking place at the electrode-electrolyte interface in batteries will play a key role in the development of next generation energy storage technologies. One of the most fundamental aspects of the electrode-electrolyte interface is the electrical double layer (EDL). Given the recent development of high spatial resolution in-situ electrochemical fluid cells for scanning transmission electron microscopy (STEM), there now exists the possibility that we can directly observe the formation and dynamics of the EDL. In this paper we predict electrolyte structure within the EDL using classical models and atomistic Molecular Dynamics (MD) simulations. Classical models are found to greatly differ from MD in predicted concentration profiles. It is thus suggested that MD must be used in order to accurately predict STEM images of the electrode-electrolyte interface. Using MD and image simulation together for a high contrast electrolyte (the high atomic number CsCl electrolyte), it is determined that, for a smooth interface, concentration profiles within the EDL should be visible experimentally. When normal experimental parameters such as rough interfaces and low-Z electrolytes (like those used in Li-ion batteries) are considered, observation of the EDL appears to be more difficult.

Accelerating electrostatic pair methods on graphical processing units to study molecules in supercritical carbon dioxide

Traditionally, electrostatic interactions are modelled using Ewald techniques, which provide a good approximation, but are poorly suited to GPU architectures. We use the GPU versions of the LAMMPS MD package to implement and assess the Wolf summation method. We compute transport and structural properties of pure carbon dioxide and mixtures of carbon dioxide with either methane or difluoromethane. The diffusion of pure carbon dioxide is indistinguishable when using the Wolf summation method instead of PPPM on GPUs. The optimum value of the potential damping parameter, α, is 0.075. We observe a decrease in accuracy when the system polarity increases, yet the method is robust for mildly polar systems. We anticipate the method can be used for a number of techniques, and applied to a variety of systems. Substitution of PPPM can yield a two-fold decrease in the wall-clock time.

Molecular dynamics study of nanoconfined water flow driven by rotating electric fields under realistic experimental conditions

In our recent work, J. Chem. Phys. 2013, 138, 154712, we demonstrated the feasibility of unidirectional pumping of water, exploiting translational-rotational momentum coupling using nonequilibrium molecular dynamics simulations. Flow can be sustained when the fluid is driven out of equilibrium by an external spatially uniform rotating electric field and confined between two planar surfaces exposing different degrees of hydrophobicity. The permanent dipole moment of water follows the rotating field, thus inducing the molecules to spin, and the torque exerted by the field is continuously injected into the fluid, enabling a steady conversion of spin angular momentum into linear momentum. The translational-rotational coupling is a sensitive function of the rotating electric field parameters. In this work, we have found that there exists a small energy dissipation region attainable when the frequency of the rotating electric field matches the inverse of the dielectric relaxation time of water and when its amplitude lies in a range just before dielectric saturation effects take place. In this region, that is, when the frequency lies in a small window of the microwave region around ∼20 GHz and amplitude ∼0.03 V Å(-1), the translational-rotational coupling is most effective, yielding fluid velocities of magnitudes of ∼2 ms(-1) with only moderate fluid heating. In this work, we also confine water to a realistic nanochannel made of graphene giving a hydrophobic surface on one side and β-cristobalite giving a hydrophilic surface on the other, reproducing slip-and-stick velocity boundary conditions, respectively. This enables us to demonstrate that in a realistic environment, the coupling can be effectively exploited to achieve noncontact pumping of water at the nanoscale. A quantitative comparison between nonequilibrium molecular dynamics and analytical solutions of the extended Navier-Stokes equations, including an external rotating electric field has been performed, showing excellent agreement when the electric field parameters match the aforementioned small energy dissipation region.

How fast does water flow in carbon nanotubes?

The purpose of this paper is threefold. First, we review the existing literature on flow rates of water in carbon nanotubes. Data for the slip length which characterizes the flow rate are scattered over 5 orders of magnitude for nanotubes of diameter 0.81-10 nm. Second, we precisely compute the slip length using equilibrium molecular dynamics (EMD) simulations, from which the interfacial friction between water and carbon nanotubes can be found, and also via external field driven non-equilibrium molecular dynamics simulations (NEMD). We discuss some of the issues in simulation studies which may be reasons for the large disagreements reported. By using the EMD method friction coefficient to determine the slip length, we overcome the limitations of NEMD simulations. In NEMD simulations, for each tube we apply a range of external fields to check the linear response of the fluid to the field and reliably extrapolate the results for the slip length to values of the field corresponding to experimentally accessible pressure gradients. Finally, we comment on several issues concerning water flow rates in carbon nanotubes which may lead to some future research directions in this area.

Dynamics of nanoscale droplets on moving surfaces

We use molecular dynamics (MD) simulations to investigate the dynamic wetting of nanoscale water droplets on moving surfaces. The density and hydrogen bonding profiles along the direction normal to the surface are reported, and the width of the water depletion layer is evaluated first for droplets on three different static surfaces: silicon, graphite, and a fictitious superhydrophobic surface. The advancing and receding contact angles, and contact angle hysteresis, are then measured as a function of capillary number on smooth moving silicon and graphite surfaces. Our results for the silicon surface show that molecular displacements at the contact line are influenced greatly by interactions with the solid surface and partly by viscous dissipation effects induced through the movement of the surface. For the graphite surface, however, both the advancing and receding contact angles values are close to the static contact angle value and are independent of the capillary number; i.e., viscous dissipation effects are negligible. This finding is in contrast with the wetting dynamics of macroscale water droplets, which show significant dependence on the capillary number.

NVU dynamics. III. Simulating molecules at constant potential energy

This is the final paper in a series that introduces geodesic molecular dynamics at constant potential energy. This dynamics is entitled NVU dynamics in analogy to standard energy-conserving Newtonian NVE dynamics. In the first two papers [T. S. Ingebrigtsen, S. Toxvaerd, O. J. Heilmann, T. B. Schrøder, and J. C. Dyre, J. Chem. Phys. 135, 104101 (2011); T. S. Ingebrigtsen, S. Toxvaerd, T. B. Schrøder, and J. C. Dyre, ibid. 135, 104102 (2011)], a numerical algorithm for simulating geodesic motion of atomic systems was developed and tested against standard algorithms. The conclusion was that the NVU algorithm has the same desirable properties as the Verlet algorithm for Newtonian NVE dynamics, i.e., it is time-reversible and symplectic. Additionally, it was concluded that NVU dynamics becomes equivalent to NVE dynamics in the thermodynamic limit. In this paper, the NVU algorithm for atomic systems is extended to be able to simulate the geodesic motion of molecules at constant potential energy. We derive an algorithm for simulating rigid bonds and test this algorithm on three different systems: an asymmetric dumbbell model, Lewis-Wahnström o-terphenyl (OTP) and rigid SPC/E water. The rigid bonds introduce additional constraints beyond that of constant potential energy for atomic systems. The rigid-bond NVU algorithm conserves potential energy, bond lengths, and step length for indefinitely long runs. The quantities probed in simulations give results identical to those of Nosé-Hoover NVT dynamics. Since Nosé-Hoover NVT dynamics is known to give results equivalent to those of NVE dynamics, the latter results show that NVU dynamics becomes equivalent to NVE dynamics in the thermodynamic limit also for molecular systems.