I am trying to reproduce some results from this article (Agrawalla, S.; van Duin, A. C. T. J. Phys. Chem. A 2011, 115, 960-972) which uses reaxff to simulate the combustion of H2/O2 to form water. I have tried both the ffield.reax.cho potential and the H2 potential include in the supplemental information of the paper, but it seems to give erroneous results when compared to the results of the paper. In particular, using the ffield.reax.cho potential seems to allow the reactions to occur at an earlier time step. I also tried the optimized potential from the paper, but those results are different as well.
1b) If I have to alter the lone pair potential function where do I do this? ( is it in the the code for the Reax program or in the potential from the supplemental information?)
Thanks in Advance,
supp_info_H2_optimized_potential.txt (10.2 KB)
in.reax.H2 (1.24 KB)
1a) In the paper they say that the H2 potential can be readily used in the
stand alone ReaxFF program, but that some alteration of the lone pair
potential function is required to run with LAMMPS. I am using the 11Jan12
version of LAMMPS, so is it possible that this particular potential
function has already been modified?
No, current pair_style reax and reax/c have not yet incorporated this change.
1b) If I have to alter the lone pair potential function where do I do this?
( is it in the the code for the Reax program or in the potential from the
It is in reax_poten.F for reax and reaxc_ffield.cpp for reax/c.
2) I have been using Packmol to generate a random cell, I then modify the
output to the correct format for the data.H2 file for LAMMPS, do you think
some form of minimization and/or equilibration should be done first, before
the actual MD run? Packmol insures that all molecules are at least 2.0
Angstroms apart and that no atoms end up on top of each other.
Yes, you should first minimize the structure.
I also attached an example of my input file. Where I had to rename
whichever potential that I used to ffield.reax.
Without the data file, nobody can run your script.
Thanks for clarifying that for me. I have attached my data file for this
run, and the field.reax(it is the potential from the paper). Using the
field.reax.cho potential gave results where the the reactions happened at
much earlier time steps( I would have produced 40 waters much sooner than
Also here is a clip from my fra.out file:
Timestep H2 O2 H H2O2 HO2 HO O O4
H2O H3O H2O4 HO3 H4O H3O3 H3O2 H2O3 H3
601150 22 7 0 0 1 9 1 0 40 0 0 0 0
0 0 0 0 0
601175 22 7 0 0 1 9 1 0 40 0 0 0 0
0 0 0 0 0
601200 22 7 0 0 1 9 1 0 40 0 0 0 0
0 0 0 0 0
So at time step 601150 = 120ps I already have 40 waters, which happens a
lot earlier than any of the trials in the paper. I have done 3 different
trials(with different starting geometries), and get similar results after
averaging them up.
In terms of changing the potential, do you have any tips or suggestions
about what would be the safest way to do this, in order to keep from
screwing things up. I would hate to change it and either it only works
for H2/O2, or it give unreliable results.
data.H2_67_30 (14.3 KB)
ffield.reax (10.1 KB)
Try with two subroutines, which incorporates the lone pair change.
You will have to use pair_style reax/c with the optimized force field.
Replace the original files with these two and re-compile.
Please let me know if this works.
reaxc_types.h (18.9 KB)
reaxc_ffield.cpp (24.3 KB)
Thanks a lot. I will let you know.