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.