Wrong Graphene/Graphite Elastic Constants with AIREBO?

Hello All,

I have been trying to do molecular mechanics simulations involving graphene sheets. To start off, I tried to find the elastic constants of graphene.

To obtain the elastic constants, I used the input script attached. This is based on the input script provided in the examples/ELASTIC folder in the lammps directory (I stripped it down quite a bit, but it essentially does the same). This gives a result of about 940 GPa for C11 and 245 GPa for C12 with LAMMPS version 13 June, 2013. But according to the reference in the LAMMPS documentation, the values should be 1150 GPa and 150 GPa respectively. These also seem to be the more established values of elastic moduli for graphene/graphite in the literature as well.

What I noticed is that immediately after change_box, LAMMPS seems to give the correct values. But upon minimization, the atoms seem to move and relax in the y-direction (the box dimension is changed in the x direction) while at the same time the box length in the y-direction is unaffected. This (atoms moving in the minimization step) does not happen in the case of the Silicon example provided. Obviously if I fix the atoms after change_box by setting the force in the y direction to zero, I get the desired result.

I have tried several permutations and combinations too with pretty much the same or equivalent result. Can anyone kindly explain what is happening here? This has been puzzling me for quite a while but it seems like I am missing something very obvious. Any help is much appreciated.

Thanks,
Narasimha.

in.elastic (2.22 KB)

Hello All,

I have been trying to do molecular mechanics simulations involving graphene
sheets. To start off, I tried to find the elastic constants of graphene.

To obtain the elastic constants, I used the input script attached. This is
based on the input script provided in the examples/ELASTIC folder in the
lammps directory (I stripped it down quite a bit, but it essentially does
the same). This gives a result of about 940 GPa for C11 and 245 GPa for C12
with LAMMPS version 13 June, 2013. But according to the reference in the
LAMMPS documentation, the values should be 1150 GPa and 150 GPa
respectively. These also seem to be the more established values of elastic
moduli for graphene/graphite in the literature as well.

What I noticed is that immediately after change_box, LAMMPS seems to give
the correct values. But upon minimization, the atoms seem to move and relax
in the y-direction (the box dimension is changed in the x direction) while
at the same time the box length in the y-direction is unaffected. This
(atoms moving in the minimization step) does not happen in the case of the
Silicon example provided. Obviously if I fix the atoms after change_box by
setting the force in the y direction to zero, I get the desired result.

I have tried several permutations and combinations too with pretty much the
same or equivalent result. Can anyone kindly explain what is happening here?
This has been puzzling me for quite a while but it seems like I am missing
something very obvious. Any help is much appreciated.

have you checked for finite size effects?
i.e. tried with different size systems and see, if there is a
(converging) trend?

axel.

Yes, I did try different size systems (as big as 100 A). I get the same results, nothing different.

Thanks,
Narasimha.

Dear Axel, Narasimha
Whether this problem is related to with in adequacy of AIREBO potential?. Recently Lindsay et al (PRB,81,205441(2010) optimized the tersoff and rebo potential parameters which yields the acoustic phonon velocities that were in good agreement with experimental values. But lammps didn’t incorporate this changes in the CH.airebo potential file.if lammps release the AIREBO potential file with the recommendation of above paper one can test and verify it.
Regards
Anees

Dear Axel, Narasimha
Whether this problem is related to with in adequacy of AIREBO potential?.
Recently Lindsay et al (PRB,81,205441(2010) optimized the tersoff and rebo
potential parameters which yields the acoustic phonon velocities that were
in good agreement with experimental values. But lammps didn't incorporate
this changes in the CH.airebo potential file.if lammps release the AIREBO
potential file with the recommendation of above paper one can test and
verify it.

two comments on that:

1) a different parameterization is just that. you can adjust
parameters and you get different results. since these are empirical
potentials, some results will match experimental data better, and some
will be worse. people have been arguing over the "right" way to
simulate water for over 30 years now and still there no single answer
to satisfy all. it is often just as well to know the deficiencies of a
(simple) model than have a model that represents a specific property
(that you are interested in) better (and be worse at others). judging
simulations only by whether they have a 1:1 match with experimental
data is seriously underestimating its power of predicting trends and
explaining processes at an atomic and molecular level.

2) LAMMPS is an MD code and not a specific parameterization. so if you
*do* want to use a different set of parameters, you are most welcome
to do so. LAMMPS ships potential files as a matter of convenience. the
bundled potential files have been contributed and validated for
specific purposes by LAMMPS contributors. in the special case of
AIREBO, this has been done by the group of steve stuart. if you want
to use a different parameterization, feel free to contact the
corresponding authors or enter (and validate!) the parameters yourself
(and don't forget to contribute the result). parameterization and
validation of potential parameter data is one of the most tedious and
thankless tasks in the simulation business and i wish more people
would show more respect for people undertaking these tasks which are
so incredibly valuable to our community and give a helping hand rather
than complaining and demanding.

axel.

Anees,

No, this has nothing to do with the older/newer versions of parameter sets. Also please note:

  1. The paper re-parameterized REBO, not AIREBO. The paper also did not discuss elastic constants.

  2. It would be difficult for pure REBO to describe graphite well since it lacks the LJ repulsion required to for a stable separation between the graphene layers. Graphite with pure REBO transforms to diamond upon compression far too easily.

I have verified what Narasimha described and obtained similar results as Narasimha did. But I could not figure out why adding a “fix setforce” after the first minimization would result in “correct” answers.

Ray

Yes, I did try different size systems (as big as 100 A). I get the same
results, nothing different.

hmmm... this pretty much reaches the boundaries of my knowledge on this subject.
hopefully somebody with more experience in this specific area can chime in.

sorry,
    axel.

Thanks for your responses, Axel, Anees and Ray.

Ray,

This is what I understand is happening, please correct me if I am wrong.

The elastic example is supposed to determine C11 by changing box dimension in the x dimension alone. Hence strain in x direction is prescribed while in the y and z directions it is kept zero. Getting the stress in the x direction and dividing it with the prescribed strain gives C11.

What seems to be happening here with graphite though is that even as “Ly” remains the same, the atoms seem to be moving around in the y-direction during minimization. What this is doing is creating a configuration wherein the strain in the y-direction is zero “globally” but not so “locally”. When I use the “fix setforce” to constrain the atoms in the y-direction, it artificially makes sure that the strain remains uniformly zero in the y-direction giving me the result I am expecting (which is about 1150 GPa from the original AIREBO paper as I pointed out in my first e-mail).

Hope that made sense.

Thanks,

Narasimha.

The elastic example is supposed to determine C11 by changing box
dimension in the x dimension alone. Hence strain in x direction is
prescribed while in the y and z directions it is kept zero. Getting
the stress in the x direction and dividing it with the prescribed
tstrain gives C11.

Comment : Yup the elastic constant can be defined by the derivative of
the stress vs strain. For C11 then sigma_xx / epsilon_xx , you must
understand the calculation of the elastic constants using this method
is very sensitive to the deformation epsilon , (just try chanding
epsilon to a value of 1.5 and you will see the difference\) soo its recommended to use small displacements 0\.05 up tp 0\.5 , This method
can be considerred computatiobal expensive because you required up to
6 deformations to find all the elastic constants. There are other
METHODS in the liteherature which use the thermal fluctuation of the
computational box or the stress ...Right now i dont remember the
reference sources

What seems to be happening here with graphite though is that even as
"Ly" remains the same, the atoms seem to be moving around in the
y-direction during minimization. What this is doing is creating a
configuration wherein the strain in the y-direction is zero "globally"
but not so "locally". When I use the "fix setforce" to constrain the
atoms in the y-direction, it artificially makes sure that the strain
remains uniformly zero in the y-direction giving me the result I am
expecting (which is about 1150 GPa from the original AIREBO paper as I
pointed out in my first e-mail).

Comment : Am , i dont either know what i am saying, but the condition
for the calculation of the elastic constants using deformations is
that the cell must already wall optimized i.e the Forces must be as
close as possible to zero …If you use setforce to set the force to
zero then may solve the problem, but i belive your cell must be
further optimized and try to check that the forces after the
minimization are close to zero (at least 10^-4) …

A Salute
Oscar G.

If one tries with the LCBOP potential for Carbon provided with the distro the same happens and the numbers compare quite OK.
Carlos