Fixing COM of each molecule, 4000 molecules in total

Dear all,

I’d like to hold the COM’s of about 4000 molecules stationary. It seems that
the only way to do this would be to define 4000 fix com’s, 1 for each molecule.

I was wondering if there’s another way of going about this. Or should I consider writing my own fix/modifying existing fix.

Regards,
Anirban

Dear all,

I'd like to hold the COM's of about 4000 molecules stationary. It seems
that

​what kind of molecules? and what is the physical justification for such a
(rather unusual) request​?

the only way to do this would be to define 4000 fix com's, 1 for each
molecule.

​i don't think that would be practical. it will be a performance disaster,
too.​

I was wondering if there's another way of going about this. Or should I
consider writing my own fix/modifying existing fix.

​it might be possible to use the */chunk features in LAMMPS:
- group the atoms into chunks by molecule id.
- compute the force on the center of each chunk through fix ave/chunk
- convert this into a per atom force that can be added to each atom with
fix addforce. this may require some LAMMPS script trickery, e.g. using a
loop.
​- if there is no initial velocity, then only forces​ inducing rotation
should remain.
- to protect against numerical noise, you may want to subtract the center
of mass velocity from each chunk regularly, too.

of course, you could just write your own integrator that disregards
per-chunk translation.

axel.

Dear Prof. Kohlmeyer,

I have a molecular crystal system with a single partial edge dislocation in it.I’m defining a cylindrical region around the dislocation line (which is along z) which is flexible, while the region
outside the cylindrical region is FIXED. Now the z direction is periodic.

I’m running Umb. Samp. simulations while holding the dislocation core in place using a bias potential,
and I’m holding the boundary region fixed in x and y, but I want the z direction to be mobile. However, this
setup is too constrained (bond atoms go missing) and thus, I want the boundary atoms to be flexible, yet
retain their average positions so that the long range elastic field (imposed by the boundary) is preserved.
If I don’t keep the boundary fixed, the edge partial might just slip out (in the absence of the field).

That is why I’m trying to keep the COM’s of the boundary region in place while allowing for some flexibility.

Regards,
Anirban

Dear Prof. Kohlmeyer,

I have a molecular crystal system with a single partial edge dislocation
in it.
I'm defining a cylindrical region around the dislocation line (which is
along z) which is flexible, while the region
outside the cylindrical region is FIXED. Now the z direction is periodic.

I'm running Umb. Samp. simulations while holding the dislocation core in
place using a bias potential,
and I'm holding the boundary region fixed in x and y, but I want the z
direction to be mobile. However, this
setup is too constrained (bond atoms go missing) and thus, I want the
boundary atoms to be flexible, yet
retain their average positions so that the long range elastic field
(imposed by the boundary) is preserved.
If I don't keep the boundary fixed, the edge partial might just slip out
(in the absence of the field).

That is why I'm trying to keep the COM's of the boundary region in place
while allowing for some flexibility.

​why don't you use fix spring/self then?
that seems a much simpler approach to me, yet should achieve pretty much
the same.

axel.​

Thank you for your comments.

Spring/self is an option I’ll try if this alternative fails. I just don’t want to affect the system energy by adding springs.

Currently, your suggestion on using chunk features seems most relevant.

I was wondering if I could try a more direct route using chunk:

(i) use chunk/atom molecule to get per-atom chunk id’s (which should be the same is mol-ids but that’s not important)
(ii) use com/chunk to get global array of all chunk com’s and store the initial state as well as current state
(iii) compute difference between initial and current states, convert values into an atom-style variable, and use fix move to adjust the coordinates of all atoms by the negative of their chunk/com displacements.

However, I’m not sure how to convert the global array values into atom style variables. It seems straightforward at first glance, but
i’m getting “illegal/invalid variable” errors:

compute cc1 all chunk/atom molecule
compute myChunk all com/chunk cc1

fix 1 all nvt temp {tempt} {tempt} 100.0

variable cid atom c_cc1
variable COMx atom c_myChunk[v_cid][1]

dump 1.3 mid custom 1 dump_midf.lammpstrj id mol type q xu yu zu v_cid v_COMx

So v_cid is storing the chunk id as I can check from the dump output, however I cant seem to use v_cid as the argument
of a global array c_myChunk.

Regards,
Anirban

Thank you for your comments.

Spring/self is an option I'll try if this alternative fails. I just don't
want to affect the system energy by adding springs.

​that is an argument that i don't agree with. since your molecules are not
freely rotating, you are essentially abusing the non-bonded interactions in
the same fashion. in fact, the more i think about what you describe as your
motivation, i believe that you are just trying to put a band aid over a not
very well designed simulation and thus the results will be tainted
regardless. i am not even certain​ that your interpretation of the
situation is correct.

Currently, your suggestion on using chunk features seems most relevant.

I was wondering if I could try a more direct route using chunk:

(i) use chunk/atom molecule to get per-atom chunk id's (which should be
the same is mol-ids but that's not important)
(ii) use com/chunk to get global array of all chunk com's and store the
initial state as well as current state
(iii) compute difference between initial and current states, convert
values into an atom-style variable, and use fix move to adjust the
coordinates of all atoms by the negative of their chunk/com displacements.

However, I'm not sure how to convert the global array values into atom
style variables. It seems straightforward at first glance, but
i'm getting "illegal/invalid variable" errors:

so that means in short: if your molecules move out of place, you just put
them back. sorry, but this sounds more like you are preparing for a career
at pixar than as a scientist. by pushing the molecules where they don't
want to be, you are adding energy to the system just as well.​ i most
certainly consider it much worse than using fix spring/self. in that case,
you have at least a measure of how much you are modifying your system
through the energy that is stored in all the springs. if that is
substantial, you need to serious review your overall simulation design.

​axel.​

But the molecules are freely rotating, their overall translation is what I intend to arrest. This is not very different from holding a set

of boundary atoms fixed in order to apply a constraint. Applying a spring/self to every atom will arrest both the translation and rotation

of molecules. Of course, if I could apply a fix spring/self to the molecule COM, that would be useful, but I don’t think that is allowed with the current commands.

Eitherway, I’ll definitely give fix spring/self a try. But I would still like to know if global array values can be converted to per-atom vectors.

Thanks,

Anirban

But the molecules are freely rotating, their overall translation is what I
intend to arrest. This is not very different from holding a set

what you are saying does not fit together: you talk about a molecular
crystal and free(!) rotation.
you are only increasing my suspicion that the main problem is rather in the
design of your simulation and that you trying to band-aid the real problem.

of boundary atoms fixed in order to apply a constraint. Applying a

spring/self to every atom will arrest both the translation and rotation
of molecules. Of course, if I could apply a fix spring/self to the
molecule COM, that would be useful, but I don't think that is allowed with
the current commands.

​you could apply spring/self to the atom closest to the COM.

Eitherway, I'll definitely give fix spring/self a try. But I would still
like to know if global array values can be converted to per-atom vectors.


you would have to create a massive formula that would match the global
array dimension with something that maps it to the atoms.
with 4000 molecules and 3 dimensions, that would be a sum over at least
1200 expressions. not something that you want to write by hand and
definitely something of a benchmark to LAMMPS scripting implementation.

anyway, it is impossible to discuss such matters in abstract terms and i
have said everything i want to say on this subject.

axel.

Dear Axel and list,

    I'd like to hold the COM's of about 4000 molecules stationary. It
    seems that

​what kind of molecules? and what is the physical justification for such a (rather unusual) request​?

The physical justification why I'm also interested in exactly this feature is that a
promising technique of coarse-graining not just static structure, but dynamics
is Mori-Zwanzig based projections. The theory can be found in C. Hijon et al.
"Mori-Zwanzig formalism as a practical computational tool".
Faraday discussions 144, 301-322 (2010) http://dx.doi.org/10.1039/B902479B

In order to measure the Langevin memory kernels required by this formalism,
MD simulations should be performed in the constrained ensemble where the
centres of mass of the CG units are kept fixed.

it might be possible to use the */chunk features in LAMMPS:
- group the atoms into chunks by molecule id.
- compute the force on the center of each chunk through fix ave/chunk
- convert this into a per atom force that can be added to each atom with fix addforce. this may require some LAMMPS script trickery, e.g. using a loop.
​- if there is no initial velocity, then only forces​ inducing rotation should remain.
- to protect against numerical noise, you may want to subtract the center of mass velocity from each chunk regularly, too.

of course, you could just write your own integrator that disregards per-chunk translation.

Thanks for the input! The last option seems the easier to implement.

Here is what I’d suggest, which has some overlap
with Axel and Anirban’s ideas.

(1) define chunks as molecule IDs
(2) define fix ave/chunk 1 1 1 chunkID fx fy fz
so that it calculates the total force per chunk (molecule)
every tilmestep
(3) add an option to fix addforce like “chunk chunkID fixID”
which adds a force to each atom based on looking
up its chunkID and indexing the per-chunk force in fixID,
actually you want the negative of that force

Only (3) requires some coding. You’d have to think
a bit about when in the tilmestep this happens, since (3)
normally happens in the middle of the tilmestep and (2)
at the end of the step. If that’s not workable, you could
just put the logic for (1), (2) in (3) (all the code added
to fix addforce), since it;s pretty simple to sum forces by
molecule into a 4000-length arrray, and MPI_Allreduce
it, which is what fix ave/chunk is doing anyway. You could
make this it’s own fix as well.

Steve

Dear Steve and list,

I was thinking of adding a capability to the variables module in lammps. Currently, as per the documentation:

Atom Values and Vectors (http://lammps.sandia.gov/doc/variable.html)
Atom values take an integer argument I from 1 to N, where I is the atom-ID, e.g. x[243], which means use the x coordinate of the atom with ID = 243. Or they can take a variable name, specified as v_name, where name is the name of the variable, like x[v_myIndex]. The variable can be of any style except atom or atom-file variables.

I was wondering what steps I need to take to add the capability so that global array arguments (defined within atom style expressions) be variables of atom/atom-file style as well. Currently, atom variable formulae can contain a large variety of expressions of atom variables, and adding array indexing capabilities using atom variables would be useful.

Sincerely,
Anirban