[lammps-users] Basic Monte Carlo moves

Dear lammps-users,

I have a system of particles with long range Coulombic interactions.
I am using LAMMPS to model the dynamics of this system.
However, for the purpose of computing the equilibrium state, I would like to add some simple Monte Carlo moves.
Specifically, when the locations of the particles are restricted to lattice sites, I would like to be able to swap arbitrary particles "i" and "j" based on Boltzmann criteria. In this case I no longer need accurate dynamics of the system.

What would be the easiest way to implement this? Add a fix in lammps?
Use SPPARKS (which seems to be missing long range Coulomb calculations)?

Best,

Nickolay Shestopalov

Dear lammps-users,

dear nikolay,

I have a system of particles with long range Coulombic interactions.
I am using LAMMPS to model the dynamics of this system.
However, for the purpose of computing the equilibrium state, I would
like to add some simple Monte Carlo moves.
Specifically, when the locations of the particles are restricted to
lattice sites, I would like to be able to swap arbitrary particles "i"
and "j" based on Boltzmann criteria. In this case I no longer need
accurate dynamics of the system.

What would be the easiest way to implement this? Add a fix in lammps?

yes, a fix is what i would try first.

however, there is one complication: you will have to check how you
handle the communication of the new coordinates of the two swapped
particles to the neighboring cells, or how you handle the case that one
atom is a local and the other is a ghost atom.

cheers,
     axel.

Would it be easier to use SPPARKS and make it use the kspace package from LAMMPS?

Thanks,
Nikolay

Axel Kohlmeyer wrote:

Dear Axel,

Thanks for the advice. I imagine the problem of communicating the new coordinates that you mentioned is a common issue in LAMMPS. Could I use some other source code from LAMMSP as a template?

Best,
Nikolay

Axel Kohlmeyer wrote:

Dear Axel,

dear nicolay,

Thanks for the advice. I imagine the problem of communicating the new
coordinates that you mentioned is a common issue in LAMMPS. Could I
use some other source code from LAMMSP as a template?

not really. this is a special case due to the monte carlo moves
that you are planning to make. LAMMPS is using Verlet integration
and advanced Verlet lists to have (near) linear scaling with system
size. during regular MD the code keeps track of how much atoms
move, and once they hit a safety margin, the a rebuild of the
neighbor lists is triggered. with your atom swaps, you have
to enforce this neighbor list rebuild. you also have to keep in
mind that you can only swap atoms that are contained within
the local and ghost atoms for a given MPI task. so if you want
to swap atoms that are fairly far apart, you may have to play
with the neighborlist build flags to have them contained on
the corresponding nodes.

the most similar "fix"es in that respect should be "fix deposit"
and "fix evaporate", so i would check what they do to trigger
a rebuild of the neighbor list when atoms are added or removed.

cheers,
    axel.

Do it in LAMMPS, not SPPARKS. To swap arbitrary
particles will require point-to-point communication, which
most options in LAMMPS don't do. But there are MPI
calls to do it - it's not difficult.

Steve

Axel, Steve, thank you for your suggestions.

I have another question. If during the simulation I use the “set” command that reassigns the charge to a particle, will that still pose the problem with the rebuilds you mentioned? In other words, I don’t need to move the pair of particles physically, I just would like to swap the charges on them.

  • Nickolay

Axel Kohlmeyer wrote:

Axel, Steve, thank you for your suggestions.

I have another question. If during the simulation I use the "set"
command that reassigns the charge to a particle, will that still pose
the problem with the rebuilds you mentioned? In other words, I don't

you cannot avoid having to communicate the information around,
regardless of whether you swap coordinates or charges.

the set command should do the communication for you, though.
not sure whether this is a better solution than triggering
the neighborlist rebuild from a fix. also, it would be easier
to track those changes in post processing via tracking the
atom id rather than having to dump the charges, too.

need to move the pair of particles physically, I just would like to
swap the charges on them.

i don't quite get what your concern is. if it is the programming,
the changes to the source code would be rather small. if it is
the communication, you cannot avoid it (except when not running
in parallel).

i am overlooking anything?

axel.

If you want to select two random particles and swap their
charges, then you'll have to write a fix to do precisely that.
The set command doesn't know about choosing two particles
or swapping their charge. I don't see where neighbor
lists are an issue. Charges aren't stored in a neighbor list.

Steve

Yes, the code I wrote swaps charges on the particles and the calls lammps that was built as a library. There is one small technicality though.

At this point I only need to switch neighboring particles. I choose one particle at random and then access its neighbor list. Sometimes the neighboring atoms have tags that are larger then the total number of atoms. Is this related to the ghost atoms? How do I obtain the actual tag of that atom?

The simulation uses periodic b.c.
So far I've been only running on a single processor.

Thanks,
Nickolay

Steve Plimpton wrote:

Yes, the code I wrote swaps charges on the particles and the calls lammps
that was built as a library. There is one small technicality though.

At this point I only need to switch neighboring particles. I choose one
particle at random and then access its neighbor list. Sometimes the
neighboring atoms have tags that are larger then the total number of atoms.
Is this related to the ghost atoms? How do I obtain the actual tag of that
atom?

no. if the neighbor list entry is larger than "nall" than
this is an "excluded" interaction for which a scaling
factor has to be applied (typically from 0.0 to 1.0).
the fraction j/nall is the index in the scaling factor
array and the modulus j % nall is the real atom index.
this can be easily inferred from looking at the compute()
method in e.g. pair_lj_cut.cpp or pair_lj_cut_coul_cut.cpp.

the existence of those topology dependent exclusions is
one of the reasons why you _have_ to enforce the
re-neighboring after a swap.

cheers,
    axel.

What Axel describes is true, but only applies if you have a system
with bonds and excluded pairwise interactions (via the special_bonds
command). If you have just an atomic system (no bonds),
then the neighbor list will still store indices of atoms (not tags which
are atom IDs), larger than the number of local owned atoms (nlocal).
This is b/c they refer to ghost atoms. You do not want to change
a property (e.g. the charge) of a ghost atom, b/c it is non-permanent.
You would only want to change the property of the real, owned atom.

Steve

Please correct me if I have wrong picture in mind:

If, for instance, the particle is close to the boundary, it would have ghost atoms representing atoms in the near image, correct? I still would like, however, to swap that particle with the real particle on the other side of the box. Is there any way to infer its index from the index of the ghost particle?

- Nickolay

Steve Plimpton wrote:

All atoms (ghost and real) have IDs. The ID for
the ghost atom is the same as its real counterpart.
The problem with what you ask is that the processor
doing the swapping would have to communicate
with the processor owning the real version of
the ghost atom.

Steve

I am confused with this description since I have been only using a single
processor.

- Nickolay

All atoms (ghost and real) have IDs. The ID for
the ghost atom is the same as its real counterpart.
The problem with what you ask is that the processor
doing the swapping would have to communicate
with the processor owning the real version of
the ghost atom.

Steve

Please correct me if I have wrong picture in mind:

If, for instance, the particle is close to the boundary, it would have
ghost
atoms representing atoms in the near image, correct? I still would

like,

however, to swap that particle with the real particle on the other side
of
the box. Is there any way to infer its index from the index of the

ghost

particle?

- Nickolay

Steve Plimpton wrote:

What Axel describes is true, but only applies if you have a system
with bonds and excluded pairwise interactions (via the special_bonds
command). If you have just an atomic system (no bonds),
then the neighbor list will still store indices of atoms (not tags

which

are atom IDs), larger than the number of local owned atoms (nlocal).
This is b/c they refer to ghost atoms. You do not want to change
a property (e.g. the charge) of a ghost atom, b/c it is non-permanent.
You would only want to change the property of the real, owned atom.

Steve

Yes, the code I wrote swaps charges on the particles and the calls
lammps
that was built as a library. There is one small technicality though.

At this point I only need to switch neighboring particles. I choose
one
particle at random and then access its neighbor list. Sometimes the
neighboring atoms have tags that are larger then the total number of
atoms.
Is this related to the ghost atoms? How do I obtain the actual tag

of

I am confused with this description since I have been only using a single
processor.

yes, but since LAMMPS is from the ground up designed as a parallel
program with a distributed memory parallelization scheme, your code
will be practically useless, if you do not factor this in. those communications
will become no-ops in the single processor case, since all atoms are
owned by that processor.

cheers,
    axel.

There are still ghost atoms on a single proc, but you can
always find the real atom corresponding to a ghost, e.g.
by using the atom->map() function. So I'm not sure
what your question is.

Steve