center of mass with periodic boundary conditions

Dear all,

I have a question about the calculation of center of mass with
periodic boundary conditions. I have seen in the code that lammps uses
unwrapped atom coordinates to perform the calculation. There seems to
be an alternative method discussed here:

http://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions

This can be easily computed in lammps, using something like that:

compute 1 all property/atom xs
variable xi atom lx/2.0/PI*cos(c_1*2.0*PI)
variable zeta atom lx/2.0/PI*sin(c_1*2.0*PI)
compute 2 all reduce ave v_xi v_zeta
variable xb equal lx/2.0/PI*(atan2(-c_2[2],-c_2[1])+PI)

I'm trying to consider a bubble in a liquid with periodic boundary
conditions (so really it is an array of bubbles), and if the liquid
flows relatively to the bubble, then it seems to me that the two
definitions of the center of mass will differ. I would like to apply
fix recenter to the system in order to maintain the bubble at a given
position, therefore using the second definition for the center of
mass, i.e. not the one used internally by lammps. What would be the
best (simplest) way to do that?

Best regards,
Laurent

Dear all,

I have a question about the calculation of center of mass with
periodic boundary conditions. I have seen in the code that lammps uses
unwrapped atom coordinates to perform the calculation. There seems to
be an alternative method discussed here:

http://en.wikipedia.org/wiki/Center_of_mass#Systems_with_periodic_boundary_conditions

This can be easily computed in lammps, using something like that:

compute 1 all property/atom xs
variable xi atom lx/2.0/PI*cos(c_1*2.0*PI)
variable zeta atom lx/2.0/PI*sin(c_1*2.0*PI)
compute 2 all reduce ave v_xi v_zeta
variable xb equal lx/2.0/PI*(atan2(-c_2[2],-c_2[1])+PI)

I'm trying to consider a bubble in a liquid with periodic boundary
conditions (so really it is an array of bubbles), and if the liquid
flows relatively to the bubble, then it seems to me that the two
definitions of the center of mass will differ. I would like to apply
fix recenter to the system in order to maintain the bubble at a given
position, therefore using the second definition for the center of
mass, i.e. not the one used internally by lammps. What would be the
best (simplest) way to do that?

hmm... this is a bit of a tricky issue.

after having thought about this for a bit, i have been wondering how
to approach this, also in a more general way, so that it would be
useful beyond your case. for starters, i would try to devise a method
to keep track of the bubble rather than infer it from the center of
mass of the rest of the system. my suggestion is the following:

put suitably sized grid for an occupancy map through the simulation
box. loop over all atoms and set all voxels that contain one or more
atoms to 1 otherwise they remain at 0. then loop over all voxels and
label consecutive zones or occupied of not occupied voxels with unique
numbers. with periodic boundaries and in running parallel, one would
then have to figure out, if such zones would need to be merged.
finally, you can do a final loop over all voxels collecting the center
of each zone and output it as a local compute. this kind of approach
would also work for detecting/quantifying clustering in self-assembly
simulations

this could be modified in a number of ways, e.g. rather than
collecting a simple occupancy (which would be efficient with STL
vectors, that can store only 0/1 maps in single bits), one could add a
weighted gaussian where width and height could be assigned by atom
type/mass and then select occupied/unoccupied voxels via a threshold
and both identify the minimum/maximum of each cluster/hole or compute
a geometrical or weighted center.

it would require a bit of programming though... not just a quick hack.

axel.

Have you considered the voronoi compute to track the bubble? I am no expert here but if the bubble size is much larges than the average distance separation between fluid units perhaps the voronoi compute can aid you construct and track the bubble somehow.
Carlos

Thank you both for your suggestions, I will try to investigate in
these directions.

Axel's approach indeed involves a bit of programming, and for
quantifying clustering I think the already implemented compute
cluster/atom could be used (alas it cannot be applied to bubbles). In
fact I had seen a similar approach in this paper:

http://pubs.acs.org/doi/abs/10.1021/jp807727p

but here again it is not trivial to implement.

Carlos, are you suggesting to track the atoms at the surface of the
bubble using compute voronoi/atom? I'm afraid that for large bubbles
there would be problems with the ghost atom cut-off, what do you
think? Anyway, it seems to me that surface atoms could also be tracked
using e.g. coord/atom, centro/atom or pe/atom.

I will also try the tracking through following the liquid center of mass...

Once / if I manage to keep track of the bubble position, what would be
the simplest way to adapt fix recenter in order to recenter the system
around the bubble position instead of the center of mass of some
group? Even if I identify the position of the bubble with the center
of mass of surface atoms, I cannot use fix recenter directly since it
uses unwraped atom coordinates to compute the com.

Best regards,
Laurent

2014-02-27 17:13 GMT+00:00 Carlos Campana <[email protected]...>:

Carlos, are you suggesting to track the atoms at the surface of the
bubble using compute voronoi/atom? I'm afraid that for large bubbles
there would be problems with the ghost atom cut-off, what do you
think? Anyway, it seems to me that surface atoms could also be tracked
using e.g. coord/atom, centro/atom or pe/atom.

Like I said I am no expert doing this stuff. Yet, (at least for a single
bubble problem) I would expect that if you compute the voronoi volume for
each atom (real atoms not ghost) the size distribution will show some atoms
having large volumes (assuming the bubble is big enough with respect to the
average inter-particle distance in the fluid) while the rest will not. The
atoms at the large size end of the distribution should pretty much be
located at the bubble surface and they may be well be a way to track the
bubble.

NB: Never used the voronoi package thus it may have some other features you
can use to better tune the localization of the bubble. In combination with
other techniques it also may help you track multiples bubbles.
Hopefully someone else more versed in this type of calculation will reply
with additional suggestions.

Carlos

Once / if I manage to keep track of the bubble position, what would be
the simplest way to adapt fix recenter in order to recenter the system
around the bubble position instead of the center of mass of some
group?

Just add a keyword to fix recenter like

com v_xcom v_ycom v_zcom

where 3 variables can be specified that calculate
what you say the current COM should be.

Those could come from the input-script formulas in your

original email. Or if that is an inefficient way to compute
it, you could write a new compute com/oddball that
calculates it, then set 3 variables to the x,y,z values

of the compute.

Do those formulas work for just a small cluster straddling a PBC,

or for a liquid (with a few bubbles) filling the box?

If you want it for the latter, why not just compute the COM
w/out the unwrapping? That would be much simpler than

those formulas. That could be done with a flag on the
compute com command.

Steve

[...]

If you want it for the latter, why not just compute the COM
w/out the unwrapping? That would be much simpler than
those formulas. That could be done with a flag on the
compute com command.

there is one problem with that. every time an atom is wrapped around
the PBC, your COM does a significant, system size dependent jump. i
don't think that is very desirable.

axel.

good point, although if you have 10K fluid atoms

it’s probably not a big deal.

Is the other method continuous when an atom crosses,
for any number of atoms?

Steve