We are interested in adding some simulation capability to LAMMPS to model aqueous solutions of paramagnetic ions using conventional point charges for the electrostatic parts and magnetic point dipoles for the spins. I’ve been looking through the LAMMPS code base and it appears that a lot of the elements for this calculation are available but not necessarily available simultaneously in a single calculation. Also, we are interested in using relaxed configurations for the magnetic dipoles, which will require iterating over the dipole configurations until they converge to a stable configuration. However, the existing force models all seem to be single pass calculations, where the short ranged, bonding and long ranged interactions are all calculated once per force calculation (at least, that appears to be what is happening in the verlet.cpp file).
A couple of questions:
Are there some examples of iteratively calculated force fields in LAMMPS? The obvious choices like polarizable atoms all seem to use methods that can evaluate the force due to polarization in a single pass but also may introduce extra time scales into the calculation.
The existing code seems to support only one Pair model and one KSpace model. Including both coulomb and spin interactions would require two models. Is this likely to be a big problem?
One solution would be to create an iterative force class, similar to KSpace and Pair and bury the spin calculation in that. Is that likely to work or would it require touching too many other parts of the code to be worthwhile?
There is charge equilibration, which is part of ReaxFF and the QEQ package.
The core-shell and Drude oscillator approaches use an extended Lagangian scheme where you have a (decoupled) fictitious dynamic to make the polarization follow the movement and ideally oscillate around the converged ground state so that errors largely cancel.
Pair styles can be combined with pair styles hybrid and hybrid/overlay. It has not been necessary so far, but you could consider doing something similar for your case or just implement a single style that does both calculations. Same for the pair style, e.g. you can use pair style lj/cut/coul/cut or pair style hybrid/overlat lj/cut … coul/cut … and the outcome is the same.
That is impossible to predict. I would be more concerned about time integration. The derivation used in LAMMPS assumes that forces are only dependent on positions and not on velocities, but that may not be a valid in your case.
You may want to look at the ELECTRODE package which does charge equilibration. In our implementation, the charge equilibration fix does an “equilibration pass” during pre_force which sets the desired charges. Subsequently, those charges are then used to calculate Coulombic forces.
I would imagine that a similar approach is possible in your case: write a fix to relax the magnetic spins during pre_force, and then calculate the force using those relaxed spins.
This looks very useful. Pre-calculating the dipole orientation and then just using a conventional dipole calculation for the forces looks like a good approach.