Communication: Truncated non-bonded potentials can yield unphysical behavior in molecular dynamics simulations of interfaces

Short- to medium-range potentials such as the Lennard-Jones¹ or the Buckingham² potential are the backbone of classical molecular dynamics (MD) simulations. They represent Pauli repulsion as well as non-directional dispersion attraction and there exist multiple flavors implemented in most MD codes under the term of non-bonded interactions. In practice, there is a need to truncate these potentials since the number of neighbors that have to be considered for each entity grows enormously, drastically increasing the computational cost for the force calculation. Truncating between r_c = 2.5σ and 3.5σ, where σ is the characteristic interaction range, is a very common practice in MD studies³ and has become the minimum standard, assuming that errors arising from this are small enough. Several studies have reported that with these settings significant problems can arise. For instance, the truncation can alter the phase diagram of the Lennard-Jones system^4,5 or yield different values for interfacial free energies.^6–10 These effects are quantitative in nature, meaning that they can in certain circumstances be analytically corrected for^11–13 or compensated for by other interaction parameters such as interaction strength or interaction range. The latter is important for the development of force fields where non-bonded potentials are often included, and the cutoff can be seen as another fitting parameter. Naturally, a parametrization with a small cutoff would be preferred to another one if they deliver equal accuracy. This however is only true in the assumption that the underlying physical characteristics that are created by truncated and longer ranging potentials are the same.

In this work, we investigated the influence of the cutoff for the interfacial phenomenon of water-wetting on a solid substrate. We found that the effect of the cutoff of the water-substrate interaction was not only unexpectedly strong but also changed the fundamental physics of the wetting transition in an unprecedented way by creating metastable wetting states that have also never been seen in experiments. We show that proposed cutoff corrections such as the use of a force switching function, a shifted potential, or a shifted force did not fix this and could even worsen the effect. This finding shows that atomistic simulations of interfaces need to be treated with great care since unphysical behavior could occur and easily remain undetected. This is particularly relevant since a large number of MD studies using truncated potentials are reported each year. Our results suggest the use of much larger-than-common cutoffs or long-range versions of non-bonded potentials in MD studies of wetting and interfacial phenomena.

We investigated two droplets composed of 3000 and 18 000 water molecules which were represented by the coarse-grained mW model,¹⁴ on top of a rigid, pristine fcc(100) surface (lattice parameter 4.15 Å). Whilst this substrate does not aim at representing any particular material, similar systems have been used to study ice nucleation^15–18 or water-metal interfaces.^19,20 The simulation cell had dimensions $17×17×11$ nm, which is enough to avoid interaction of the water molecules with their periodic images for all wetting states. Even though the liquid is rather non-volatile even at the highest temperature considered, we employed a reflective wall at the top of the cell to avoid evaporation and mimic experimental conditions. Our simulations were performed with the LAMMPS code,²¹ integrating the equations of motion with a time step of 10 fs. This rather large time step is commonly used in combination with the mW model and is acceptable for our system since during NVE simulations the total energy drift was found to be only about $2×1 0 − 9$ eV per water molecule per ps. In addition, we verified that we obtain the same results using standard protocols for updating the neighbor lists compared with unconditionally updating them every time step. All production simulations were performed in the NVT ensemble with constant temperature maintained by a ten-fold Nosé-Hoover chain²² with a relaxation time of 1 ps. The substrate-water interaction was given by a distance (r) dependent Lennard-Jones potential

with $𝜖$ = 29.5 meV and σ = 2.5 Å truncated at a cutoff r_c. This resulted in a maximum interaction energy of 154 meV for an adsorbed water monomer (weakly depending on the cutoff). Additionally we performed well-tempered metadynamics simulations^23,24 for the smaller droplet with the PLUMED2 code.²⁵ In these simulations, the Gaussian height, width, bias-factor, and deposition stride were 2.16 meV, 0.15 Å, 20, and 20 ps, respectively. Metadynamics is usually applied to drive rare events such as nucleation^26–29 or protein folding.^30,31 In our systems, this method helped to uncover the underlying free energy profile of wetting.

We studied the wetting behavior of the larger droplet by performing standard MD runs at different temperatures first. As starting configurations, we chose either a flat water film in direct contact or a spherical droplet placed above the substrate. Within at most 5 ns, the simulation was equilibrated and a seemingly stable configuration was reached, where the water is either wetting (contact angle $θ= 0 °$⁠) or partially wetting (⁠$0 ° <θ<18 0 °$⁠). An illustration of the two wetting states can be found in Fig. 1(a). Initially we employed a radial cutoff at r_c = 3.0σ for the water-substrate interaction. With this setting, we found that interestingly a wetting transition happened at a finite angle $θ 0 ≈2 3 °$⁠, i.e., a smaller non-zero contact angle was not possible. This behavior cannot be explained by the standard Young’s equation.

However, upon increasing the cutoff, we found that the wetting behavior drastically changed. First, the wetting temperature T_w at which the wetting transition took place increased as we increased the cutoff [Fig. 1(b)]. Whilst T_w shows a clear convergence behavior with r_c, it is unexpectedly slow. A reasonably converged wetting temperature T₀ is only reached for $r c >7σ$⁠. Second, we noticed that for an increasing cutoff, the minimum possible contact angle $θ$₀ got smaller and eventually vanished. Most importantly, we also found that for temperatures around T_w, the stable configuration that was reached after 5 ns could depend on the starting configuration for smaller cutoffs, while for larger r_c it always reached the same state. This suggests that for small r_c, we actually found metastable wetting states that are absent for large r_c. This also means that T_w cannot naively be defined through the visual analysis of trajectories at different temperatures but needs to be defined by the free energy of wetting. For a first order phase transition, we define T_w to be the temperature where the two basins (corresponding to wetting and partial wetting) have the same free energy. For a continuous phase transition, T_w is the temperature where the single basin represents a contact angle of $θ= 0 °$ for $T< T w$ and $θ> 0 °$ for $T> T w$⁠.

Understanding the character of these wetting states with standard MD can prove difficult as the dependence on the starting configuration always leaves doubt on the outcome of the equilibrated configuration obtained from it. To clarify, we show the results from the metadynamics simulations in Fig. 2. As a collective variable, we chose the z-component of the center of mass of the water droplet (COM_z), where z is the surface normal direction. While this choice is not equivalent to the contact angle (as they are related in a non-linear manner), it is clear that significantly different values for COM_z correspond to different contact angles and can therefore distinguish the different wetting states. For the smallest cutoff at T_w and around, we found that two basins coexist, one being the flat film (⁠$COM z ≈$ 4 Å) and the other being a droplet with a certain contact angle (⁠$COM z ≳5$ Å). These two states are separated by a significant barrier larger than 20 k_BT, which explains why we observed metastable states in the unbiased simulations for small r_c. This corresponds to a first-order phase transition between the wetting states. The occurrence of a minimum possible contact angle $θ$₀ is explained by the existence of the second basin, which does not approach the wetting basin, but rather becomes less stable as the temperature changes. However, this character faded as we increased r_c. The barrier became smaller and the distance between the basins got smaller. For the largest cutoff investigated (8σ), we clearly see that only a single basin exists that changes its position with temperature. As a result, no metastable wetting states exist and the phase transition is continuous. We note that in this case, the estimate of T_w is more difficult than for the first order transitions; however, in this work, we aim at presenting qualitative results and from Fig. 2 it is clear that T_w is higher than for the smaller cutoffs.

Only the results for the largest cutoff are in agreement with the fact that water wetting transitions are generally continuous when probed in experiments^32,33 and finite-angle wetting transitions have, to the best of our knowledge, never been observed experimentally. Therefore, the correct qualitative wetting behavior in our system is not achieved with standard cutoffs and if undetected could potentially lead to false conclusions. Differences between short- and long-ranged interactions have been highlighted for other interfacial phenomena, such as drying³⁴ or grain boundary melting.³⁵ 

We further study the effect of the most commonly used correction schemes to cutoffs as follows:

1. A shifted potential (sp) which ensures that the value of the potential energy U does not jump at the cutoff distance, given by

   $U s p ( r ) = U L J ( r ) − U L J ( r c ) .$

   (2)

   The corresponding force F remains unaltered,

   $F s p ( r ) = F L J ( r ) .$

   (3)
2. A switching function (switch) which brings the force to zero between an inner cutoff r_c,1 and an outer cutoff r_c,2 (we chose 3σ and 4σ),

   $F s w i t c h ( r ) = F L J ( r ) r ≤ r c , 1 , F s w i t c h ( r ) = ∑ k = 0 3 C k ( r − r c , 1 ) k , r c , 1 < r ≤ r c , 2 ,$

   (4)

   where C_k are constants determined to ensure a smooth behavior.²¹ 
3. A shifted-force potential (sf), which ensures that force and potential do not jump,

   $U s f ( r ) = U L J ( r ) − U L J ( r c ) − ( r − r c ) F L J ( r c ) , F s f ( r ) = F L J ( r ) − F L J ( r c ) .$

   (5)

The latter approach was found to give good results for a homogeneous system and even allowed for a reduction of the cutoff.³⁶ Our results for these three corrections can be found in Fig. 3. By definition and thus unsurprisingly, the shifted potential does not yield any significant difference (where the remaining minor deviations are due to the metadynamics sampling) over the plain cutoff since the forces remain unaltered. The smooth cutoff via the switching function seems to improve the situation; however, the fact that the transition temperature lies between the ones we found for a plain cutoff at 3σ and 4σ suggests that the improvement stems from the effectively increased interaction range rather than the fact that the force vanishes smoothly. Interestingly, the shifted force with the same cutoff performs worst out of all the candidates as the barrier increases by a factor of two, which increases the likelihood that simulations are performed in the metastable state without realizing it. The fact that none of the considered correction schemes significantly improved the character of the wetting free energy profile leads us to conclude that it is not the way in which the cutting is done that matters most, but rather the effective cutoff distance as well as the overall interaction strength at that distance.

As an initial attempt to understand the results obtained, we looked at the potential energies of the various systems with the different cutoffs considered. This, however, did not reveal any obvious explanation. One possible interpretation for the creation of metastable states in our systems with a shorter cutoff can be obtained by considering the droplet state (not assuming anything about the stability relative to the film state). For a transition towards the film state, there needs to be thermal fluctuations of water molecules that are above the contact layer in the downward direction (the fact that COM_z has proven a good reaction coordinate supports this statement). With an infinite interaction range, all molecules that are loosing height contribute to these fluctuations since they have an interaction with the substrate. Therefore we expect the interaction energy to change monotonically and the free energy to follow monotonically either up or down depending on the balance of the interfacial free energies (see Fig. 2, r_c = 8σ). But if the interaction range is finite, not all molecules contribute to an increased interaction with the substrate even if they decrease their height (and subsequently weaken the water-water interaction of the system by leading to deviations from a perfect spherical droplet). In other words, there is a minimum distance from the substrate that has to be surpassed by a molecule for it to contribute to a fluctuation increasing the interaction energy, otherwise it will (on average) actually decrease the total interaction energy. This minimum fluctuation for a single molecule translates into the macroscopic states (droplet and film) being connected by a barrier shaped free energy profile rather than a monotonic one (see Fig. 2, r_c = 3σ). The entropic contributions to the free energy are unlikely to change this since they are essentially dominated by the environment a molecule is in (quasi-static contact layer or quasi-liquid water on top). The entropic change between these two states will be monotonic for a single water molecule and therefore also for the whole droplet.

Finding a general recipe for how to avoid such unphysical wetting states is difficult. Other aspects like, e.g., the substrate density or the liquid-liquid interaction strength will have an influence on how strongly the fluctuations in the droplet state are affected by r_c. Generally, cutoffs that are deemed acceptable from the inter-molecular perspective do not necessarily mean that the interaction between macroscopic states such as a film/droplet and a substrate is sufficiently captured. This is especially important in an interfacial simulation setting such as a slab, where a cutoff-caused change in the interaction from the substrate side is not compensated by an equal change from the vacuum side. Consequently, only employing much larger cutoffs or techniques to calculate the long-range part of the dispersion force^37–39 can ensure that unphysical effects are avoided. A minimal sanity check for future wetting studies could be to start simulations from both a wetting film and a spherical liquid snapshot. If both of them end up in the same configuration, the existence of an unphysical metastable wetting state is unlikely.

In light of the vast amount of work that is done in the MD community using similar interactions, our findings urge extreme caution when dealing with truncated non-bonded potentials in simulations of interfacial phenomena. We have seen both quantitative and qualitative differences for the wetting transition. The former could be accounted for by changing other interaction parameters to reproduce the transition at the right temperature T₀. This assumption is fundamental to fitting force fields with truncated potentials to obtain quantitative agreement with, e.g., experimental values. But it does not hold for the character of the transition because it arises purely from the value of the cutoff itself. If the resulting metastability of states remains undetected, the use of truncated interaction potentials could lead to wrong inferences about physical properties being made. While this conclusion has resulted from a simulation of wetting, similar implications could hold for other interfacial phenomena such as capillary flow,^40,41 evaporation/condensation,^42,43 mixtures,^44–46 or heterogeneous nucleation^47–51 where it is commonplace to use truncated interactions.

This work was supported by the European Research Council under the European Union’s Seventh Framework Programme (No. FP/2007-2013)/ERC Grant Agreement No. 616121 (HeteroIce project). A.M. is supported by the Royal Society through a Royal Society Wolfson Research Merit Award. We are grateful for computational resources provided by the London Centre for Nanotechnology and the Materials Chemistry Consortium through the EPSRC Grant No. EP/L000202. L.J. is supported by the French Ministry of Defense through the Project DGA ERE No. 2013.60.0013 and by the LABEX iMUST (ANR-10-LABX-0064) of Université de Lyon, within the program “Investissements d’Avenir” (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). M.M. is supported by the Thousand Young Talent Program from the Organization Department of the CPC Central Committee.