Comparison of AIREBO and Tersoff potentials in water–graphene simulations and ambiguity in using Tersoff

[Declan.Harris]

Hello all,

I have implemented the AIREBO potential in my simulations and obtained good agreement with a reference paper regarding Kapitza resistance and heat flux. However, **I’m facing issues when switching to the Tersoff potential: graphene layers deform significantly (bend or warp) during the simulation.Snapshots taken toward the end of both the AIREBO and Tersoff simulations are attached for visual comparison.

The simulation setup follows the LAMMPS manual, relevant discussions on this forum, and Professor Axel Kohlmeyer’s explanations.

Simulation setup:

10 graphene layers on each side, ~4 nm water in between.

Each graphene layer is defined as a separate atom type (only in the Tersoff simulation).

Hot and cold thermostats are applied to the outermost graphene layers to induce a heat flux for Kapitza resistance calculation.

The outermost graphene layers on the left and right sides are fixed in place using fix setforce.

Potential definitions in Tersoff simulation (example for 2 layers):

Blockquote

pair_style hybrid tersoff tersoff lj/cut 10.2 lj/cut/coul/long 12.0

Carbon intralayer

pair_coeff * * tersoff 1 SiC_1994.tersoff NULL NULL C NULL
pair_coeff * * tersoff 2 SiC_1994.tersoff NULL NULL NULL C

Carbon interlayer

pair_coeff 3 4 lj/cut 0.00239721633 3.35

Water-water

pair_coeff 1 1 lj/cut/coul/long 0.0 0.0
pair_coeff 1 2 lj/cut/coul/long 0.0 0.0
pair_coeff 2 2 lj/cut/coul/long 0.00673977503 3.165

Carbon–H

pair_coeff 1 3 lj/cut/coul/long 0.0 0.0
pair_coeff 1 4 lj/cut/coul/long 0.0 0.0

Carbon–O

pair_coeff 2 3 lj/cut/coul/long 0.0040627276 3.190
pair_coeff 2 4 lj/cut/coul/long 0.0040627276 3.190

Blockquote

This pattern is extended accordingly for each of the 20 graphene layers in the full simulation.

To prevent deformation of the graphene sheets, I have separately tested both fix spring and fix spring/self. The reference paper suggests fixing the center of mass (COM) of each graphene layer, so I used fix spring for this purpose. However, it didn’t prevent warping.

Alternatively, I tried fix spring/self, which ties each atom to its initial position. This reduced deformation and matched AIREBO’s behavior more closely. But since this approach constrains individual atoms rather than the COM, it seems unphysical and also alters the heat flux depending on the spring constant.
Also, in the Tersoff simulation without using fix spring/self, the Kapitza resistance is about 13% lower than the expected (reference) value.

Is this kind of deformation expected behavior from the Tersoff potential itself?
If not, could it be due to an issue in the way Tersoff interactions are defined in my input script?
If the Tersoff part is correctly implemented, is the use of fix spring/self a recommended workaround in this context?
What would be the best way to make the behavior of graphene in the Tersoff simulation more similar to the AIREBO case?

In summary, with the Tersoff potential I do not reach the expected results either physically or numerically. I am searching for the cause of this issue—whether it is due to the Tersoff coding itself or if other modifications, such as adding auxiliary commands like fix spring/self, are needed to obtain correct behavior.

I have searched extensively through the LAMMPS manual, related papers, forum threads, and every other source I could think of, but I have not yet found a definitive answer to this ambiguity.
Any help would be highly appreciated. Even a small hint or a keyword to further guide my search would be extremely helpful and valuable.

[akohlmey]

@Declan.Harris Please note that what you are asking is effectively that somebody here explains your research to you. That is not the purpose of this forum, but rather something to discuss with the people that advise and tutor you.

@Declan.Harris

Is this kind of deformation expected behavior from the Tersoff potential itself?
If not, could it be due to an issue in the way Tersoff interactions are defined in my input script?

I cannot answer this in general and with great confidence, since there is simply too little tangible information in your post, but I can speculate: have you considered that AIREBO and your Tersoff parameterization may have a different equilibrium C-C distance? Your snapshot image looks like the the sheets are buckling because they cannot relax the way they “want” and “be straight” due to too small box lengths in the periodic dimensions.

@Declan.Harris

I have searched extensively through the LAMMPS manual, related papers, forum threads

I don’t think that this is a LAMMPS issue per se and thus you cannot expect to find a direct answer. Instead, you will have to apply ‘the scientific method’ and try to reason this out. But first and foremost, you need talk to people that care about your results. They may then want to reach out to other experts for additional suggestions.

@Declan.Harris

adding auxiliary commands like fix spring/self, are needed to obtain correct behavior.

This makes no sense. It would just take your simulation and make it a computer animation. Again, please note that this is outside my area of expertise, but note that in LAMMPS there is the INTERLAYER package with potentials for the expressed purpose to improve Tersoff and alike for multi-layer configurations. What kind of impact those would have is something that you need to learn from the corresponding publications. I don’t think that can solve the issue of buckling sheets.

@Declan.Harris

but I have not yet found a definitive answer to this ambiguity

If there were simple and easy to find answers, then this would not be called research but school.

[Declan.Harris]

Thank you for the comments you shared. I will follow up on the points you mentioned. You are right, and I am aware of the research-based nature of the work. This topic was not about the nature of the research itself, but rather for consultation and guidance regarding the correct implementation of the Tersoff potential in the simulation, and related aspects that may have been overlooked by me.

I will continue searching into the correct implementation of the Tersoff potential and will keep your guidance in mind. Once again, thank you for the points you shared and I appreciate them.

[akohlmey]

@Declan.Harris

but rather for consultation and guidance regarding the correct implementation of the Tersoff potential in the simulation

Sorry, but as far as LAMMPS is concerned, it only matters whether you are using the correct syntax. Everything else is not a LAMMPS issue, but a research issue. LAMMPS does not care whether you feed it random numbers, it will take them apply the algorithms that are implemented and then show the outcome. To decide whether that is meaningful, is entirely up to people that are experts in that particular area of research and also independent of the simulation software.