Dear LAMMPS developers,
I would like to ask if there is any compute command available that outputs the number of molecular collisions. I took a look at the compute coord/atom and contact/atom, however as the name suggests these commands count the number of atom-atom interactions. I was wondering if there is any option (by specifying a chunk id) to compute the chunk-chunk collisions based on the atomic interactions of each chunk.
Thank you very much,
Best regards,
Efstratios
Dear LAMMPS developers,
I would like to ask if there is any compute command available that outputs the number of molecular collisions.
not that I know of. I would not even know how to cleanly define what that would be in a general way.
I took a look at the compute coord/atom and contact/atom, however as the name suggests these commands count the number of atom-atom interactions. I was wondering if there is any option (by specifying a chunk id) to compute the chunk-chunk collisions based on the atomic interactions of each chunk.
computer logic does not work like human logic and thus you cannot just compute something and then throw it at the computer and tell it to “figure it out”.
this seems like a job that is best done in post-processing, possibly based on storing per-atom contact information in the trajectory for use during the post-processing. you may want to check out using compute pair/local as well.
it all depends on what you actually want to learn and what information to extract.
What is a molecular collision?
Steve
Dear Axel,
Dear Steve,
I apologize for the misunderstanding. What I would like to measure in some ReaxFF combustion simulations would be the number of collisions between molecules. To further clarify what I mean, consider an H2 and O2 domain. If between two O2 molecules an O–O distance was found less than a certain value (assuming for example the Van der Waals radius that value would be ~3 Ang.) then the O2 molecules are considered a collision pair. Similarly, between an H2 and O2 molecule, if an H atom was found at a certain distance from an O atom (assuming for example the Van der Waals radius that value would be ~2.6 Ang.) then the O2–H2 molecules collided.
It would be very helpful to have a per-atom compute that considering pairs of neighbouring colliding molecules, counts one collision for all atoms of both molecules. Therefore, a collision between atom-atom of different molecules would be broadcasted in all atoms of those involved molecules, even if the atoms are outside this cutoff distance. Another possible characteristic of that command would be: when considering the next atom of the same molecule, which is within that cutoff distance of “interaction”, a collision should not be counted twice with the same colliding molecule.
Hope I was able to make myself clearer now.
Best regards,
Efstratios
there are multiple problems here:
all of that points to a) this being very specific for the problem you are looking at and b) this being best done in post processing as you need to be tracking pairs of atoms for multiple timesteps to determine atomic collisions and then apply the bond order information as to which groups of atoms are currently molecules. implementing this in LAMMPS with proper parallelization would be a massive undertaking for a rather peculiar application.
Axel.
Two ideas that might help. The dump local pair command uses the
compute pair/local and compute property/local commands for
output. This means you can write a dump file for I,J pairs of
atoms and their distance at any timestep. You could also use
the dump_modify thresh command to limit the output to certain
types of atoms. If you post-process that kind of output you could
find nearby pairs of atoms.
It would be nice if the thresh option allowed you to select on pair
distances, but I can’t think of an easy way to do that.
You can also use the output of compute pair/local (the dist column) as input to
the compute reduce command or fix ave/histo to generate some stats
on the set of distances it computes.
Steve
Dear Steve and Axel,
Thank you for your suggestions and time to help me. After our discussion I tried to write a new compute command that computes the collisions between molecules. I believe that post processing that information would be more cumbersome, so I tried to go with that compute command. Concerning the validity of that “method” of identifying colliding molecular pairs, I agree with Axel that you cannot be certain a collision occurs with those standards, but it still gives you an estimate of what is happening in your simulation. Moreover, I have seen that method being used in a few papers (such as 10.1021/acs.jpca.8b11686). I agree that due to the high KE of atoms in a combustion simulation you can be confident that such close interaction is happening due to a collision.
I attach the compute command in case you would like to consider integrating it; or if you find it redundant someone else could benefit. There is a detail syntax of this compute command in the header file. I should note that in the .cpp file there is a draft solution for the chunk recognition with ReaxFF that I have made, because the reax/c/species per-atom array (which contains the molecule ID of each atom) erroneously outputs the information from the previous timestep. That eventually results in wrong data if multiprocessors change the local atomic ids. I have addressed that issue in a previous question. Nevertheless, that part of the .cpp file can be easily substituted with the vector of the compute chunk/atom command. I have quickly tested the collision/contact command with multiprocessing as well. There is an issue when a molecule is in-between processors, and a collision takes place near the communication boarder. In that case a collision may be counted twice. I am sure this command could be further improved with your level of computer programming and knowledge, but I found it a good solution for my problem.
As I mentioned in my previous email, this command counts molecular collisions based on intramolecular atomic distances. Atomic radii are inputted by the user for each atom type. All atoms of the same molecule acquire the information that a collision occurred. A collision is not counted twice even if two atoms of the same molecule are close enough with that neighbouring molecule. Collision events are summed up per molecule. Therefore, if a C12H26 molecule collides with 3 O2 molecules all atoms of the C12H26 molecule will have a count of collisions equal to 3, and all 3 O2 molecules will have a count of 1. The output is a per-atom vector although it could be probably more appropriate to output a per-chunk information. Nevertheless, I thought it is redundant to memory destroy and create arrays based on when the number of molecules change (with ReaxFF this should happen quite often). I have not tested the efficiency of the parallelization, but it does follow the LAMMPS standards. If I run into bugs, I could inform you further if this is something of interest.
Thank you again for your help,
Best regards,
Efstratios
compute_collision_contact.cpp (8.38 KB)
compute_collision_contact.h (3.72 KB)