# Direct calculation of COM RDF

Hello dear lammps users,

I’d like to calculate the radial distribution function in my simulation of organic solid.

I know the compute rdf command can compute rdf between paris of atoms, but I’m looking for a simple way to calculate RDF for center of mass of molecules.

I use the ReaxFF pairstyle, so I don’t have molecules defined in the strict sense, but just groups of atoms.

I also know there is a possibility to calculate RDF within VMD, but then again, it only lets me calculate RDF between two selections (atom types).

Is there a rather simple way to calculate RDF between molecules (center of mass of molecules)?

Thanks

David Furman, MSc. student
The Fritz Haber Research Center for Molecular Dynamics
Institute of Chemistry
The Hebrew University of Jerusalem
Jerusalem, ISRAEL

Phone: +972-8656-8909

Hello dear lammps users,

I'd like to calculate the radial distribution function in my simulation of
organic solid.

I know the compute rdf command can compute rdf between paris of atoms, but
I'm looking for a simple way to calculate RDF for center of mass of
molecules.

I use the ReaxFF pairstyle, so I don't have molecules defined in the strict
sense, but just groups of atoms.

so your _first_ problem is to program a feature in LAMMPS
that teaches it how to determine what _is_ a molecule.

if possible, i would just pick an atom that would be close
to the COM. that would usually be good enough.

I also know there is a possibility to calculate RDF within VMD, but then
again, it only lets me calculate RDF between two selections (atom types).

again, the same problem. you first have to indicate what is a molecule.
VMD, like LAMMPS, assumes properties to be static unless there is
a handler/script/fix that manipulates it.

Is there a rather simple way to calculate RDF between molecules (center of
mass of molecules)?

i would probably want to write a postprocessing tool
that converts your atom-based trajectory into one that
lists only the centers of masses (as atoms) and then
compute the RDF (either through my own code, or using VMD).

cheers,
axel.

Thanks Axel !

I’ve written a simple tcl script to calculate the center of mass coordinates of each molecule in a custom LAMMPS dump file. Then used the g® tool in VMD to get the RDF.

Another question I have regarding VMD is each time I load a dump trajectory file from lammps into VMD, I see that VMD connects atoms across the two sides of the simulation box. So bonds between edge atoms appear really long (all over the simulation box)… How can I set it to show only bonds between atoms inside the box, and not connect atoms if they are on the edges ? Is this a LAMMPS or a VMD issue ? I hope this description is clear enough…

david,

Thanks Axel !

I've written a simple tcl script to calculate the center of mass coordinates
of each molecule in a custom LAMMPS dump file. Then used the g(R) tool in
VMD to get the RDF.

excellent!

Another question I have regarding VMD is each time I load a dump trajectory
file from lammps into VMD, I see that VMD connects atoms across the two
sides of the simulation box. So bonds between edge atoms appear really long
(all over the simulation box)... How can I set it to show only bonds between
atoms inside the box, and not connect atoms if they are on the edges ? Is
this a LAMMPS or a VMD issue ? I hope this description is clear enough...

yes. it is mostly a VMD issue. VMD computes the bond information from
heuristics (unless you first read a data file through topotools) and maintains
this bonding pattern over the trajectory. if your trajectory uses "wrapped"
coordinates (the default), then atoms will be replaced by periodic images
without respect for the bonds. potential remedies are:

a) don't visualize bonds
b) use the dynamic bonds representation
c) use "dump_modify <dumpid> unwrap yes" in LAMMPS
and then use the "wrap" function in the pbctools plugin
in VMD to wrap coordinates back into the principal unit
cell, but use the option '-component fragment' to have
entire fragments (= molecules) wrapped
d) program an option into LAMMPS' dumps to unwrap only
bonds that stretch across periodic boundaries

cheers,
axel.

p.s.: FYI the VMD mailing list is at [email protected]... and
cross-posting between mostly unrelated mailing lists is
usually not well received unless most subscribers and
particularly those that would answer are subscribed to
both lists.