It seems the fix reax/c/species command taking the atom type information
from the ffield.reax file where
atom serial is C, H, O,N. Atom type 1 & 2 in the fix command are pointing C
and H respectively in the ffield.reax file not H & O that I am assigning.
No, it does not take anything from the force field file. Default
elements are in the sequence of C, H, O, N, but you can change that
using the element keyword. Please see the doc page for more details.
I have another issue about reax/c/species.
There are 418 H2O molecules in the water model that I am relaxing at 300 K
using NVT.
VMD animation shows, no H2O molecules dissociation.
But reax/c/species output files shows variation of H2O molecules with time
(MD step). This is about 400-418.
How can I make this range more narrow?
For H-O bond, bond order cutoff is about 0.33.
I have tried with different cutoff values around 0.33. But couldn't overcome
the problem.
In VMD, I tried with BOND, CPK,VDW drawing method.
None of those will show bond breaking. VMD knows nothing about what model you ran and assumes that you do classical MD with a fixed bond topology. Only the dynamic bonds representation can show bond breaking, and it only uses a simple.distance criterion.
I am still facing problem with 'reax/c/species' and corresponding VMD
visualization.
Their outputs are not in concert.
I am conducting relaxation of 64 GPS (C9H20O5Si) molecules at NVT for 100
ps. Following commands are used for
reax/c/species and VMD dumping.
fix 410 all reax/c/species 1 1 100 species-water.out element C H O Si
dump 2 all xyz 5000 dump.silica-water-coord-xyz-format
dump_modify 2 element C H O Si
'reax/c/species' shows there are 64 GPS molecules after the relaxation while
VMD (with DynamicBonds method) shows molecule dissociation.
as i have already explained, VMD does not now anything about physics
and chemistry (and very deliberately so) and thus you only get to see
what you tell VMD to show you using information you provide. please
have a closer look at the VMD user's guide. with dynamic bonds the
sole parameter to decide over bonding is a single distance cutoff.
that is surely unsuitable for your system. in fact, the question of
what is a bond and what not, is not that easy to answer, since the
physics doesn't care, it just has interactions of (continuously)
varying strength between atoms. the concept of bonds is an
interpretation and simplification based on observations.
with that in mind, any expectation that a simple visualization tool is
able to correctly classify bonded and non-bonded interactions is
unfounded.
Thanks for your comments.
Yes, as you said changing the distance cutoff value in VMD, bond
connectivity network can be changed.
Also changing the bond order value in LAMMPS (though 0.3 is standard),
different sets of species can be obtained.
So it seems there is no straightforward answer about the bond formation and
exact quantification of the types and number of molecules (reactants,
products) while studying reaction using reactive potential. Right?
I need another suggestion. Using the ' reax/c/species', I can get
information about different species (i.e., molecules type, number and COM).
I need to know the xyz coordinate or atom ID in a particular species. How
can I get these in LAMMPS domain?
Thanks for your comments.
Yes, as you said changing the distance cutoff value in VMD, bond
connectivity network can be changed.
Also changing the bond order value in LAMMPS (though 0.3 is standard),
different sets of species can be obtained.
So it seems there is no straightforward answer about the bond formation and
exact quantification of the types and number of molecules (reactants,
products) while studying reaction using reactive potential. Right?
Right. It is also not always straightforward to define the concept of "molecule" in a reactive system.
Keep in mind also, that as Axel already said, you are unlikely to obtain the same results with reax/c/species and visualization in VMD. And you definitely shouldn't try to achieve that by tuning reax/c/species threshold bond orders. ReaxFF gives much more consistent bonds than simple distance criterion used in VMD.
I need another suggestion. Using the ' reax/c/species', I can get
information about different species (i.e., molecules type, number and COM).
I need to know the xyz coordinate or atom ID in a particular species. How
can I get these in LAMMPS domain?
Thanks for your comments.
Yes, as you said changing the distance cutoff value in VMD, bond
connectivity network can be changed.
Also changing the bond order value in LAMMPS (though 0.3 is standard),
different sets of species can be obtained.
So it seems there is no straightforward answer about the bond formation and
exact quantification of the types and number of molecules (reactants,
products) while studying reaction using reactive potential. Right?
Right. It is also not always straightforward to define the concept of "molecule" in a reactive system.
Keep in mind also, that as Axel already said, you are unlikely to obtain the same results with reax/c/species and visualization in VMD. And you definitely shouldn't try to achieve that by tuning reax/c/species threshold bond orders. ReaxFF gives much more consistent bonds than simple distance criterion used in VMD.
I need another suggestion. Using the ' reax/c/species', I can get
information about different species (i.e., molecules type, number and COM).
I need to know the xyz coordinate or atom ID in a particular species. How
can I get these in LAMMPS domain?
I think you'll need to write some code.
...and even better would be not just "some code", but rather a plugin
for VMD that would read and store the output for reax/c/species in
some internal data structures and then dynamically apply it to a
trajectory and update bonds and visualization parameters dynamically.
Thanks for your comments.
Yes, as you said changing the distance cutoff value in VMD, bond
connectivity network can be changed.
Also changing the bond order value in LAMMPS (though 0.3 is standard),
different sets of species can be obtained.
So it seems there is no straightforward answer about the bond formation and
exact quantification of the types and number of molecules (reactants,
products) while studying reaction using reactive potential. Right?
Right. It is also not always straightforward to define the concept of "molecule" in a reactive system.
I will just add that there is no true bond or molecule concept in any
system described by a pair_style. Since pair_style is an atomic style
description and atoms interacts with whoever is in the cutoff -
knowing nothing about bonds or molecules. This is true for both
reactive and non-reactive pair styles. So you can only "approximate"
a bond using an atomic style pair_style, and there are many way to do
that approximation.