I am working on molecular dynamics simulation of a crystal system, starting with neutral Li2SiO3. My goal is to study the effect on Li ion mobility when substituting Si with P, creating Li2PO3, which introduces a net charge in the system. I need a way to add a uniform background charge to maintain neutrality without affecting the dynamics, particularly the mobility of the Li ions.
In LAMMPS, what would be the best approach to introduce a neutralizing background that does not interfere with the dynamics of Li ions? Are there specific commands or methods that can achieve this in molecular dynamics simulations?
First of all, having a charged simulation cell with periodic boundaries is creating a system that is by definition unstable and has an infinite energy due to the repulsion of the cells.
That said, a uniform charge background by definition has no effect on forces since its derivative in every direction is zero.
Many papers, when having a charged unitcell and using a long-range lattice solver (like Ewald summation, PME, or PPPM or variants of it) often claim that they apply a uniform charge jellium, but that is not exactly what happens. Rather they act as if the system was neutral. The lattice sum would then result in a diverging term for the k=0 reciprocal vector and that is simply ignored. This is effectively equivalent to applying a uniform background charge which is equal but of opposite sign than the diverging term.
However, this does have some subtle effects on the overall dielectric behavior and how the lattice summation influences the dielectric properties of the simulated medium.
Thus the bottom line is, that if you care this much about charge neutrality, you have to make your system cell neutral through adding particles with opposite charge. Otherwise, you will have to accept some subtle and indirect consequences that are hard to predict and will be the larger the large the system is removed from neutrality. Adding counter charges usually requires using a very large system.
Bottom line, there is no simple choice, you have to pick the evil that you can tolerate best.
I still stand by my conclusion that such simulations are often useless. In this particular case, the material Li2PO3 doesn’t exist as a stable bulk material (since it’s non-neutral!), and calculating the diffusivity of lithium in a non-existent material may not be the best use of computer time.
You would be better served either calculating for an actually stable material (so that experimental verification is possible at least in principle), or working out a more physically valid way to calculate whatever else you may be wanting to calculate.
I understand your points, but in my case, I have fixed the P and O ions at their initial positions throughout the simulation to avoid any structural instability. What I’m particularly interested in is the uniform jellium background that is mentioned in some papers. Does this background require any extra command in LAMMPS, or is it implicitly handled by the Ewald summation?
For my setup, I plan to run the simulation in an NVT ensemble, with the P and O ions fixed and only Li ions allowed to move. Can I simply proceed with this simulation without any additional commands for the background charge, and expect to obtain practical results regarding the effect of P on the Li dynamics?
Thank you for the detailed insights! I’ve reviewed the points, and they clarify a lot. I also understand the context provided about the non neutral system. However, in my specific case, I need to proceed with this non-physical setup intentionally, mainly to isolate the effect of phosphorus on lithium mobility.
Given that I’m constraining P and O atoms in place to maintain structural stability, could you guide me on implementing a jellium background effectively in this context? I’m looking for any specific commands or tweaks that might be necessary to mimic a uniform background charge without skewing the lithium dynamics too much. If simply using Ewald summation with these fixed constraints suffices, please confirm if this is a viable and realistic approach.
As Yoda would say: do, or do not; there is no “dear Internet volunteer, please confirm a mathematically correct way to simulate a physically impossible system”.
If you can find a particular prior paper that has taken a similar approach (simulating a non-physical bulk charged crystal to measure … something), and you are not sure how to duplicate a particular step in their methodology, then we can take a look and determine what LAMMPS settings would replicate their results.
If you can’t find a prior paper, then you’re trying to do something that nobody else has done before, and you will simply have to wait to obtain the results and verify them somehow. I’m not saying this to dissuade you: methods development is extremely rewarding for any academic, especially early-career, since anyone using a method you’ve pioneered and benchmarked will have to cite your paper no matter what new thing they’re trying to do. But you need a really firm grip on theory to develop methods, and you’ll have to know how to verify theoretical soundness and experimental usefulness when you’re literally the first person in the world to have tried something.