Issue with LAMMPS parameters in bilayer graphene configuration

Dear all LAMMPS Users,

I am writing to share an issue that I have encountered while configuring a bilayer graphene system in LAMMPS, and to request your insights and suggestions.

As you may know, bilayer graphene has a distance of 3.34Å between two layers. However, when I configured my system and minimized it in LAMMPS, I found that the distance between the two layers was about 4.57Å. I suspect that this discrepancy may be due to an error in the parameters that I am using, but I have not been able to identify the source of the problem. I uploaded my reference below this post. I modified the parameter to use proper value. Because in LAMMPS, pair_style lj/cut command used 12-6 LJ potential, but the paper used 6-12 type.

I am hoping that some of you may have encountered similar issues or have insights into what could be causing this discrepancy. If you have any suggestions or recommendations on how I can resolve this issue, I would be very grateful.

Thank you in advance for your help and insights.

LJ potential para.pdf (1.3 MB)

This is my parameter modified version, type 1 is carbon atom, type 2 is Si atom.

pair_coeff 1 1 0.105 4.622601348
pair_coeff 1 2 0.205450724 4.571787923 16.0012577305
pair_coeff 2 2 0.402 4.820974497

bond_coeff 1 469 1.4

angle_coeff 1 63 120

dihedral_coeff 1 0 7.25 0 0

improper_coeff 1 5 180

Additionally, this is minimization part of my LAMMPS inputfile,

---------- Initialization ----------

units real
atom_style full
pair_style lj/cut 14
boundary p p p

bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style harmonic

special_bonds lj 0.0 0.0 0.5

---------- System definition ----------

read_data /home/jungwan/lj_project/bilayer/

include /home/jungwan/lj_project/parameters/PARM_Si_graphene.lammps
region supported_graphene block INF INF INF INF -0.5 0.5
region sample_graphene block INF INF INF INF 2.85 3.85

group substrate region supported_graphene
group sample region sample_graphene

group Si type 2

---------- System Relax ----------

fix freeze substrate setforce 0.0 0.0 0.0
dump 1 all atom 1 minimized.lammpstrj

minimize 1.0e-10 1.0e-10 200000 200000
write_data pair ij

That is just a nomenclature issue. A Lennard-Jones potential where the repulsive term has a smaller power than the attractive term makes no sense. If you look at the sign in formula 20, you can see that the terms are swapped around. A more serious issue is the factor of 2 for the attractive term. This means, that this paper uses a different convention for sigma than what LAMMPS uses. LAMMPS expects sigma to be the distance where the potential energy is zero. When the r^6 term has a prefactor of 2, then this means it assumes sigma to be the minimum of the potential, i.e. the sigma values are too large for use in LAMMPS by a factor of 2^{\frac{1}{6}}. This could be a contributing factor to getting too large a distance.

You are using a “generic” force field. The more generic a force field is, the less accurate it is. This may give OK bonds, but I doubt it will be good for the interlayer distances.

I suggest you have a look at the INTERLAYER package in LAMMPS and check out the various provided styles and their corresponding publications for accurate representation of layered materials.

If you prefer a “generic” force field, I would look at Amber (or perhaps CHARMM) and use the parameters for benzene carbon atoms. Those are likely much more accurate since they are parameterized not for a generic carbon atom, but specifically for carbon in an aromatic 6-ring.

As you can see from the picture, the distance coordinate of the minimum point of the 6-12 type matches LAMMPS’ expected sigma value. Therefore, I did not adjust the value as stated in the paper and ran the simulation. However, I obtained a distance of 3.78-4Å instead of the expected 3.34Å. Is it too large yet to use?

This is my parameter modified version.

pair_coeff 1 1 0.105 3.851
pair_coeff 1 2 0.2535 4.066945414 14.23430895
pair_coeff 2 2 0.402 4.295

I’m sorry, but I didn’t quite understand what you meant by ‘intent.’ Could you please provide me with more information or clarify your request?


When I configured only bilayer graphene with pair_coeff 1 1 0.105 3.851 and minimized the system, I obtained this result. I’d like to understand the cause of this result.

Before minimization

스크린샷 2023-02-22 16-20-55

After minimization

No, it does not. You have just proven my point. For the red curve the potential energy U is -1\epsilon at r = 1.0\sigma, but LAMMPS expects the blue curve with U = 0. That proves that you must divide all \sigma values from the paper by 2^{\frac{1}{6}}.

I think you need to discuss with your adviser/tutor. We are discussing the use of LAMMPS here, not how to interpret the science of your system setup and results of your simulations.

I appreciate your help more than words can express. Thank you for pointing out that I missed that point. Your help has been so valuable to me. I’m looking forward to delving deeper and exploring more complex configurations with topotools. Thank you for your guidance and support.