Query in fix reax/c/species command

Dear Lammps users,

I am trying to calculate molecular composition using fix reax/c/species.

I am simulating water system where H – atom type 1, O – atom type 2 and I am using the following fix command

fix 42 all reax/c/species 1 1 100 species.out cutoff 1 2 0.35

In the output file, I am getting C2H instead of H2O.

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.

I am using lammps Jan 13, 2014 version.

I am also attaching my input files herewith.

I will appreciate your suggestion/comments in this issue.

Thanks,

Sanjib

water-model-418-molecules.txt (93.2 KB)

ffield.reax (20.4 KB)

in.water-reaxff-01.txt (2.97 KB)

Sanjib,

[...]

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.

Ray

Hi Ray,

Thanks for your comments. Now it is working.

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.

Thanks,
Sanjib

Are you sure you are using VMD correctly?

What representation and what parameters do you use?

Axel.

You can try decreasing the BO cutoff to 0.2 or even 0.1.

Ray

In VMD, I tried with BOND, CPK,VDW drawing method.

Thanks,

sanjib

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.

Axel.

Dear Axel and Ray,

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.

I didn’t understand why this difference is occurring.

I have attached herewith the lammps input script.

I will appreciate your comments/suggestions in this issue.

Thank you.

Kind regards,

Sanjib

in.silica-silane-reaxff-01.txt (4.6 KB)

Dear Axel and Ray,

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.

axel.

Dear Axel,

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?

Thank you.
Kind regards,
Sanjib

Sanjib,

22.02.2014, 00:43, "Sanjib Chowdhury" <[email protected]>:

Dear Axel,

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.

Regards,
Oleg

Sanjib,

22.02.2014, 00:43, "Sanjib Chowdhury" <[email protected]>:

Dear Axel,

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.

almost 10 years ago i wrote and maintained a webpage about these kind
of things. the web page has since been taken down, but its content is
still available on my homepage (and many components have found their
way into VMD directly). check out.
https://docs.google.com/viewer?a=v&pid=sites&srcid=ZGVmYXVsdGRvbWFpbnxha29obG1leXxneDo0MzMwYzUwMmIyODU5NDFk
or
https://sites.google.com/site/akohlmey/redirect/cpmd-vmd.pdf?attredirects=0&d=1
if you'd rather download the PDF file.

axel.

Sanjib,

22.02.2014, 00:43, "Sanjib Chowdhury" <[email protected]>:

Dear Axel,

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.

Ray