Binning by molecule COM?

Hi All,

I have recently been working on a shock code. One of the things I would like to compute is the intermolecular temperature (based on the COM KE of molecules) as well as the intramolecular temperatures (based of net KE of the constituent atoms inside of a molecule).While these can be done in post-processing using the molecular COM coordinates and VCM, I am trying to figure out how to do this on the fly in LAMMPS. Essentially, I would like to do a 1D spatial chunk of the molecules in a system, not the individual atoms. My current approach is to use a series of compute reduce/atom/chunk commands to spread the per molecule quantities to all of the atoms comprising it, but this causes problems when computing the molecular COM KE of the molecules (since different parts of the same molecule, which all have the same COM KE, are in different spatial bins, which have different KE). Perhaps some code will make the error clearer:

compute		molec Rest chunk/atom molecule ids once nchunk once compress yes
# chunking atoms by molecule

compute 	thermT Rest temp/chunk molec temp com yes adof 2.857
# computing intramolecular temperature
compute		comke Rest temp/chunk molec kecom
# computing COM KE, which can be converted into intermolecular temperature

compute		ldz Rest chunk/atom bin/1d z 20 3 ids once nchunk once compress yes units box
compute		edz Rest chunk/atom bin/1d z 20 6 nchunk once compress yes units box
#Making two bins of atoms, one eulerian (edz, for density) and one lagrangian (ldz, for all other props)

compute		bcomke Rest temp/chunk ldz kecom
compute		bcomkes Rest chunk/spread/atom ldz c_bcomke[1]
compute		chsize Rest property/chunk ldz count
compute		schsize Rest chunk/spread/atom ldz c_chsize

compute		ctemp Rest temp/chunk ldz kecom com yes

compute		comT Rest chunk/spread/atom molec c_comke[1]
compute		intT Rest chunk/spread/atom molec c_thermT[1]
compute		molsze Rest property/chunk molec count
compute		smolsze Rest chunk/spread/atom molec c_molsze
# these commands assign the previous computes to all of the atoms in the molecule, so we can now use compute ave/chunk on them
# bin COM KE must be subtracted out to get good temperature for the comT

variable	comTvar atom (2*4184*(c_comT*(c_schsize/c_smolsze)-c_bcomkes))/(6.022e23*3*1.38e-23*(c_schsize/c_smolsze))
#converting COMKE to COM Temp, multiplies comKE of molecule by # of molecules in the bin to get a total comKE of all 
# the molecules in the bin, then subtracts COM KE of the bin itself and divides by number of molecules to get COMKE per molecule.

fix comtempprof Rest ave/chunk 10 5 500 ldz v_comTvar file comtempprofile.data

# The problem here is that different atoms in the same molecule will have different bin KE's subtract, resulting in bad results at bin boundaries
# This could be fixed if the bins were done based on molecular COM coordinates. Is it possible to bin things this way?

Please let me know if this is unclear.

I don’t think this is currently possible in LAMMPS. If a molecule straddles a processor boundary it would also require communication to get it right. And a long-chain alkane could span multiple processor domains.

Instead of the center of mass I would just pick the atom that is (typically) closest to the center of mass. This as arbitrary a choice as would be the center of mass for linear molecules and usually a very good approximation for small or globular molecules. If needed you can assign this atom a new atom type so it is easier to identify and put into a group.

Since you are collecting per molecule data, you must not overinterpret the distribution data anyway. If you bin all atoms data is too spread out, if you bin only one point it is not spread out enough.

1 Like

What about:

compute   mols_compos       Rest com/chunk molec  # calculate COMs from molec chunks
compute   atom_compos       Rest chunk/spread/atom  molec c_mols_compos[*]  # now each atom knows its molecular COM
variable  atom_composz_bin  atom floor(50*(c_atom_compos[2]-zlo)/lz)  # 50 bins for molecular COM z coordinate
compute   bin_by_composz    Rest chunk/atom v_atom_composz_bin
1 Like

This is a brilliant suggestion. I had no idea that chunk/atom could be used with a variable, so this opens up a lot of doors for future LAMMPS use. I am curious about how exactly the variable command works here. Suppose we are considering some atom with a COM_z coordinate of 25. If we make this command congruent with the other bins, (i.e. zlo = 20, lz = 6), it seems that this variable would be something around 0. So this is what assigns the bin number to each atom. My question is where does the integer multiple in front come in? It seems like this would be 95 for atoms with COM_z 20-26, 95*2 for 26-32, etc. I must admit this is a little bit confusing and I’m not sure how this would translate into a set of binned atoms.

Perhaps this is a lack of understanding of how the chunk/atom command parses per-atom variables to chunks. Would just like to understand this seemingly powerful technique. A quick plot of the variable function shows a straight line. I guess I can’t really see how the floor function would either give the number of chunks or assign atoms to these chunks

Hi! It pays to check documentation – I did not know how to use chunk/atom with a variable until half an hour ago either.

zlo and lz are so-called thermo variables or thermodynamic keywords (see thermo_style documentation for details under “custom keywords”) which store the low z-boundary and the z-length respectively of the simulation box. I don’t think chunk/atom takes similar keywords (which seems to be what you’d mean with something like “zlo = 20, lz = 6”) and even if it does, this is a completely different concept.

What I’m doing is calculating the bin number by brute force. Something like scale_z = (v_z-zlo)/lz returns a scaled z from 0 (v_z = zlo) to 1 (v_z = zhi). Calling bin_z = ceil(50*scale_z) then returns a bin number from 1 to 50 (LAMMPS reserves bin 0 for special cases so you should use ceil, not floor from my earlier entry).

EDIT: As per documentation, compute chunk/atom does the rounding for you, so you don’t need to use floor or ceil (you still have to check for off-by-one errors carefully).

1 Like

Excellent. I suppose I should have read more closely. Thanks for your insight!

By the way – my colleagues and I are working separately on building more integrated per-molecule capabilities into LAMMPS. It would be quite useful to us to know your use case (we can’t promise we’d prioritize it, but it would be nice if we were helping the wider community at the same time as ourselves). Feel free to PM me to get in touch.