# LAMMPS to build a magnetic nanoparticle molecular dynamic simulation?

Hello

I am new to LAMMPS and was hoping I could ask a calibration question. Basically, I would like to know if LAMMPS can be used to build what I need + how heavy of a lift such a build might be (I am experienced with Matlab, and I don’t have any aversion to learning new syntax).

I would like to create a molecular dynamic simulation of interacting magnetic nanoparticles moving through a barrier.
(The possibly unusual part is that there is a coupling between the particles, the magnetization of one affects magnetization of another, item 3## below.)

The features of the simulation I need are:

1. I would track N magnetic nanoparticles, each particle a rigid spherical body (e.g. 100 nm in size). So the state vector would be R = [R1;R2;…;RN] where Rj=(xj,yj,zj,wx,wy,wz), the 3D location and 3D orientation for the jth magnetic paricle.

2. Each particle contains Kj iron oxide cores inside it. An iron oxide core is usually about 5 nm in diameter. So each rigid sphere contains a random placement of iron cores inside it, randomly placed. As the sphere translates and turns, the cores are displaced just by rigid body translation and rotation of that nanoparticle sphere.

3. There are magnetic forces on each on iron oxide core. Two types of magnetic forces: a) an force on each core due to the applied external field; b) dipole-dipole magnetic interactions between the cores. (##) Part b creates a nonlinear time-varying magnetic coupling between the cores. If core p is magnetized, that creates a surrounding magnetic field that can change the magnetization of cores q,r,s. The end result is that at each moment in time there is a linear matrix equation ( [ A® ] M = B) that needs to be solved to find the magnetizations of all the core (M = [m1x,m1y,m1z, m2x,m2y,m2z,…] ). Once the magnetizations are known, then the force on each core can be calculated. Once force are known, rotation is just sum of the torques on each sphere.

4. The applied magnetic field and magnetic gradient is prescribed and is simple, is just set by an analytical formula. I don’t need to solve Maxwell’s equations.

5. The magnetic particles will have diffusion and steric hindrance (so they don’t overlap each other). It’s ok to make the steric hindrance smoothed out over a long length scale, e.g. over 20% of the particle diameter, so that the governing equations are less stiff. The real important part is in the core to core interactions and how the forces that creates causes the magnetic particle to aggregate and chain up. If they aggregate with 20 nm gaps instead of 1 nm gaps between them, that’s ok (and I can correct for this by choosing particle size and core loading density appropriately).

6. There will be viscous drag on the particles, both translational and rotational drag. It’s ok to keep this part simple (there are modifications to Stokes drag when particles are nearby, but it is not the central feature I am trying to simulate so is ok to keep this part as simple as possible.)

7. If helpful, it is ok to ignore inertia (both translational and rotational inertia). These are nanoparticles with very tiny mass. So m x’’ + mu x’ = Fmagnetic + Fsteric + dW could be reduced to x’ = (Fmagnetic + Fsteric + dW)/mu where m is the mass, mu is the viscosity coefficient, ’ denote d/dt, and dW is Brownian noise. Whether LAMMPS allows a drop from second to first order dynamics, I don’t know, but it’s ok from a physics point of view.

8. There will be a barrier that the particles have to cross. The key issue is to better understand the physics of that crossing. The barrier could be a) a region of increased viscosity, or b) pores (meaning steric hindrance at a wall except for gaps in that wall). I need to see how the particles cross that barrier as the magnetic field moves and aggregates them.

If at all possible, please let me know if this is something that is sensible/feasible to set up in LAMMPS. I don’t know how easy/hard it is to customize forces in LAMMPS. I would have to set up the kind of rigid-body multi-core particles I am interested in. I am especially concerned about how to put in the core-to-core magnetic coupling mentioned in item 3## above because that involves solving an additional (linear at each moment) system of equations at each time step within LAMMPS.

It is my first time posting to this ListServe. Any thoughts would be much appreciated.

My email is BenjaminShapiro.Shapiro@…24…

Thank you
Ben

A quick comment: I think magnetic dipole interactions are not your biggest concern since that problem has been solved for electric dipoles already. As long as there is no nett charge it should be doable. However, it would be wise to incrementally build up such a complex system and check and double check at each expansion.

I would like to create a molecular dynamic simulation of interacting magnetic nanoparticles moving through a barrier.
(The possibly unusual part is that there is a coupling between the particles, the magnetization of one affects magnetization of another, item 3## below.)
The features of the simulation I need are:

1. I would track N magnetic nanoparticles, each particle a rigid spherical body (e.g. 100 nm in size). So the state vector would be R = [R1;R2;…;RN] where Rj=(xj,yj,zj,wx,wy,wz), the 3D location and 3D orientation for the jth magnetic paricle.

(see the comment below regarding “fix rigid” and DATA files.)

1. Each particle contains Kj iron oxide cores inside it. An iron oxide core is usually about 5 nm in diameter. So each rigid sphere contains a random placement of iron cores inside it, randomly placed. As the sphere translates and turns, the cores are displaced just by rigid body translation and rotation of that nanoparticle sphere.

There are many different versions of “fix rigid … molecule” which should do what you want.
Create a data file and assign each rigid sphere a different molecule-ID number
https://lammps.sandia.gov/doc/fix_rigid.html

As for how to create DATA files, you can either write a program to build them yourself,
https://lammps.sandia.gov/doc/2001/data_format.html
https://lammps.sandia.gov/doc/atom_style.html
…There is also a list of DATA-file builder tools at
https://lammps.sandia.gov/prepost.html
Either moltemplate or topotools should be general enough to build these kinds of particles. (Other molecule builders tend to be more specialized.) Topotools uses TCL. Moltemplate uses its own “object-oriented” style language.

1. There are magnetic forces on each on iron oxide core. Two types of magnetic forces:
a) an force on each core due to the applied external field;

First, I found the following PDF file which looks vaguely relevant, but what they did is not as sophisticated as what you are trying to do:
http://www.rsc.org/suppdata/sm/c0/c0sm01424a/c0s

The closest thing that LAMMPS currently has to implementing an external B field is “fix addforce”, but I don’t think you can use it to apply torque to a point dipole. However if you represent each dipole by constructing it from two discrete particles, then it should be possible. Take a look at these (related) web pages:

https://lammps.sandia.gov/doc/variable.html

If that is not flexible enough for your needs, it sounds like you will have to write your own “fix”. In LAMMPS, a “fix” refers to some code which modifies the equations of motions of the atoms.

Incidentally, you might not have to do this. Tim Sirk posted that he had written a “fix bfield” a couple years ago. He still checks the mailing list, but your more likely to get a reply from him if you mail him directly.

A (slopily written) introduction to writing fix was posted a few days ago
https://sourceforge.net/p/lammps/mailman/message/36443505/

…and you should also read the official lammps documentation for customizing the code here:
https://lammps.sandia.gov/doc/Developer.pdf

b) dipole-dipole magnetic interactions between the cores.

Fixed point dipole-dipole interactions are already implemented in lammps:
https://lammps.sandia.gov/doc/pair_dipole.html
If you want to use point dipoles, then this will require that you use “atom_style hybrid full dipole” or something similar:
https://lammps.sandia.gov/doc/atom_style.html

…however, if you need polarizable dipoles, then you cannot (currently) use point dipoles, so please ignore what I just said (and the links above).

(##) Part b creates a nonlinear time-varying magnetic coupling between the cores.
If core p is magnetized, that creates a surrounding magnetic field that can change
the magnetization of cores q,r,s. The end result is that at each moment in time
there is a linear matrix equation ( [ A® ] M = B) that needs to be solved to find
the magnetizations of all the core (M = [m1x,m1y,m1z, m2x,m2y,m2z,…] ).
Once the magnetizations are known, then the force on each core can be
calculated. Once force are known, rotation is just sum of the torques on
each sphere.

An approach similar to the one used for polarizable electric dipoles could be used. (See below.) But I suspect this will need to be modified in a non-trivial way if the matrix is not diagonal. So it sounds like you are going to have to write your own fix to implement this feature.

I suggest you copy and modify the code for “fix drude”. You can obtain this code by downloading the lammps code, cd to the “src/” subdirectory, and typing “make yes-user-drude”. This will copy the “fix_drude.cpp” and “fix_drude.h” files into the “src/” directory where you can copy and edit them. Your fix will probably be more complicated. Anyway, here is the documentation for “fix drude”

1. The applied magnetic field and magnetic gradient is prescribed and is simple, is just set by an analytical formula. I don’t need to solve Maxwell’s equations.

Good.

1. The magnetic particles will have diffusion and steric hindrance (so they don’t overlap each other). It’s ok to make the steric hindrance smoothed out over a long length scale, e.g. over 20% of the particle diameter, so that the governing equations are less stiff. The real important part is in the core to core interactions and how the forces that creates causes the magnetic particle to aggregate and chain up. If they aggregate with 20 nm gaps instead of 1 nm gaps between them, that’s ok (and I can correct for this by choosing particle size and core loading density appropriately).

You can use Lennard-Jones particles with large “sigma” parameters to implement as much steric hindrance as you need.
If you need to have the Lennard-Jones barrier centered exactly half-way between the two particles in the dipole (assuming you are approximating dipoles as two particles), then you can add a third particle between the other two and connect all three together with (stiff?) springs.

If you decide you don’t need polarizability and are using point-dipoles, then you can use “pair_style hybrid/overlay dipole lj/cut …” (or something similar):

https://lammps.sandia.gov/doc/pair_hybrid.html

1. There will be viscous drag on the particles, both translational and rotational drag. It’s ok to keep this part simple (there are modifications to Stokes drag when particles are nearby, but it is not the central feature I am trying to simulate so is ok to keep this part as simple as possible.)

There are many ways to add viscous damping and thermostatting to LAMMPS. You can even add hydrodynamic interactions.

1. If helpful, it is ok to ignore inertia (both translational and rotational inertia). These are nanoparticles with very tiny mass. So m x’’ + mu x’ = Fmagnetic + Fsteric + dW could be reduced to x’ = (Fmagnetic + Fsteric + dW)/mu where m is the mass, mu is the viscosity coefficient, ’ denote d/dt, and dW is Brownian noise. Whether LAMMPS allows a drop from second to first order dynamics, I don’t know, but it’s ok from a physics point of view.

For what it’s worth, in my own experience, I was able to get simulations to move much more quickly and efficiently using second order dynamics (with fix langevin) than writing my own code that implemented Brownian dynamics. I ranted about this topic here:
For the timescales you probably care about, it does not matter if you use fix langevin.
Docs below:

https://lammps.sandia.gov/doc/fix_langevin.html

https://lammps.sandia.gov/doc/fix_nve.html

The following web sites are probably not relevant, but I mention them anyway:
https://lammps.sandia.gov/doc/pair_dpd.html

https://lammps.sandia.gov/workshops/Aug13/Denniston/Denniston_LAMMPS2013.pdf

https://lammps.sandia.gov/doc/fix_lb_fluid.html

https://lammps.sandia.gov/doc/PDF/SPH_LAMMPS_userguide.pdf

http://web.math.ucsb.edu/~atzberg/pmwiki_intranet/index.php?n=SoftwareMango.Homepage?setskin=mango_selm_Layout1

https://www.research-collection.ethz.ch/bitstream/handle/20.500.11850/153451/eth-5603-01.pdf

1. There will be a barrier that the particles have to cross. The key issue is to better understand the physics of that crossing. The barrier could be a) a region of increased viscosity, or b) pores (meaning steric hindrance at a wall except for gaps in that wall). I need to see how the particles cross that barrier as the magnetic field moves and aggregates them.

You could either use “fix addforce” (as I mentioned above)

https://lammps.sandia.gov/doc/variable.html

or create a wall by placing atoms at locations at specific locations in the barrier, and preventing them from moving. For an example of this, see:

http://moltemplate.org/visual_examples.html#translocation

http://moltemplate.org/visual_examples.html#nanotube+water

(If you don’t use “moltemplate” then you can ignore the “.lt” files and just look at the “run” files in these examples. The “run.in.nve” and “run.in.npt” files for both examples have the instructions for specifying groups of atoms which you wish to remain immobile.)

If at all possible, please let me know if this is something that is sensible/feasible to set up in LAMMPS. I don’t know how easy/hard it is to customize forces in LAMMPS. I would have to set up the kind of rigid-body multi-core particles I am interested in. I am especially concerned about how to put in the core-to-core magnetic coupling mentioned in item 3## above because that involves solving an additional (linear at each moment) system of equations at each time step within LAMMPS.

It is my first time posting to this ListServe. Any thoughts would be much appreciated.

My email is BenjaminShapiro.Shapiro@…24…

Thank you
Ben

I don’t know if this helps.
What you are doing is probably possible in LAMMPS if you are willing to write some C++ code (and/or email Tim Sirk, see link above). In that case, it’s probably easier to start with LAMMPS than writing your own MD engine from scratch. (Gone are those days, hopefully forever.)
Good luck.

Andrew

In my last post, I mentioned embedding masses connected by springs
into a rigid spherical nanoparticle.

...I forgot to mention how you might implement this. It's okay to
attach mobile particles (via bonds) to other particles which are part
of a larger rigid structure. When you use "fix rigid" you have an
opportunity to specify what group of particles you want to include.
Just omit the particles attached to the springs from this group.
https://lammps.sandia.gov/doc/group.html

But it sounds like you will need to add some way to restrain their
motion to preferentially lie along some vector relative to the
permanent dipole by creating your own fix for this purpose. (Perhaps
some customized version of "fix drude", as I mentioned earlier.)

Another thing to consider is periodicity:
If you are running a condensed-phase periodic simulation, you may run
into trouble because the periodic long-range electrostatics kspace
solver used by LAMMPS considers the electric charge of each atom. I
don't think there's an easy way to override that to ignore the
electric charge and substitute some kind of magnetic monopole or
dipole there. If you don't have both electric charges and magnetic
dipoles (or if you don't need to run a periodic simulation), then I
guess this is not an issue.

Have fun

Andrew

While not directly applicable to rigid bodies with a magnetic moment, there is plenty to build up from in the SPIN package for your needs. Please refer to the links below about a package in LAMMPS for simulating per-atom magnetic spins.

https://lammps.sandia.gov/doc/Howto_spins.html

https://www.sciencedirect.com/science/article/pii/S0021999118304200

Long range interactions of magnetic spins is under current development here at Sandia, it will be added soon (if not already). If you are eager to adapt this code for your needs, please let me or Julien Tranchida (jtranch@…3…) know and we may be able to help.

-Mitch Wood

P.S. Thanks Andrew for a detail response, these are great to see on the mailing list.

Hi Ben,

According to your requirements in modeling magnetic particle, the work below may to some extent help you. This work implements a new fix to handle the magnetic force field exerting on the nanoparticle based on the existed fix efield which solve the electric field prolem. So, my suggestion is to write a new fix for the magnetic field and integrate the motion of nanoparticle using fix rigid like class.