Dear lammps users,
I'm computing elastic constants of diamond using ELASTIC example in lammps distribution (May 5, 2012 distribution). I try to employ two different interaction models - AIREBO and reax/c. For doing this I change the interaction in the potential.mod file
for AIREBO:
pair_style airebo 3.0
pair_coeff * * CH.airebo C
for reax/c
pair_style reax/c lmp_control
pair_coeff * * ffield.reax.cho 1
fix 2 all qeq/reax 1 0.0 10.0 1e-6 param.qeq
The system consists of 1000 carbon atoms in diamond lattice with the lattice constant 3.57A. The initial coordinates are read from the same file in both cases.
However, the elastic constants are very different in both cases. For example, C11all = 1118.890519 GPa for AIREBO while C11all = 2749339.15 GPa for reax. The value for AIREBO is in good correspondence with literature data and reax/c gives horribly wrong results.
Does it mean that reax is not appropriate for elastic constant calculations? In the publication where the reax force field for hydrocarbon oxidation is introduced ( J. Phys. Chem. A 2008, 112, 1040-1053) there is nothing about elastic properties. As I understand from this article only quantum calculations of energies were used for fitting the model, so have any one tried to see if elastic properties are reproduced well with it?
Thanks in advance for help!
Kind regards,
Yury Fomin
Hi Yury,
You would need to use default control settings. Make sure in your
lmp_control control file all settings are set to the default values,
as given in the pair reax/c doc page:
http://lammps.sandia.gov/doc/pair_reax_c.html.
Another choice would be to update to the most recent distribution and
do not use any control file.
pair styles airebo and reax/c re in different LAMMPS units. Make sure
the units are properly changed.
fix qeq/reax used with pair_style reax or reax/c should not need to
assign parameter file. Use "reax/c" instead of "param.qeq"
Lastly, if ECs are not published, maybe you need to consult the
authors of the paper for their values.
Dear Ray,
Thank you for your reply. Of course, I took care about the default settings in lmp_control file and about units convertion. I also tried to use the latest lammps version and it also didn't help. May be it's a problem of the model itself - can be that it reproduces well chemical reactions (for what it was designed), but not elastic moduli.
Kind regards,
Yury Fomin
25.07.2012, 19:24, "Ray Shan" <[email protected]...>:
Why don't you ask Adri's group what they think of your results?
Steve
The elastic constants are off by a factor of 1000. This is too large to be
caused by the potential, but is probably due to the elastic constant
calculation not converging. You should get a second opinion by computing
an energy/volume curve for one of the deformation modes.
Aidan
Hi,
Calculation of the elastic constant is very senstive to the percentage of deformation . i.e You have to be sure that you’re in the linear regime and Hook’s law works, so the smaller the strain the better chance that you’ll not be mistaken. BUT! Small values of strain means less accuracy. Since you calculate cijkl =dsigma_ij/depsilon_kl. just for Fun, You should try different percentage of deformation (between 0.2-1.2%)
and compare results
Oscar G.