There are several ways to include long-range interactions when modeling granular spheres in LAMMPS. For instance, one can use atom_style hybrid charge sphere, where the center of the sphere coincides with the charge center. However, as far as I know, there is no built-in option to define an offset between the charge and the sphere center.
My question is: is there any strategy to work around this using only LAMMPS commands? Or would it be necessary to implement a new pair_style?
For example, you could define two points, one with only the granular interactions and charge 0 and the second with only the charge, but no significant diameter. If you give each of these pairs a different molecule ID, you can use fix rigid/small for time integration and that should implement the model you describe.
This kind of system can be conveniently created using the molecule and the create_atoms command, unless you have specific needs and then you need to create a suitable data file with an external tool.
Additionally to @akohlmey’s suggestion of using fix rigid/small – if your systems are sufficiently simple you may also use fix shake to hold the charge centres at a specified distance from the sphere centres, by adding a bond between the two centres.
The main requirement of this method is that every bond of the same “type” has to be of equal length across the system – whereas fix rigid takes an arbitrary arrangement of initial bond lengths (and angles, dihedrals, etc).
I don’t think this can work, since this will require the point charge particle to have a significant mass and it also can freely rotate around the center of mass of the extended particle.
To me, the situation is somewhat comparable to implementing TIP4P. You either need to use a rigid fix to make the entire group of particles rigid and move as one unit or a dedicated pair style where you can handle the impact of the point charge on the rotation of the extended particle and you would derive from the orientation. And that may not possible unless you use the asphere atom style instead of sphere. For long-range Coulomb, also the corresponding long-rage solvers would be needed to complement the custom pair style.
Funny you should mention that – I have been working on some explicit-charge TIP4P simulations recently.
I faced significant problems equilibrating from a rather packed starting configuration, but eventually found that:
fix rigid/smallcan integrate the simulation from relatively unfavourable starting configurations (PACKMOL-packed, so there weren’t any initial clashes < 2 Angstrom), but slows down unbelievably badly.
fix shake with the usual very small virtual particle mass crashes outright on any reasonable timestep, even down to 0.01 fs. (Frankly as a practitioner I much rather my system crash outright, than slow to a trickle without much feedback.)
Running fix shake with a massive virtual particle (8 amu transferred from the O to the charge site) allows fast equilibration of the worst clashes and subsequently fix rigid/small runs smoothly.
I haven’t tried fix shake on the “cured” configurations, but now I’m curious.
Additionally – some of the solution-finding routines in fix shake can be analytically improved, although I’d probably do so in derived fixes for compatibility with old runs.
How would you keep the M point in position with fix shake without violating its restrictions?
For a three point water, basically you end up with three bond constraints (2x O-H, and H-H which is converted from the angle constraint). If you add the M point you would have to add a bond from M to each of the three atoms and constrain its length. This would require a new constraint sub-function to solve the constraint equation for six bonds.
@Danilo_Borges while this isn’t exactly your field of research, you may be interested in the following prior implementations of off-center charges for modelling polarisable surfaces:
Full disclosure: I was not a co-author on any of these papers, but I work in the general field of MD with fast polarisable metallic surfaces and any attention I can bring to the field makes me happy
LAMMPS allows for the following kinds of clusters to be constrained: one central atom bonded to 1 or 2 or 3 atoms, or one central atom bonded to 2 others and the angle between the three atoms also constrained.
and assumed that a three-bond structure could also have three angles constrained (which it can’t, not with the current solver).
Where the arguments of the style coul/long/offcentre (and coul/cut/offcentre) are:
cutoff, as in the regular Coulomb pair style.
nsites Types of off-centre charges in the system.
type molFrameSiteX molFrameSiteY molFrameSiteY molFrameCharge: atom-type, XYZ location, and charge, in the molecular frame of reference. This sequence must be repeated nsites times.
The same extra parameters are needed for kspace_style pppm/offcentre. In the example above, three charges are defined for a spheroid with atom type 1.
The original code has since been morphed into a new pair style incorporating the Gay-Berne potential. The original code is still available in this fork of LAMMPS, version 30Oct19. You are welcome to port the original code to the latest version of LAMMPS, or use the old one as is.
For @srtee: we added off-centre charges to finite-size spheres to create a kind of coarse-grained model of TIP4P water. The results are published here. In a nutshell: the force field is exactly the same as the corresponding atomistic model of TIP4P. The torques are slightly different, as the charges act on the centre of the spheroid instead of the COM. Also, the inertia moment of the CG-TIP4P is bigger than that of the atomistic model, resulting in slower dynamics.
