weird species when using fix nphug to simulate shock

Dear All,
After I use fix nphug (with ReaxFF-lg) to simulate the effect of shock wave on HMX molecules , I found the generated species are quite different from experiment because there are many big molecules. I used new cutoff value instead of the default one. Now I am wondering where is the problem? I’ll appreciate your help if anyone can offer some advice. My input script is as follow:

fix 2 all qeq/reax 1 0.0 10.0 1e-6 reax/c
fix 3 all reax/c/species 1 10 100 species.out cutoff 1 4 0.3 cutoff 1 1 0.55 cutoff 1 3 0.65 cutoff 1 2 0.4 cutoff 3 3 0.65 cutoff 4 3 0.4 cutoff 3 2 0.4 cutoff 2 2 0.55 cutoff 2 4 0.55 cutoff 4 4 0.55

thermo 10
thermo_style custom step vol press temp etotal lx ly lz pxx pyy pzz
#thermo_modify lost ignore flush yes

#timestep 0.2

dump 1 all custom 1000 dump.reax.HMX id type q x y z

#run 10000
#write_restart HMX.afternvt.restart

#unfix 1

fix 1 all nphug temp 1.0 1.0 100.0 x 680000 680000 1000.0
#fix 1 all npt temp 3500.0 3500.0 100.0 x 395000.0 395000.0 1000.0
#fix 3 all deform 1

timestep 0.02
run 200000

The generated species is as follow:

799400 71 27 1 1 1 1 1 1 9 3 1 1 4 1 1 3 2 1 1 1 1 1 3 21 5 1 1 1 3

Timestep No_Moles No_Specs C8H16O17N17 C8H16O15N11 C4H7O8N8 C12H24O21N18 C24H48O45N40 C16H31O26N19 C4H8O8N8 C4H8O6N5 C12H23O23N22 C4H8O9N10 C4H9O7N5 C4H8O8N7 C9H17O14N11 C47H94O92N91 C8H16O14N12 C4H9O8N8 C4H8O7N6 C4H8O8N9 C4H8O8N6 C4H7O7N7 C8H15O15N12 C4H8O9N11 C4H8O7N5 HN HON2 ON2 ON3 ON4 O2N4 N2

799500 68 30 1 1 1 1 1 1 7 1 1 1 1 4 1 1 1 1 4 1 2 1 1 1 1 1 3 19 4 1 1 3

Timestep No_Moles No_Specs C36H72O71N68 C8H16O15N11 C4H7O8N8 C12H24O21N18 C16H31O26N19 C4H8O8N8 C4H8O6N5 C47H93O92N90 C4H8O9N10 C4H9O7N5 C4H8O8N7 C9H17O14N11 C8H16O15N15 C8H16O13N10 C4H9O8N8 C4H8O7N6 C4H8O8N9 C4H7O7N7 C8H15O16N14 C4H8O7N4 C4H8O7N5 HN HON2 ON2 ON3 ON4 O2N4 N2

799600 66 28 1 1 1 1 1 9 1 1 1 1 4 1 1 1 2 3 1 1 1 1 1 1 2 20 4 1 1 2

Timestep No_Moles No_Specs C24H48O45N41 C8H17O15N12 C4H7O8N8 C12H24O21N18 C16H31O26N19 C4H8O8N8 C8H16O16N17 C52H103O99N97 C4H9O7N5 C4H8O8N7 C9H17O14N11 C4H8O7N5 C4H8O9N9 C4H8O7N6 C8H16O12N8 C3H6O7N6 C4H9O8N8 C4H8O8N9 C4H7O7N7 C8H15O16N14 C4H8O7N4 C4H8O9N11 HON2 ON2 ON3 ON4 O2N4 N2

799700 70 28 1 1 1 1 1 8 1 1 1 4 1 2 1 4 1 1 2 1 1 1 1 1 2 23 4 1 1 2

Timestep No_Moles No_Specs C16H32O29N24 C8H17O15N12 C4H7O8N8 C10H20O17N15 C16H31O26N19 C4H8O8N8 C8H16O16N17 C52H103O99N97 C8H16O15N15 C4H9O7N5 C4H8O8N7 C9H17O14N11 C4H8O9N9 C4H8O7N6 C8H16O13N10 C3H6O7N6 C4H9O8N8 C4H8O7N5 C4H8O8N9 C4H7O7N7 C8H14O14N10 C4H8O7N4 C4H8O6N6 C4H8O9N11 C2H4O4N3 HON2 HO2N4 ON2 ON3 ON4 O2N4 O2N N2

799800 75 33 1 1 1 1 1 7 1 1 1 1 4 1 1 5 1 1 1 2 1 1 1 1 1 1 1 3 1 22 5 1 1 1 2

Timestep No_Moles No_Specs C12H24O22N19 C8H17O15N12 C4H7O7N6 C12H24O21N17 C20H39O32N24 C4H8O8N8 C44H89O85N86 C12H23O25N26 C8H16O15N15 C4H9O7N5 C4H8O8N7 C4H8O7N5 C9H17O14N11 C4H8O7N6 C4H8O9N9 C8H16O13N10 C3H6O7N6 C4H9O8N8 C4H8O8N9 C4H8O8N6 C4H7O6N5 C8H15O16N14 C4H8O7N4 C8H16O16N16 HON2 ON2 ON3 ON4 O2N4 O2N N2 N

799900 76 32 1 1 1 1 1 7 1 1 1 1 3 2 1 6 1 1 1 1 1 1 1 1 1 1 2 24 4 1 1 1 3 2

Timestep No_Moles No_Specs C20H40O38N35 C8H17O15N12 C4H7O7N6 C12H24O21N18 C4H8O6N4 C4H8O8N8 C28H55O50N47 C12H23O25N26 C8H16O16N17 C4H9O7N5 C4H8O8N7 C16H31O26N21 C4H8O7N5 C12H23O20N16 C4H8O7N6 C16H34O35N39 C4H8O9N9 C8H16O12N8 C4H9O8N8 C4H8O8N6 C4H7O7N7 C8H15O16N14 C4H8O7N4 HON2 ON2 ON ON3 ON4 O2N4 O2N N2

800000 74 31 1 1 1 1 1 8 1 1 1 1 3 1 1 1 7 1 1 1 1 1 1 1 1 2 23 1 4 1 1 1 3

And where is the problem?

Oleg

WangGuangyu <[email protected]…8…> 31 октября 2014 г. 23:17:54 написал:

Oleg,
Check out the generated species. There are many big molecules (like C16H32O29N24). It is weird because explosive will decompose after detonation. Thus it is weird to have so many big molecules.
Guangyu Wang

As far as I can see, your calculation is 800000 * 0.02 fs = 16 ps long. I’m not sure it is enough for all the chemistry taking place during “detonation”.

Second, there is no guarantee that any parametrization gives the experimental kinetics or correct reaction mechanisms.

Besides, the set of “threshold” bond orders you use is still arbitrary. You may try to tune it further, and if you obtain good agreement with the experimental data, this new set will be very interesting for many people, I suppose.

In other words, this problem does not have the one correct solution.

Oleg

WangGuangyu <[email protected]…8…> 31 октября 2014 г. 23:39:50 написал:

Oleg,
Thanks a lot for your advice. I got the cut off value from another paper. Now I am trying to extend the simulation time to see if there is any difference. I have another question. In the setting of cutoff value, does 1 mean C, 2 mean H, 3 mean O, 4 mean N according to the ReaxFF file or the data file determining the structure of molecules?Because I defined the position of atoms in the order C H N O in my input file. Thanks for reply.
Guangyu Wang

Guangyu,

the correspondence between type number, mass and positions of atoms of this type is set while reading data file. Types in pair_coeff command in input script should then be placed in the same order (see rdx example in lammps distribution).

Oleg

01.11.2014, 00:52, "WangGuangyu" <[email protected]>: