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.