Dear LAMMPs users,
I’m writing a custom fix and need to know the number of atoms of each molecule I define in the data file. I’ve searched through the documentation and source files and I found 3 arrays (molecule, molatom and molindex).
I am aware that the first one holds the molecule ID for each atom but I couldn’t find what the other 2 are, I wonder if they have the information I am looking for, but when I try to access their elements I get segmentation faults, and when I access them as if they were variables ( printf(%d, molatom), don’t ask me why I tried this) I get = 0 for both of them.
I was hoping this information was already available in some pointers. If not I could loop over all atoms and get it myself, if this is the case I wonder where I should do this: in the constructor, in init() or setup()?
Thank you very much,
The identity of molecules is completely irrelevant to the kind of force fields implemented in LAMMPS.
All it cares about are atoms and - if provided - bonds, angles, dihedrals, and impropers.
Thus the molecule ID is completely arbitrary and can be used for all kinds of different purposes very few of them related to the interactions (some INTERLAYER potentials use molecule IDs to identify different layers, rigid fixes can use it to identify rigid bodies which may or may not coincide with molecules).
Because of that, LAMMPS does not keep track of molecules. If you can depend on you molecule IDs to be consistent, you can use compute chunk/atom to identify molecules by molecule IDs as “chunks”.
When writing a custom fix you can create the compute as part of the constructor (or require an argument with a suitable compute ID), then you can trigger the compute during execution and process the information from it for your needs.