Hi there,
I want to simulate the pyrolysis of polyethylene using ReaxFF potential in LAMMPS. In Materials Studio, I built the amorphous PE structure of 20 chains with each chain having 64 carbon atoms, performed a geometric optimization using Forcite module, and finally converted it into a data file for use in LAMMPS. After importing into LAMMPS, it was firstly equilibrated for 3 ns under CVFF potential field to get a better structure.
However, when I conducted a pyrolysis simulation with it, lots of warnings like “Fix qeq/reaxff CG convergence failed after 200 iterations at step 1 (src/REAXFF/fix_qeq_reaxff.cpp:774)” appeared, and finally the simulation stopped with the error “Non-numeric atom coords - simulation unstable (src/OPENMP/domain_omp.cpp:58)”.
I know that such errors usually occur due to bad geometric configuration, so I equilibrate my system for 3 ns to obtain a better structure. But it still failed. I also try to use the PE structure in the example directory of Materials Studio, but the same error still occurred.
Could you please give me some suggestions on how to solve this problem. The attachments are my input file, log file and data file that is converted from the Materials Studio example PE structure.
MS-PE.data (1.4 MB)
pyrolysis.in (736 Bytes)
log.lammps (16.4 KB)
Please see these warnings in your log files:
WARNING: Bonds are defined but no bond style is set (src/force.cpp:192)
WARNING: Likewise 1-2 special neighbor interactions != 1.0 (src/force.cpp:194)
WARNING: Angles are defined but no angle style is set (src/force.cpp:197)
WARNING: Likewise 1-3 special neighbor interactions != 1.0 (src/force.cpp:199)
WARNING: Dihedrals are defined but no dihedral style is set (src/force.cpp:202)
WARNING: Likewise 1-4 special neighbor interactions != 1.0 (src/force.cpp:204)
WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (src/comm_brick.cpp:209)
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
They are a sign of serious trouble. The ReaxFF potential models bonds, angles, and dihedrals not explicitly from the topology, but implicitly on the fly from the bond order. So there should not be any bonds defined in the bond topology. If you do have them defined (they are easily deleted from the data file with a text editor, BTW), then you must ensure that all pair wise non-bonded interactions are considered. That is controlled with the special_bonds command. In your case you need a setting of special_bonds lj/coul 1.0 1.0 1.0 while the default setting is to exclude all pairwise non-bonded interactions for 1-2, 1-3, and 1-4 pairs that are involved in bonds.
in addition:
fix 2 all nve
fix 3 all nvt temp 1000.0 1000.0 25
These two lines are an error. Please read the documentation about fix nvt with more care and pay attention to what it says about time integration of atoms.
Thank you very much!
I’ve addressed the problem according to your kind suggestion. And modify the commands for thermostatting.
I do appreciate your help!