Hi! I’m simulating a system with ozone and organic molecules using the ReaxFF potential [1]. But my results show that the bonds between the carbons of the aromatic ring are completely broken, forming a system with several dispersed molecules, instead of an ozone attack on the molecule, which was expected.
I assembled the xyz system using packmol and used OVITO to generate the system data file.
I will be placing the input, data and potential files here.
Thanks for all the help!
INPUT
units real
atom_style full
boundary p p p
read_data gua_ozo.data
mass 1 15.9994 # O
mass 2 12.0107 # C
mass 3 1.00794 # H
pair_style reaxff NULL safezone 3.0 mincap 150
pair_coeff * * CHO-2016.ff C H O
compute reax all pair reaxff
fix myqeq all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff maxiter 400
group gO type 1
group gC type 2
group gH type 3
variable qO equal charge(gO)/count(gO)
variable qC equal charge(gC)/count(gC)
variable qH equal charge(gH)/count(gH)
thermo 10
thermo_style custom step temp etotal press vol v_qO v_qC v_qH
dump dmp all custom 50 dump.lammpstrj id type q x y z
fix myspec all reaxff/species 5 1 5 species.log element O C H
minimize 1.0e-5 1.0e-7 10000 100000
#thermo 50
velocity all create 300 4928459 mom yes rot yes dist gaussian
#fix 1 all setforce 0 0 0
#fix mynpt all npt temp 300.0 300.0 100 iso 1.0 1.0 1000
fix mynvt all nvt temp 300.0 300.0 100
timestep 0.01
run 10000
write_data system_relaxed.data
DATA
LAMMPS data file written by OVITO Basic 3.11.3
2350 atoms
3 atom types
-5 65 xlo xhi
-5 65 ylo yhi
-5 65 zlo zhi
Masses
1 15.9994 # O
2 12.0107 # C
3 1.00794 # H
Atoms # full
Here goes all atoms locations
Reactive MD-force field: Combustion C/H/O force field - September 7, 2016
39 ! Number of general parameters
50.0000 !p(boc1)
9.5469 !p(boc2)
26.5405 !p(coa2)
0.6863 !p(trip4)
2.7295 !p(trip3)
70.0000 !kc2
1.0588 !p(ovun6)
4.1262 !p(trip2)
12.1176 !p(ovun7)
13.3056 !p(ovun8)
-68.9784 !p(trip1)
0.0000 !Lower Taper-radius (swa)
10.0000 !Upper Taper-radius (swb)
0.0000 !not used
33.8667 !p(val7)
6.0891 !p(lp1)
1.0563 !p(val9)
2.0384 !p(val10)
6.1431 !not used
6.9290 !p(pen2)
0.3989 !p(pen3)
3.9954 !p(pen4)
0.0000 !not used
5.7796 !p(tor2)
10.0000 !p(tor3)
1.9487 !p(tor4)
0.0000 !not used
2.1645 !p(cot2)
1.5591 !p(vdW1)
0.1000 !Cutoff for bond order*100 (cutoff)
2.1365 !p(coa4)
0.6991 !p(ovun4)
50.0000 !p(ovun3)
1.8512 !p(val8)
0.0000 !not used
0.0000 !not used
0.0000 !not used
1.0000 !not used
2.6962 !p(coa3)
3 ! Nr of atoms; atomID;ro(sigma); Val;atom mass;Rvdw;Dij;gamma
alfa;gamma(w);Val(angle);p(ovun5);n.u.;chiEEM;etaEEM;n.u.
ro(pipi);p(lp2);Heat increment;p(boc4);p(boc3);p(boc5),n.u.;n.u.
p(ovun2);p(val3);n.u.;Val(boc);p(val5);n.u.;n.u.;n.u.
C 1.3674 4.0000 12.0000 2.0453 0.1444 0.9500 1.1706 4.0000
9.0000 1.5000 4.0000 27.5134 79.5548 5.0191 7.0500 0.0000
1.1168 0.0000 181.0000 14.2732 24.4406 6.7313 0.8563 0.0000
-4.1021 5.0000 1.0564 4.0000 2.9663 0.0000 0.0000 0.0000
H 0.9479 1.0000 1.0080 1.1364 0.0232 0.9900 -0.1000 1.0000
9.0643 4.7746 1.0000 0.0000 121.1250 4.7757 9.7732 1.0000
-0.1000 0.0000 62.4879 2.5194 2.3785 0.2223 1.0698 0.0000
-15.7683 2.1488 1.0338 1.0000 2.8793 0.0000 0.0000 0.0000
O 1.1939 2.0000 15.9990 1.9289 0.1201 0.9900 1.0981 6.0000
10.4842 8.2916 4.0000 28.8967 116.0768 7.9703 7.0500 2.0000
1.0479 20.0000 60.8726 10.0338 2.2024 0.9942 0.9745 0.0000
-3.6141 2.7025 1.0493 4.0000 2.9225 0.0000 0.0000 0.0000
6 ! Nr of bonds; at1;at2;De(sigma);De(pi);De(pipi);p(be1);p(b
p(be2);p(bo3);p(bo4);n.u.;p(bo1);p(bo2)
1 1 80.8865 107.9944 52.0636 0.5218 -0.3636 1.0000 34.9876 0.7769
6.1244 -0.1693 8.0804 1.0000 -0.0586 8.1850 1.0000 0.0000
1 2 179.5195 0.0000 0.0000 -0.5242 0.0000 1.0000 6.0000 0.7187
5.4740 1.0000 0.0000 1.0000 -0.1144 6.7029 0.0000 0.0000
2 2 113.9232 0.0000 0.0000 -0.5971 0.0000 1.0000 6.0000 0.9093
1.7152 1.0000 0.0000 1.0000 -0.0450 6.0710 0.0000 0.0000
1 3 136.4945 164.1201 5.5000 -0.9159 -0.1075 1.0000 10.6519 0.8644
0.6858 -0.4602 9.5754 1.0000 -0.1745 4.5987 0.0000 0.0000
3 3 148.0798 155.2406 20.1160 -1.0000 -0.1254 1.0000 33.0027 0.7790
0.7673 -0.1697 7.0028 1.0000 -0.1300 5.1959 1.0000 0.0000
2 3 169.1351 0.0000 0.0000 -0.8810 0.0000 1.0000 6.0000 0.5757
1.5482 1.0000 0.0000 1.0000 -0.1788 4.6622 0.0000 0.0000
3 ! Nr of off-diagonal terms. at1;at2;Dij;RvdW;alfa;ro(sigma);r
1 2 0.1253 1.5717 9.9736 1.2057 -1.0000 -1.0000
2 3 0.1125 1.6311 8.7528 1.0929 -1.0000 -1.0000
1 3 0.0953 1.7397 8.8986 1.4256 1.1067 1.1265
18 ! Nr of angles. at1;at2;at3;Thetao,o;p(val1);p(val2);p(coa1);
1 1 1 76.1370 34.6920 1.1328 0.0000 0.0050 0.3556 1.8065
1 1 2 68.0572 9.9461 4.7000 0.0000 0.4566 0.0000 1.8532
2 1 2 65.6815 35.0000 1.8622 0.0000 0.0490 0.0000 1.0937
1 2 2 0.0000 4.0000 7.2043 0.0000 0.0000 0.0000 1.0728
1 2 1 0.0000 3.4110 7.7350 0.0000 0.0000 0.0000 1.0400
2 2 2 0.0000 30.0000 5.6235 0.0000 0.0000 0.0000 1.0400
1 1 3 78.3624 13.0773 9.0480 0.0000 0.1270 52.1129 2.3964
3 1 3 76.7101 24.3833 5.8613 -21.8559 2.6395 -32.6534 3.6179
2 1 3 79.1288 30.0000 1.4632 0.0000 0.2065 0.0000 2.0000
1 3 1 80.7352 16.4130 4.9987 0.0000 0.0843 0.0000 1.0137
1 3 3 85.4436 14.4937 3.9928 0.0000 1.4350 44.5320 1.1348
3 3 3 89.9282 32.1199 2.7181 0.0000 0.3323 57.6122 1.0000
1 3 2 82.9640 32.4874 0.8777 0.0000 0.9627 0.0000 1.0010
2 3 3 85.7838 17.3139 1.9157 0.0000 3.6306 0.0000 2.1858
2 3 2 84.2527 33.1226 0.6730 0.0000 0.7238 0.0000 2.4348
1 2 3 0.0000 14.4588 3.1507 0.0000 3.4571 0.0000 1.0149
3 2 3 0.0000 0.9696 3.6303 0.0000 0.0000 0.0000 1.6987
2 2 3 0.0000 0.5797 1.9739 0.0000 0.0000 0.0000 2.4494
21 ! Nr of torsions. at1;at2;at3;at4;;V1;V2;V3;p(tor1);p(cot1);n
1 1 1 1 2.0474 32.6719 0.5282 -9.0000 -2.6449 0.0000 0.0000
1 1 1 2 1.6328 78.4995 -0.1514 -6.9161 -2.9986 0.0000 0.0000
2 1 1 2 2.4142 78.7025 0.3506 -8.8640 -6.9283 0.0000 0.0000
1 1 1 3 -0.7104 22.6038 0.5309 -2.0000 -0.6614 0.0000 0.0000
2 1 1 3 1.9323 52.9368 0.6554 -8.8118 -3.9854 0.0000 0.0000
3 1 1 3 -1.2500 1.1248 -0.1230 -9.9453 -3.9000 0.0000 0.0000
1 1 3 1 -0.6848 56.7751 -1.2733 -2.2937 -4.0000 0.0000 0.0000
1 1 3 2 -1.4557 78.6279 0.9945 -3.2742 -2.4240 0.0000 0.0000
2 1 3 1 0.6928 78.1546 0.5608 -3.1713 -3.7301 0.0000 0.0000
2 1 3 2 -1.4343 77.0699 0.9875 -3.4139 -1.4053 0.0000 0.0000
1 1 3 3 0.5153 2.1584 0.2000 -6.5859 -3.0000 0.0000 0.0000
2 1 3 3 0.2018 80.0000 0.3778 -2.5000 -2.8750 0.0000 0.0000
3 1 3 1 -1.9875 79.2591 1.0000 -2.4206 -3.9342 0.0000 0.0000
3 1 3 2 -1.1000 78.8002 -1.0000 -2.6282 -4.0000 0.0000 0.0000
3 1 3 3 -1.0000 83.5323 4.3660 -2.6805 -1.2938 0.0000 0.0000
1 3 3 1 3.4682 0.0781 0.9887 -6.1195 -0.5004 0.0000 0.0000
1 3 3 2 1.0000 16.5478 -1.0313 -2.0000 -2.6888 0.0000 0.0000
2 3 3 2 4.0818 -3.2744 -0.9664 -7.1634 -3.0000 0.0000 0.0000
1 3 3 3 4.2014 -10.0642 1.8690 -2.4805 -2.5000 0.0000 0.0000
2 3 3 3 1.0000 -10.0500 -1.0000 -2.1946 -0.5300 0.0000 0.0000
3 3 3 3 1.0000 1.6871 3.0000 -6.2660 -0.5500 0.0000 0.0000
1 ! Nr of hydrogen bonds. at1;at2;at3;r(hb);p(hb1);p(hb2);p(hb3
3 2 3 1.8130 -3.5409 2.3815 21.9463
This looks like you are trying to use a combustion force field for a room temperature molecular system. Not likely to work.
Thanks! I didn’t pay attention to that
I found the problem, it was the order of the atoms after the file name of the reaxff potential