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?
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.
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
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).
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.