Hi all, I am interested in computing molecular properties for molecules which use bonds. In looking at atom.cpp::data_bonds and neigh_bond.cpp::bond_all, it looks like all bonds must be between real atoms on the same processor, i.e., not processor or periodic ghosts. Is this correct? Otherwise, what happens when map(bond_atom[i][j]) < 0? While I realize molecules are not required to be collocated on a processor, it would help a great deal to understand how bonds (and angles and dihedrals) work in this mode.
Thanks
Jeremy
hi jeremy,
Hi all, I am interested in computing molecular properties for molecules which use bonds. In looking at atom.cpp::data_bonds and neigh_bond.cpp::bond_all, it looks like all bonds must be between real atoms on the same processor, i.e., not processor or periodic ghosts. Is this correct? Otherwise, what happens when map(bond_atom[i][j]) < 0? While I realize molecules are not required to be collocated on a processor, it would help a great deal to understand how bonds (and angles and dihedrals) work in this mode.
as far as i understand the atom mapping,
you will get map(i) < 0 only for atoms
that are neither local or ghost atoms.
thus an atom in a bond has to be either
a local or a ghost atom.
how the bond list is compiled now depends
on whether newton_bond is active or not.
in the first case, it is made sure that forces
from a single bond are only calculated once,
the forces affecting ghost atoms are communicated
after all forces are computed.
in the second case, forces between local atoms
are computed as before, but for bonds between
local and ghost atoms, only the local contribution
is kept. this procedure is extended in a similar
fashion for angles, dihedrals and impropers.
as far as i understand the code, the molecule tag
is available only to ease bookkeeping but has no
direct relation to bonds and so on.
HTH,
axel.
Hi Axel, thanks for the information. I have one more question: is the bond list only symmetric if newton_bond is off? If so, can the fix adjust this parameter like it can the regular neighbor list?
Thanks,
Jeremy
Hi Axel, thanks for the information. I have one more question: is the bond list only symmetric if newton_bond is off? If so, can the fix adjust this parameter like it can the regular neighbor list?
No. There is no infrastructure to have multiple bondlist requests.
Neighbor lists are special in this respect. But you can test whether
newton_bond is active or not and refuse to run if it is on. Note that
can be set independently from the setting for nonbonded interactions.
Axel.
All of Axel's feedback is good. I would also suggest
you look at any of the compute foo/molecule code
which already has the logic to do things on a per-molecule
basis.
Steve
Thanks Axel and Steve. I was trying to be slick and make some approximations for small molecules, i.e., of small enough diameter that all atoms in a molecule will be either real or ghosts on a processor. This way we can avoid a lot of the MPI communication and reduce the storage overhead because we are doing these computations possibly every timestep. Most of the molecules we are looking at now are small, e.g., water. But for the general case, the approach in the computes is the only way I can think of.
Jeremy
Hi Steve and Axel, I've made some progress on this and now I have questions regarding periodic boundaries. I saw how in the compute com molecule routine, the atoms' coordinates are corrected for periodicity. I noticed they are not corrected for a triclinic domain. Is this intentional, and if so, why? Also, how does this avoid having molecules with atoms crossing periodic boundaries have a consist center of mass because one atom could be off by the periodic distance? For large molecules in particular, it seems like half the atoms could be on one side, half on the other, so the center of mass would be exactly in the center of the domain which would be incorrect. Finally, after I compute the molecule's center of mass, I want to make sure that location is within the computational box. How do I transform a specific coordinate without using the image variable which is done per atom?
Thanks (yet again),
Jeremy
Hi guys, so after a bit of digging, I think I found most of the answers in fix_rigid. Basically what we have to do is set the image variable of all the atoms in a molecule to be consistent with how it crosses a periodic boundary, but only if we start them crossing periodic boundaries. I've also got some code to do the periodic mapping, so my only remaining question is why the computes don't consider the triclinic mapping (but please let me know if what I've said above is incorrect).
thanks,
Jeremy
There are probably several places in the code where triclinic effects
are not accounted for, due
to the triiclinic ideas being added later. Your discussion of
periodic boundaries seems
incorrect, at least for orthogonal systems. If you have given your
initial molecule correct
image flags on each atom, then it will be unwrapped correctly wrt the
periodic box,
and the correct COM calculated. That COM is not mapped back into the box,
b/c that isn't the way it would typically be useful. Instead you want
to know how
fhat molecule has drifted over time, which is the unwrapped COM.
There is a function domain::remap(x,image) or domain::remap(x)
that will remap a x[3] from outside the box back into it, but its not used
by the COM computes.
Steve