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.
I was able to write a fix module using Claude and it appears to work for timestep 0. The calculation blows up after that. I’ve been trying to use the pair_style hybrid command to combine the lj/cut/coul/long and spin/dipole/long pair styles without having to modify code, but I can’t get the code to evaluate the long range part of the interactions for both the coulomb and spin dipole interactions. For starters I’ve just been using the Ewald sum. I’ve put print statements in the compute commands and all I’m seeing is the pair calculations for both pair styles and the ewald calculation for the spin dipole. There doesn’t seem to be a corresponding hybrid option for the kspace_style command. Is there a way to do this using the existing command interface or will I need to create a new kspace module that combines both calculations?
Isn’t that inconsistent? Shouldn’t there be dipole - charge interactions, too?
There cannot be. Kspace is a “everything with everything” kind of interaction and thus cannot be partitioned like pair-wise interactions. That means, that for the pair-wise interactions that include the short-range part of a long-range interaction, those parts have to be consistent, i.e. you cannot have two different kinds.
There are two separate interactions that don’t mix. The first is point coulomb charge interactions that are part of the lj/cut/coul/long interactions and the second are magnetic spin dipole interactions that are part of the spin/dipole/long interactions. It should be possible to calculate the long range part of these separately and then add them together at the end, although it might be more efficient to calculate them both simultaneously.
Of course Axel, I was following the conversation and planning to comment at some point ;-).
I agree with you Bruce, but I am not really sure this is possible within the current implementation of kspace. In the meantime, could you try something like this:
The idea is to only use the kspace calculation for the coulombic interactions for now. For the spin-spin dipolar interactions, we could try first with a direct real space calculation only.
One quick question: I am guessing your atoms are defined “hybrid charge spin”. Is it correct? Also, the ‘spin’ atomic style will only work within the “metals” units for now. I understand this might not very very convenient for ‘macro’ spins.
Also, note that adding the interaction between the spins and a magnetic field should be straightforward with:
There is no direct interaction between the point charges defined as part of the lj/cut/coul/long pair style and the magnetic dipole defined by the spin/dipole/long pair style. The cross-terms that would couple those together don’t exist.
Julien,
We are using the atom style command
atom_style. hybrid full spin
Replacing full with charge leads to an error
Bond_style command when no bonds allowed
We are trying to simulate an aqueous solution of paramagnetic ions, so we need bonds and angles to model the water molecules.
To replace the spin/dipole/long pair style with spin/dipole/cut will take a bit of programming, since the current fix command only uses the long ranged dipole interactions. We did try using a short ranged cutoff for the coulomb interactions (lj/cut/coul/cut), but we get an error
Coulomb styles of pair hybrid sub-styles do not match (src/pair_hybrid.cpp:1147)
The calculation does appear to get to the first step if I make everything short ranged (both coulomb and spin dipole) but it blows up immediately. The forces on water at t=0 appear to be extreme (order 10^6-10^7).
I haven’t looked at the precession/spin command, but we were planning on introducing the external field by including it in our fix command. The external field combines with the local field due to the other dipoles in the system to bias the orientation of the dipoles and thereby effect the interactions between paramagnetic ions. We aren’t at all sure this is enough to reproduce the behavior of these solutions, but this is where we are planning on starting.
Ah indeed, you need bonds so that you cannot use simple charge. That being said, ‘hybrid full spin’ should work fine.
The calculation does appear to get to the first step if I make everything short ranged (both coulomb and spin dipole) but it blows up immediately.
It may be surprinzing, but I see this as rather encouraging. At least the calculation seems to be running. There could be a lot of reasons for those very large forces, such as unit issues for example, or cell size / initial configuration of the particles.
Would it be possible, in your script, to include a spin and a lattice minimizations at the beginning of the calculation? For example, by setting something like this:
Another question that could help me understand the potential issues: did you try running a ‘simpler’ calculation with what is already in LAMMPS? Something like this:
atom_style hybrid full spin
# set box
kspace_style ewald 1.0e-4
pair_style hybrid/overlay spin/dipole/cut 20.0 lj/long/coul/long 10.0
# set interactions
min_style cg
fix 2 all box/relax iso 0.0 vmax 0.001
minimize 0.0 1.0e-12 1000 1000
min_style spin
minimize 0.0 1.0e-12 10000 1000
I am sorry to insist Bruce, I understand this is quite far from your goal, but this is already rather ambitious as I don’t believe this has been tested before. And that could help understanding why large forces may occur.
Thanks for the response. I will give these a try. Before I do, I want to verify that I can run the configuration I have using just the lennard-jones and coulomb interactions and leaving the spins out. I’m getting the same blowup when I do that so I need to see if there is some problem with how I’m setting the SHAKE constraints or whether or not there might be an issue with the units. Our initial configuration comes from an equilibrated simulation of an ionic solution, so it should run if you leave out the spins. Once I get that sorted out, I will try your suggestions.
I was interested to learn that you have a CG routine that mimizes the spin configuration, Once I get some other things working, it probably would be worth seeing if I can replace the iterative solver in our fix routine with the CG minimization.
MBX includes CG-minimised dipoles interacting with Coulombic charges, including monopole-dipole concentric sites. They use Thole damping on all dipoles, which seems sensible to me – since dipoles have extra reciprocal scaling, dipoles strong enough to influence macro structure must be exceptionally susceptible to over-polarisation without short-range shielding.
You may be able to get lots of bang for your buck if you can somehow specify your paramagnetic ions in their force field format and thus utilise their framework!
We tried the coordinate mimimization and that seemed to fix the problem of the simulation blowing up on the first step. The minimization of the spins seems to be coupled to using the nvt/spin algorithm for the time integration (it looks like it wants to set a time step for some reason). We don’t want to add additional time dynamics to the simulation beyond minimizing the spins before each time step (at least for the time being).
The fact that the configuration blows up on the first step is a little confusing since we are starting from an equilibrated configuration. If we just use “atom_style full” and “pair_style lj/cut/coul/long”, along with the ewald sum for long range interactions, the simulation starts up no problem but if I add the “hybrid” keyword to the atom_style and just use full, it blows up and the problem persists if I try adding the spin interactions.
I’m currently debugging our fix command, which doesn’t look like it is actually modifying the spin interactions, but I did notice something that looks like a bug in ewald_dipole_spin.cpp file in directory KSPACE.
In line 451, it looks like the last index on tk[i] should be a 2 instead of a 3.
Shern,
I will take a look at the MBX code but it sounds like their modeling is a bit orthogonal to what we are doing. The coulombic charges and magnetic dipoles are not interacting, so there is no dipole-point charge interaction to optimize. We would be interested in optimizing the spin orientation calculation so there might be something to look at there.