Hi everyone,
Could someone explain why we typically neutralize the system by adding ions when using periodic boundary conditions and Particle Mesh Ewald (PME) in molecular dynamics simulations? What potential issues might arise if we don’t neutralize the system?
Thanks!
First off, you posted to the wrong category. This is a question for the “Science Talk” category. I’ve moved your post there.
If your system is not neutral, the lattice sum in the Ewald summation part will diverge. The energy of a crystal would be infinite if a unit cell is charged. It would just explode.
Now what is done in about every Ewald sum (and similar algorithms) implementation is that the (infinite) residual term at k=0 is ignored because it - technically speaking - will not result it a force.
So, you neutralize a system to avoid simulating an unphysical system.
As with most (serious) questions in science, the answer is: it depends! Technically speaking, the balance between the real space and reciprocal space force calculation is no longer given for a charged unit cell. If the residual charge is small, the error is small and may be smaller than the systematic error for the choice of force field or from using a classical model (instead of quantum chemistry). Or smaller than the error from having less trajectory when simulating a neutral system (which may require simulating a much larger system) for a fixed amount of available computing resources.
Bottom line, in most cases the simulation results are significantly tainted, if the system is not neutral with a lattice sum based long-range Coulomb solver.
Thanks for your response but what if we want to simulate a charged system? Does this impact our results for simulating a charged system?
There is always a counter charge somewhere…
A charged system in real live can only exist when you separate charges.
The problem is really a problem of using a lattice sum to improve the accuracy of the Coulumb forces. If you use a cutoff coulomb or multi-level summation, then you are not doing calculations assuming an infinite lattice and performing part of your calculation in reciprocal space.
I have already answered that question.
You can find a recent, straightforward introduction to proper treatment of periodic charged systems in Ballenegger, Cerda and Arnold’s 2009 paper on charged slabs (Simulations of non-neutral slab systems with long-range electrostatic interactions in two-dimensional periodic boundary conditions | The Journal of Chemical Physics | AIP Publishing ; this paper is cited in the LAMMPS kspace modify documentation). Depending on your level of analysis, the results of such simulations are either mostly right, mostly wrong, or mostly useless.
Charged periodic simulations can be mostly right because you can simply imagine smearing a counter charge density uniformly throughout space. You don’t even need extra calculations because such a space charge exerts no forces. (As renowned physicist Syndrome might say, “when all space is charged … no space is charged.”) It simply amounts to a static added self-energy term for the counter-charge.
But such simulations can be mostly wrong if there is a real counter charge with physically significant behaviour. A simple, intuitive example is lipid bilayer membrane simulations. Suppose you simulate a lipid bilayer with explicit water and ions, but the ions have a net charge. The uniform counter-charge described earlier exists even within the lipid tail volume. This is non-physical – the dielectric constant is much lower in the lipid tail volume, and in a real system all charged and polar molecules will rarely enter that space. So, relative to a realistic simulation, in a charged simulation water molecules and ions will penetrate far too often into the membrane space.
Which is why charged simulations are often useless – if your system volume is so small (relative to the macroscopic system of interest) that you cannot afford to put in the counter-charge, it is almost certainly so small that you will have other significant finite size effects that destroy any realism in the simulation results. My default judgement would be to doubt such a simulation until proven otherwise.
To give another example of how charged systems are tricky, check out this paper (which is one of my all-time favorite) from Wong-Ekkabut and Karttunen.
The case studies are other examples of how quickly wrong simulation setup, especially electrostatic forces computation, can lead you to wrong results and conclusions, even with “simple” systems.
Thanks – this is indeed an excellent resource. Having said that, the particular paper is slightly dated – the choice of reaction field electrostatics is far less common in 2024 than in 2016 and prior.
More importantly, many of the examples arise from the peculiarities of a software that is, shall we say, Generally Respectable – Occasionally Making Absolutely Catastrophic Splats. The piling up of several individually acceptable shortcuts in neighbourlisting, in particular, makes for some eye-popping parameter sensitivity – see this paper: Neighbor List Artifacts in Molecular Dynamics Simulations - PMC from 2023!!
LAMMPS’s avoidance of charge groups (except, technically, somewhat in the TIP4P styles) and more conservative neighbourlisting choices do avoid some of the worse algorithmic errors. The basic point of the paper is useful and well-noted though: user, know your results! This graphic in particular is lovely:
