Hi developers,
I am modifying fix_addforce.cpp to apply a force to a molecule even if only some atoms of the molecule are within the specified region where a force is applied. As I understand, it is possible to find out the molecule index to which an atom belongs to by using the array atom->molecule. My question is can the inverse be found, i.e., given a molecule index, can I find all the atoms belonging to the molecule?
Thank you,
Marcus
Hi developers,
I am modifying fix_addforce.cpp to apply a force to a molecule even if only
some atoms of the molecule are within the specified region where a force is
applied. As I understand, it is possible to find out the molecule index to
which an atom belongs to by using the array atom->molecule. My question is
can the inverse be found, i.e., given a molecule index, can I find all the
atoms belonging to the molecule?
you either have to build a corresponding mapping and
store it in your custom fix or do an (inefficient) loop over
all (local and ghost) atoms.
axel
Thanks Axel. Right now, I am just trying to do this for water molecules. If
the oxygen atom is located in the specified region, then a force is applied
to the whole water molecule. My plan is to loop over the local atoms, find
out whether an atom is an oxygen atom. If it is an oxygen atom and it lies
in the region, then apply a force to the whole molecule. I am having trouble
with finding a hydrogen atom bonded to the oxygen atom. The array I am
looking at is bond_atom. Does bond_atom[i][m] mean the global index of atoms
bonded to atom i? If so, how do I determine the local index? Thanks.
Marcus
You could probably do most of the work using the LAMMPS compute class
ComputePropertyMolecule....
The most efficient way to do this is how the TIP4P style works
where it requires the IDs of the 2 H atoms connected to the O
atom be consecutive IDs. E.g. O = ID 100, then H = 101, 102.
If this is the case then you can "lookup" the local ID of the
the 2 H atoms from the O atom ID, directly, using atom->map().
Steve