How to output the van der waals energy in graphene with AIREBO potential

Dear lammps users.

Now I am trying to get the van der waals energy in a double-layer graphene. I can get a value of “evdwl” by modifing the thermo-style as “thermo_style custom step temp etotal pe ke ebond evdwl pzz”, but I found that the the values for etotal, pe, evdwl are the same while the value of ebond is zero. this is kind of weird since AIREBO consists of REBO terms which is bond-order potential, how could I get a zero ebond and a large evdwl?

I also checked the “pair-airebo.cpp” file and found that the evdwl is still calculated in the REBO part. it seems that AIREBO calculate REBO energy into evdwl?

does any body know how to obtain the accurate van der waals energy in double-layer graphene?

Best

Ming

my code is simplified as follows:

atom_style atomic
units metal
boundary f f f

read_data data.dat
mass * 12

lattice sc 1
region 3 block INF INF -1.5 1.5 INF INF #central lines of graphene

#initial the velocity environment and ignore the thermal effect

velocity all create 0.0001 234567

#define the used potential

pair_style airebo 3.0 1 0
pair_coeff * * CH.airebo C C C C C

neighbor 3 bin
timestep 0.001

fix 1 all nvt temp 0.001 0.001 0.001 drag 0.2

#define output style

#trajectory
dump 1 all atom 10000 cnt.lammpstrj
dump_modify 1 format “%7d %3d %10.5e %10.5e %10.5e” scale yes
thermo_style custom step temp etotal pe ke ebond evdwl pzz
thermo 1000

run 500000000

Dear lammps users.

Now I am trying to get the van der waals energy in a double-layer graphene. I can get a value of “evdwl” by modifing the thermo-style as "thermo_style custom step temp etotal pe ke ebond

“evdwl” in LAMMPS includes all charge-independent short range potential energy calculated by the pair_style. For AIREBO this includes everything since AIREBO is essentially a pure short range, charge-independent potential, such that pe is always the same as evdwl.

evdwl pzz", but I found that the the values for etotal, pe, evdwl are the same while the value of

etotal is only equal to pe when there is no kinetic energy, i.e. T is zero.

ebond is zero. this is kind of weird since AIREBO consists of REBO terms which is bond-order

AIREBO is an atomic style potential, i.e., it contains only atoms, not bonds. ebond is only not zero for bond_style potentials.

potential, how could I get a zero ebond and a large evdwl?

AIREBO, being a pair_style, has no true bonds; therefore ebond must equal to zero.

I also checked the “pair-airebo.cpp” file and found that the evdwl is still calculated in the REBO part. it seems that AIREBO calculate REBO energy into evdwl?

does any body know how to obtain the accurate van der waals energy in double-layer graphene?

You can’t obtain van der waals energy from AIREBO on-the-fly. You may use the rerun command, with the LJ_flag on and off to get the difference.

Ray

Dear Ray and lammps-users,

Thx very much for your comments, which help me a lot.

in your letter, you suggest that I can obtain the van der waals energy by turning on and off the LJ option and get the difference. actually I have tried this method, but the result is not right.

for a circular double-layer graphene with radius 50 angstrom and 3007 atoms in each layer, we assume E1 as the total energy of the double-layer graphene when LJ flag is on while E2 is the total energy of the double-layer graphene when LJ option is off. in my simulations, E1=46168.206 eV, E2=45978.228 eV, the per-atom van der Waals energy can be defined as (E1-E2)/3007=0.063 eV.

also we try another method and assume E3 as the total energy for the single-layer graphene, the LJ energy for the double-layer graphene can be defined as (E2-2E3)/atom_num. since E3 for the single-layer graphene is 23028.788 eV, the per-atom LJ energy can be calculated as (46168.206-223028.788)/3007=0.036791, this value agree well with the reference paper (Figure 5a in Nanotechnology 18 (2007) 395703).

I attached my in and data files here.the data file for the single-layer graphene can be obtained by removing the atoms beyond 3007.

thus, by turning on and off the LJ option in the AIREBO potential, we can not get the right van der waals energy. but the second method is time consuming, I want a simple way to get the accurate van der waals energy when using the AIREBO potential.

any comments?

thx very much,

Best

Ming

data.dat (365 KB)

in.cnt (895 Bytes)

Well, I have to say you are comparing apples to oranges. The reference you provided does not even use AIREBO potential, and the LJ parameters used were totally different from AIREBO published values. How can LJ energy be the same with different parameters?

From single point calculations, the first and second methods give the exact same energy difference - this is only because the two layers are ideally 3.4 angstroms apart.

Please note that the first method gives you the LJ energy, while the second, depending on inter-layer separation, may give you something else. You can not rule out other energy terms if the layers are too close, which is what might have happened after 500 million steps in your original script (which I have changed to 1 of course).

My suggestion is to read more MD and AIREBO papers.

Ray