PDMS theramal decomposition

Hello!
I would like to simulate thermal decomposition in LAMMPS but this software is new to me. I am thinking about making a simulation that would involve several chains of PDMS in the environment of oxygen and water at high temperatures (about 380K). Additionally I would like to compress that system with a piston (with the pressure of, in example 20MPa). Now i would like to put words into action but I have several questions:

  1. I have positions of PDMS chain in .xyz file and I am wondering if is it possible to attach somehow this file in my input script.
  2. Is it possible to replicate PDMS chain (for example ten times) in the simulation box?
  3. Is it possible to simulate pressure increase with a moving piston?
    I would be very greatful for all hints.
    Best regards
    Andrzej Battek

Hello!
I would like to simulate thermal decomposition in LAMMPS but this software

hmmm... thermal decomposition hints at bond breaking.
LAMMPS is a classical MD code, where bonds are usually
considered unbreakable. you have to either use a monte carlo
like approach (see fix bond/break and fix bond/create), but
those usually considered for bead-spring model polymers.
or you need to use a bond-order potential that computes
bonded interactions on the fly (like REAXFF), for which you
would need a force field parameterization specific enough
to model your environment well. this is something that you
should check carefully, before wondering about technicalities.

is new to me. I am thinking about making a simulation that would involve
several chains of PDMS in the environment of oxygen and water at high

what is a PDMS?

temperatures (about 380K). Additionally I would like to compress that system
with a piston (with the pressure of, in example 20MPa). Now i would like to
put words into action but I have several questions:
1. I have positions of PDMS chain in .xyz file and I am wondering if is it
possible to attach somehow this file in my input script.

no. not really. best you create a data file. there are
multiple approaches that can help there. in any case
you will have to assign atom types somehow and tell
LAMMPS how to identify bonds and so on (if needed).
a .xyz file does not contain that information.

2. Is it possible to replicate PDMS chain (for example ten times) in the
simulation box?

not as such, but you can replicate the entire simulation
box with all that is in it and that is for all practical purposes
the same thing. also some of the preprocessing tools
that can output data files support replication of molecules.

3. Is it possible to simulate pressure increase with a moving piston?

yes. you don't need a piston, just do a variable cell calculation
and crank up the target pressure as needed, if you want to
set the pressure. or us fix deform to compress your simulation
box as desired, if you want to pre-set the volume rather than
the pressure.

I would be very greatful for all hints.

the best hint i can give is to study the provided examples
and spend some time just reading through the manual and
getting familiar with all the various features and potentials.
making good choices depends on being well informed. :wink:

cheers,
     axel.

I'll just add that there is a fix wall/piston command
that may be of use.

Steve

I’ll just add that there is a fix wall/piston command
that may be of use.

Steve

Thanks Steve. I guess this command will be very useful in a while.

well, there is more than just this one line format.
the read_data page lists the entire format. it is also
dependent on the atom style you choose.

That is correct Axel but I have read that when using ReaxFF atom_style should be ‘charge’.
I have the position of atoms in my chain but I lack the values of charges. How can I calculate this values and using what software? I have never wrote that kind of script so I don’t know what to do. Maybe there are some tutorials that show how to do that?

That is correct Axel but I have read that when using ReaxFF atom_style
should be 'charge'.

topotools can output data files correctly formatted for atom_style charge.

I have the position of atoms in my chain but I lack the values of charges.

just use some reasonable approximate initial guess...

in reax the charges will be recomputed by a charge equilibration
procedure, so you don't have to have them perfectly from the onset.
you may even get away with keeping them at 0.0 initially, but it is
probably better to at least break the symmetry and initialize
the trends (positive/negative) as expected.

How can I calculate this values and using what software? I have never wrote
that kind of script so I don't know what to do. Maybe there are some
tutorials that show how to do that?

topotools has some tutorial fragments. the rest is basic VMD
and Tcl scripting knowledge that can be obtained from the
VMD user's guide and Tcl tutorials.
https://sites.google.com/site/akohlmey/software/topotools

axel

How can I calculate this values and using what software? I have never wrote
that kind of script so I don’t know what to do. Maybe there are some
tutorials that show how to do that?

topotools has some tutorial fragments. the rest is basic VMD
and Tcl scripting knowledge that can be obtained from the
VMD user’s guide and Tcl tutorials.
https://sites.google.com/site/akohlmey/software/topotools

Dear Axel,
I have downloaded Topotools from your website (topotools-1.3-20120215.tar.gz) but I don’t know how to make it work. Is there any step-by-step instruction describing what to do after unpacking .tar.gz file? I have gone through tutorial and it seems that everyting what I need to do to convert my PDMS.xyz file to .data is there.
Best regards
Andrzej Battek

Dear Axel,
I have downloaded Topotools from your website
(topotools-1.3-20120215.tar.gz) but I don't know how to make it work. Is

topotools is also part of VMD and requires VMD. so all you have to do
is to find the location of the bundled topotools package and add the
updated plugin or use a recent alpha test version of VMD, where the
update is already applied.

there any step-by-step instruction describing what to do after unpacking
.tar.gz file?

there is only one step. move the directory with the updated plugin
into your VMD tree. see INSTALLATION. i just updated it for
topotools 1.3 and VMD 1.9.1.

axel.

Hello,

I have obtained “.data” (data2.pdms) file of my PDMS chain (PDMS.xyz) using Axel Kohlmeyer’s Topotools. I have set atom style to charge because I would like to simulate thermal decomposition of PDMS using ReaxFF reactive force field. I have been wrtiting that kind of script for the first time, so I would really appreciate if someone check my work. Also I have some questions:
I would like to put some H20 and O2 particles in simulation box and watch what will happen. Should I add them to my Lammps input file? Or maybe I can do it Lammps script?
Now i don’t have any bonds in my chain? Will ReaxFF make bonds by computing them from interatomic distance? Does this force field compute distance between every pair of atoms?
Best regards
Andrzej Battek

PDMS.xyz (3.19 KB)

data2.pdms (5.26 KB)

pdmsxyztodata1.tcl (1.64 KB)

It is probably easier to create a data file containing PDMS, H2O and
O2 for the input script to read.

No, you do not need bonds for ReaxFF simulations. ReaxFF will
determine them for you based on bond orders, which depends on bond
distances, angles, dihedrals and charges.

The tricky part is, however, do you have a good enough ReaxFF
parameterization for PDMS and more importantly one that describes PDMS
interacting with H2O and O2?

Ray

Hello,

I have obtained ".data" (data2.pdms) file of my PDMS chain (PDMS.xyz) using
Axel Kohlmeyer's Topotools. I have set atom style to charge because I would
like to simulate thermal decomposition of PDMS using ReaxFF reactive force
field. I have been wrtiting that kind of script for the first time, so I
would really appreciate if someone check my work. Also I have some
questions:
I would like to put some H20 and O2 particles in simulation box and watch
what will happen. Should I add them to my Lammps input file? Or maybe I can
do it Lammps script?

it would be difficult to add them after you have read in the data file,
but combining multiple molecules and then creating a data file or
merging multiple data files is fairly straightforward to do with VMD
and topotools scripting.

cheers,
     axel.

Dear Ray,

2012/8/25 Ray Shan <[email protected]…24…>

It is probably easier to create a data file containing PDMS, H2O and
O2 for the input script to read.

No, you do not need bonds for ReaxFF simulations. ReaxFF will
determine them for you based on bond orders, which depends on bond
distances, angles, dihedrals and charges.

The tricky part is, however, do you have a good enough ReaxFF
parameterization for PDMS

I don’t know exacly what do you mean. How/Where can I find if my parametrization for PDMS is good enough? I have a force field file that describes PDMS, maybe that is what you are talking about? (I have attached it to this mail)

and more importantly one that describes PDMS
interacting with H2O and O2?

I have thought, if all the atoms involved are described by the force field, that means they will interact with each other. Correct me if I am wrong.

I have also a question to all users. I have obtained PDMS “.xyz” file from my professor, so I didn’t have to model the chain. But now I want to add the H2O and O2 molecules to my system but I lack the “.xyz” files of these molecules. Could you advise me which software I should use to create that kind of files? Will that software do a energy minimization, so I will obtain an accurate model of, for example, water molecule (with an accurate bond distance and angle between bonds)?

Best regards
Andrzej Battek

I have also a question to all users. I have obtained PDMS ".xyz" file from
my professor, so I didn't have to model the chain. But now I want to add the
H2O and O2 molecules to my system but I lack the ".xyz" files of these
molecules. Could you advise me which software I should use to create that
kind of files? Will that software do a energy minimization, so I will obtain

come on. you don't know the structure of O2 and H2O?

an accurate model of, for example, water molecule (with an accurate bond
distance and angle between bonds)?

what does it matter? you are doing a simulation with your
force field, right? so all you need is some ballpark numbers.
once you run with the force field, the bond length will be
determined by the force field in any case. just start with
a little energy minimization.

here is a short suggestion for a procedure:
- write an .xyz file for a single O2 and H2O each (or pdb)
- use the packmol program to add a number of each of
  those to your system (according to your liking and
  what is reasonable).
- make a data file out of the remaining compound file.
- minimize/equilibrate in LAMMPS and rejoice.

+ alternatively, you can just write a script that generates
   a mixture of several O2 and H2O molecules and
  combine that with the existing data file and remove
  overlaps (see the corresponding tutorial section on
  the topotools homepage).

axel.

Hi Andrzej,

Replies below, thanks.

Ray

The tricky part is, however, do you have a good enough ReaxFF
parameterization for PDMS

I don't know exacly what do you mean. How/Where can I find if my
parametrization for PDMS is good enough? I have a force field file that
describes PDMS, maybe that is what you are talking about? (I have attached
it to this mail)

I don't see the attachment. If you have a reactive force field for
PDMS that has been published, it is probably good enough. However, to
make sure yourself, you should try to reproduce some of the properties
of PDMS that you consider important, such as bond lengths, angles,
stress-strain relations of bonds, heat of formation and/or bond
dissociation barriers.

and more importantly one that describes PDMS
interacting with H2O and O2?

I have thought, if all the atoms involved are described by the force field,
that means they will interact with each other. Correct me if I am wrong.

A force field that works with PDMS + H2O + O2 is one thing, one that
reasonably describes their interactions so that important physical and
chemical properties are captured is another; what I meant was the
latter.

2012/8/28 Ray Shan <[email protected]…1125…>

Hi Andrzej,

Replies below, thanks.

Ray

The tricky part is, however, do you have a good enough ReaxFF
parameterization for PDMS

I don’t know exacly what do you mean. How/Where can I find if my
parametrization for PDMS is good enough? I have a force field file that
describes PDMS, maybe that is what you are talking about? (I have attached
it to this mail)

I don’t see the attachment.

Sorry, I’ve forgotten to attach the file. Here it is.
Best regards
Andrzej Battek

ffield (21.2 KB)

Dear Lammps Users,

I have got problem with my simulation. I have run the script pdmsmix.in and Lammps have ploted ploted following informations:

[email protected]…3678…:~/Pulpit/lammps$ ./lmp_g++<pdmsmix.in
LAMMPS (12 Oct 2011)
Reading data file …
orthogonal box = (-19.9718 -20.0002 -20.077) to (80.0282 79.9998 79.923)
1 by 1 by 1 processor grid
2816 atoms
WARNING: Resetting reneighboring criteria during minimization (min.cpp:160)
Setting up minimization …
Memory usage per processor = 5.8669 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -223951.58 0 -223951.58 -2109.4451
3 0 -228485.11 0 -228485.11 139.05638
Loop time of 2.11996 on 1 procs for 3 steps with 2816 atoms

Minimization stats:
Stopping criterion = energy tolerance
Energy initial, next-to-last, final =
-223951.579018 -228359.555674 -228485.109062
Force two-norm initial, final = 5469.42 983.39
Force max component initial, final = 198.911 197.928
Final line search alpha, max atom move = 9.80582e-05 0.0194084
Iterations, force evaluations = 3 6

Pair time () = 2.10659 (99.3691) Neigh time () = 0 (0)
Comm time () = 0.000197777 (0.00932924) Outpt time () = 0 (0)
Other time (%) = 0.0131767 (0.621553)

Nlocal: 2816 ave 2816 max 2816 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 160911 ave 160911 max 160911 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 160911
Ave neighs/atom = 57.1417
Neighbor list builds = 0
Dangerous builds = 0
Setting up run …
Memory usage per processor = 4.98756 Mbytes
Step Temp E_pair E_mol TotEng Press
3 0 -228485.11 0 -228485.11 139.05638
203 4387.9468 -184014.95 0 -147195.78 880.41593
Loop time of 69.8521 on 1 procs for 200 steps with 2816 atoms

Pair time () = 68.4062 (97.9301) Neigh time () = 0.426353 (0.610364)
Comm time () = 0.00869792 (0.0124519) Outpt time () = 0.946552 (1.35508)
Other time (%) = 0.0642965 (0.0920466)

Nlocal: 2816 ave 2816 max 2816 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 35 ave 35 max 35 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 152857 ave 152857 max 152857 min
Histogram: 1 0 0 0 0 0 0 0 0 0

Total # of neighbors = 152857
Ave neighs/atom = 54.2816
Neighbor list builds = 19
Dangerous builds = 18

I have got couple of questions to this simulation and maybe someone will know anwsers:

  1. From the 1st step of simulation there are bonds between H atoms in water molecules. How is that possible?
  2. In contrast, there are no bonds between Si and C atoms in PDMS chain.
  3. Why when the time increases some of the atoms “fly away”?
  4. What does this warning means “WARNING: Resetting reneighboring criteria during minimization (min.cpp:160)”?
  5. Am I doing right this minimization?

I want to minimize the system and then observe what will happen when pressure and temperature will increase, but as you see I have problems with the minimization so I haven’t changed conditions. I would be very greatful for any comments, tips and help.
I think all necessary files are attached.
Best regards
Andrzej Battej

pdmsmixdump.zip (3.31 MB)

ffield (21.2 KB)

pdmsmix.data (139 KB)

pdmsmix.in (378 Bytes)

Replies and comments below, thanks.

Ray

1. From the 1st step of simulation there are bonds between H atoms in water
molecules. How is that possible?

I added dump cfg to your script and saw no bonded H atoms in water
molecules. You might be deceived by the visualization software you
used (i.e. by default cutoff).

2. In contrast, there are no bonds between Si and C atoms in PDMS chain.

Same as above: I did see Si-C bonds.

3. Why when the time increases some of the atoms "fly away"?

Because you had bad dynamics. I tried several NVE runs and it seemed
the force field you used only conserves with a rather small timestep
size. Default timestep is too big, and you can try yourself to
determine a good timestep.

4. What does this warning means "WARNING: Resetting reneighboring criteria
during minimization (min.cpp:160)"?

Minimization requires "neigh_modify every 1 delay 0 check yes". If
your script is not set to the above, "min" will reset it for you,
hence the warning message.

5. Am I doing right this minimization?

Your minimization is probably okay, what is not is the dynamics.

Dear Lammps Users,

I have got problem with my simulation. I have run the script pdmsmix.in and
Lammps have ploted ploted following informations:

[email protected]:~/Pulpit/lammps$ ./lmp_g++<pdmsmix.in
LAMMPS (12 Oct 2011)

ouch. that is a rather old version of LAMMPS. particularly for using
REAXFF you should update to a current version. there were quite a
few bugfixes and significant improvements, particularly to the USER-REAXC
package, which would perform much better.

[...]

Setting up run ...
Memory usage per processor = 4.98756 Mbytes
Step Temp E_pair E_mol TotEng Press
       3 0 -228485.11 0 -228485.11 139.05638
     203 4387.9468 -184014.95 0 -147195.78 880.41593
Loop time of 69.8521 on 1 procs for 200 steps with 2816 atoms

[...]

I have got couple of questions to this simulation and maybe someone will
know anwsers:
1. From the 1st step of simulation there are bonds between H atoms in water
molecules. How is that possible?

there are no bonds. this is what the visualization program shows
you due to a failed heuristic. you dump file has no information about
the atom elements or names, just numbers for the atom types. how
should a viz program know what is a hydrogen or oxygen or silicon
or carbon?? furthermore, the dump file also does not contain any
bonding information (for reax there isn't anything to dump in the first place).

2. In contrast, there are no bonds between Si and C atoms in PDMS chain.

same reason.

3. Why when the time increases some of the atoms "fly away"?

because you "created" over 4000K worth of kinetic energy.
*everything decomposes" at that temperature. in your case
the cause is obvious. you didn't choose your simulation
parameters like the time step properly. running a reax based
simulation with hydrogens in them with the default timestep
for real units will result in a totally bogus simulation. you
need to use a much shorter timestep. what is a good time
step depends on many factors. please refer to your favorite
text book on MD for a detailed explanation.

4. What does this warning means "WARNING: Resetting reneighboring criteria
during minimization (min.cpp:160)"?

for a minimization the neighborlist build rules have to be changed.
don't worry about it. lammps will restore the original settings afterwards.
check the documentation.

5. Am I doing right this minimization?

yes. the problem is later.

I want to minimize the system and then observe what will happen when
pressure and temperature will increase, but as you see I have problems with
the minimization so I haven't changed conditions. I would be very greatful
for any comments, tips and help.

there is no problem with the minimization.

axel.

Dear Axel and Ray,

My comment and questions are below. Thanks very much for the previous advices!

ouch. that is a rather old version of LAMMPS. […]

Well, I’m aware of that but I’m working on universisty and this is the Lammp’s version that is currently installed. I don’t have permission to change anything since I’m only student but I’ll talk to the administrator about upgrading lammps when he will come back from holidays. For this moment I’m forced to work on this version of Lammps.

Setting up run …
Memory usage per processor = 4.98756 Mbytes
Step Temp E_pair E_mol TotEng Press
3 0 -228485.11 0 -228485.11 139.05638
203 4387.9468 -184014.95 0 -147195.78 880.41593
Loop time of 69.8521 on 1 procs for 200 steps with 2816 atoms

[…]

I have got couple of questions to this simulation and maybe someone will
know anwsers:

  1. From the 1st step of simulation there are bonds between H atoms in water
    molecules. How is that possible?

there are no bonds. this is what the visualization program shows
you due to a failed heuristic. you dump file has no information about
the atom elements or names, just numbers for the atom types. how
should a viz program know what is a hydrogen or oxygen or silicon
or carbon??

I still have problems with visualization. I tried with this line:
“dump snap all cfg 2 pdmsmixdump.txt.*cfg id type xs ys zs”
but I find it annoying that there are several dozens of files. Is there a way to make a one file which will be understandable for AtomEye? And in addition in AtomEye all my atoms are “C” type and i don’t know why.

I tried also with this:
“dump 1 all custom 10 xxx.txt id type xs ys zs”
but VMD still have problems with understanding this.

furthermore, the dump file also does not contain any
bonding information (for reax there isn’t anything to dump in the first place). […]

Well, I undestand that in the begining there are no bonds in the system but after running a simulation, ReaxFF should compute bond orders and bonds should appeare. Maybe I don’t attach the information about bonds into the dump file? If that’s true, how to include this type of information? Now when I load the date to VMD I have bonds from the 1st second of simulation, which is weird for me becous I thought the bonds will appeare later. Correct me if I’m wrong. I would be grateful for any kind of advices connected with “dumping” the correct (and understandable for visualization software) data.

because you “created” over 4000K worth of kinetic energy.
*everything decomposes" at that temperature. in your case
the cause is obvious. you didn’t choose your simulation
parameters like the time step properly. running a reax based
simulation with hydrogens in them with the default timestep
for real units will result in a totally bogus simulation. you
need to use a much shorter timestep. what is a good time
step depends on many factors. please refer to your favorite
text book on MD for a detailed explanation.

Thanks for the advice. I’ve changed the timestep to 0.2 and now everything looks good.

Step Temp E_pair E_mol TotEng Press
3 0 -228485.11 0 -228485.11 139.05638
203 299.18188 -231017.52 0 -228507.09 -235.05011

As you see the temperature is now circa 26C (almost room temperature). I think that would be correct.

Comments below.

Best,
Ray

Well, I'm aware of that but I'm working on universisty and this is the
Lammp's version that is currently installed. I don't have permission to
change anything since I'm only student but I'll talk to the administrator
about upgrading lammps when he will come back from holidays. For this moment
I'm forced to work on this version of Lammps.

Why not install an up-to-date version of LAMMPS on your own machine/laptop?

I still have problems with visualization. I tried with this line:
"dump snap all cfg 2 pdmsmixdump.txt.*cfg id type xs ys zs"
but I find it annoying that there are several dozens of files. Is there a
way to make a one file which will be understandable for AtomEye? And in
addition in AtomEye all my atoms are "C" type and i don't know why.

Use dump_modify to change element types. Having multiple files is
require by AtomEye; I don't think there is a way out.

Well, I undestand that in the begining there are no bonds in the system but
after running a simulation, ReaxFF should compute bond orders and bonds
should appeare. Maybe I don't attach the information about bonds into the
dump file? If that's true, how to include this type of information? Now when
I load the date to VMD I have bonds from the 1st second of simulation, which
is weird for me becous I thought the bonds will appeare later. Correct me if
I'm wrong. I would be grateful for any kind of advices connected with
"dumping" the correct (and understandable for visualization software) data.

ReaxFF is a pure atomic force field, hence there are no "bonds".
ReaxFF does computes bond orders and considers bonded interactions
between atoms, but "bonds" are just atoms with bond orders larger than
a specific bond order cutoff or with distances smaller than a specific
distance cutoff. You have to tell your vis software how to draw bonds
between atoms if you think there should exist a bond. There is
absolutely no information on "bonds" that you should/can write to a
dump file. Look at the AtomEye or VMD instructions on how to set
distance cutoffs for visualizing atoms.

Thanks for the advice. I've changed the timestep to 0.2 and now everything
looks good.

Step Temp E_pair E_mol TotEng Press
       3 0 -228485.11 0 -228485.11 139.05638
     203 299.18188 -231017.52 0 -228507.09 -235.05011

As you see the temperature is now circa 26C (almost room temperature). I
think that would be correct.

A final step temperature after only 200 steps of NVE really tells you
nothing. You have to run longer, or with NVT to obtain correct
statistics.

Dear Axel and Ray,

My comment and questions are below. Thanks very much for the previous
advices!

ouch. that is a rather old version of LAMMPS. [...]

Well, I'm aware of that but I'm working on universisty and this is the
Lammp's version that is currently installed. I don't have permission to
change anything since I'm only student but I'll talk to the administrator
about upgrading lammps when he will come back from holidays. For this moment
I'm forced to work on this version of Lammps.

this is complete nonsense. lammps doesn't need anything that
requires administrator privileges. such software can be compiled
by yourself and you can place the executable where you can
store your data. i am doing it all the time and on many different
machines (including some of the largest supercomputers available
to researchers) and never needed an administrator to do anything
about that.

[...]

I still have problems with visualization. I tried with this line:
"dump snap all cfg 2 pdmsmixdump.txt.*cfg id type xs ys zs"
but I find it annoying that there are several dozens of files. Is there a
way to make a one file which will be understandable for AtomEye? And in
addition in AtomEye all my atoms are "C" type and i don't know why.

this is how atomeye works. you would have to use dump_modify
element to set the element names as needed.

I tried also with this:
"dump 1 all custom 10 xxx.txt id type xs ys zs"
but VMD still have problems with understanding this.

that is creating essentially the same file
as the atom style dump. the issue is that
there is no information about the elements.
type is just a number. VMD has no psychic
powers and knows what "type" is supposed
to stand for.

*however*, what you *can* do is to create/use
a psf file that matches your element names.
your data file, since it was created with VMD
has almost all the necessary information.

running the following commands on the VMD
console creates a .psf file that can be used
in combination with the dump files to convey
the necessary information to VMD:

topo readlammpsdata pdmsmix.data charge
set sel [atomselect top all]
$sel set segname PDMS
animate write psf pdmsmix.psf
quit

furthermore, the dump file also does not contain any
bonding information (for reax there isn't anything to dump in the first
place). [...]

Well, I undestand that in the begining there are no bonds in the system but
after running a simulation, ReaxFF should compute bond orders and bonds
should appeare. Maybe I don't attach the information about bonds into the
dump file? If that's true, how to include this type of information? Now when

yes, reax computes those properties. but those are not "bonds" in the
LAMMPS sense, i.e. defined as explicit force field interactions.
you can output the bond information that is computed by reax/c using this fix.
http://lammps.sandia.gov/doc/fix_reaxc_bonds.html

I load the date to VMD I have bonds from the 1st second of simulation, which
is weird for me becous I thought the bonds will appeare later. Correct me if
I'm wrong. I would be grateful for any kind of advices connected with
"dumping" the correct (and understandable for visualization software) data.

you have to spend some time with the VMD tutorials.
a visualization tool shows bond not because of any knowledge
of the bonding situation, but based on some rather simple
heuristics. also, in the case of VMD, bonding is considered
as in forcefield calculations, i.e. it usually imports bonds
from a .psf file. short of that, it uses a heuristic guess and
then *keeps* those bonds. the entire topology information
is only stored once, as for a regular force field calculation,
this doesn't change and hence a lot of storage can be
saved to allow for storing longer trajectories.

in short, there is no simple way. you can use the simple
heuristics of a viz program (but you have to feed it the
missing information, e.g. what types of atoms you have,
and work around cases where it gets it wrong) or use
your own heuristics or take the data stored by lammps
and do some scripting. you will have to spend some time
on this and develop some tools. i can provide some
suggestions for VMD, but there exists no perfect solution,
since reax is a new concept compared to traditional
force field calculations and there are only a small number
of people that use this. just do a little research on
the problem that people have with visualization data
from ab initio simulations. i had an entire section
of my home page dedicated to it, when i was still
working on/with some of those codes. you can find
the remnants here:

https://sites.google.com/site/akohlmey/redirect/cpmd-vmd.pdf?attredirects=0&d=1

because you "created" over 4000K worth of kinetic energy.
*everything decomposes" at that temperature. in your case
the cause is obvious. you didn't choose your simulation
parameters like the time step properly. running a reax based
simulation with hydrogens in them with the default timestep
for real units will result in a totally bogus simulation. you
need to use a much shorter timestep. what is a good time
step depends on many factors. please refer to your favorite
text book on MD for a detailed explanation.

Thanks for the advice. I've changed the timestep to 0.2 and now everything
looks good.

Step Temp E_pair E_mol TotEng Press
       3 0 -228485.11 0 -228485.11 139.05638
     203 299.18188 -231017.52 0 -228507.09 -235.05011

As you see the temperature is now circa 26C (almost room temperature). I
think that would be correct.

well, i guess that is by chance, since you were using NVE integration.
and what you see is just the result of the potential energy from
changing to reax

cheers,
    axel.