Recently I want to try a Tersoff + LJ model for graphite and do GK calculations for thermal conductivity. Before I did that, for a single layer graphene with 576 atoms, I tried:
For input #1, I only used 1 type of atoms, with tersoff potential only. For input #2, I defined 2 types of atoms in data file, however, all the atoms are of type 1 (no type 2 atoms in data file). Hybrid tersoff + LJ potential is used so that tersoff only acts for type 1-1, 2-2 while LJ only acts for type 1-2 interactions.
I suppose the results from both inputs should be identical since only tersoff potential is active in single layer graphene by definition. However, for the same random seed and same MD system, I'm getting quite different results in the final thermal conductivity. If I look at the log file, the temperature and energy evolution are also different.
I don't know what part did I get wrong, could you please help? Thanks!
Is the thermo output different immediately?
If so, please post sample files for a small problem
that reproduces it.
Also, this is not possible:
Hybrid tersoff + LJ
potential is used so that tersoff only acts for type 1-1, 2-2 while LJ
only acts for type 1-2 interactions.
Tersoff is a many-body potential. When you define it
for atom types 1 and 2, then all Tersoff many-body interactions involving
any atoms of type 1 or 2 are calculated.
So try removing the pair_coeff 1 2 none statement.
The thermo outputs are different immediately. I have attached the input, data files, and sample output files to this email.
About the Tersoff potential, I'm a little bit confused, you mentioned "then all Tersoff many-body interactions involving any atoms of type 1 or 2 are calculated.", do you mean "type 1 OR 2" or "type 1 AND 2"? If it's "OR", does that mean when I use "pair_coeff 1 2 none", ALL interactions are excluded? If that's the case, there will be no interaction in the system, how could they evolve with time? Thanks!
The thermo outputs are different immediately. I have attached the input,
data files, and sample output files to this email.
actually, they only differ in the second line of output (i.e. after
10,000 steps),
not the first. and if you change the thermo output to be more frequent,
you'll see that they will be identical for quite a while. with NPT integration
and the fact that the pressure tensor is much more sensitive to noise
resulting from truncation errors due to using floating-point numbers,
this kind of divergence after 10,000 steps is to be expected. up to about
4000 steps numbers are pretty comparable.
About the Tersoff potential, I'm a little bit confused, you mentioned "then
all Tersoff many-body interactions involving any atoms of type 1 or 2 are
calculated.", do you mean "type 1 OR 2" or "type 1 AND 2"? If it's "OR",
does that mean when I use "pair_coeff 1 2 none", ALL interactions are
excluded? If that's the case, there will be no interaction in the system,
how could they evolve with time? Thanks!
tersoff is not a pairwise potential. you select the interactions that it
operates on by selecting atom type not by pairs of atom types. this is the
reason why you have to always specify pair_coeff * * for many-body potentials.
if you include types 1 and 2, then interactions 1-2, 2-1, 1-1-1, 1-1-2, 1-2-1,
1-2-2, 2-1-1, 2-1-2, 2-2-1, and 2-2-2 will be handled by tersoff.
Thanks for the clarification, just want to ask one more question:
If I want to study graphite, and I artificially assign different layers with different atom types and say Tersoff only acts within a layer but never stretches outside (even for many body interactions). Do you think the below command lines can do what I hope them to? Thanks!
************** layer 1, type 1
************** layer 2, type 2
pair_style hybrid tersoff
pair_coeff * * tersoff SiC.tersoff C C
pair_coeff 1 2 none
Thanks for the clarification, just want to ask one more question:
If I want to study graphite, and I artificially assign different layers with
different atom types and say Tersoff only acts within a layer but never
stretches outside (even for many body interactions). Do you think the below
command lines can do what I hope them to? Thanks!
no. this would only work for pairwise additive potentials.
and it doesn't make any sense in the first place.
if there are no interactions between the two layers,
you can just simulate one layer. the outcome will
be the same as simulating two layers that don't
"see" each other.
in order to do this (meaningless) calculation, you
would need a use a different pair style for each layer
(which could all compute the same functional form...)