NPT, NVT Graphene-SiO2

Dear Lammps users,

Hi, I am calculating thermal conductivity of Graphene on amorphous SiO2 substrate using both EMD (Green-Kubo) and NEMD (MP) methods. I created an initial structure set up with VMD topotools (nanotube & inorganic material builder), which is contained in a data file. I am using airebo for C, tersoff for Si and O, and vdW (lj/cut) for interactions between C-Si and C-O.

For both EMD and NEMD simulations, problems occurred during the structure relaxation. During NPT or NVT run, the graphene sheet on the amorphous SiO2 moved far away from the substrate, eventually reached the box top, and broke. I have set the boundary condition to be (p p f) where graphene is on xy plane. I first though it might have been the vdW repulsion-attraction phenomena and checked the parameters for lj/cut interactions, but they seemed to be correct (parameter values extracted from literature). I wonder if anyone has any clue on why this is happening.

One more question. It does not seem to be very meaningful to relax amorphous structure though, I just tried to relax both Graphene and SiO2. When I did that, I observed some atoms were lost (or moved out of the box). I set the simulation to ignore the lost atoms, so the simulation could run nonetheless though, is there anyway to prevent the atoms to move that much?

Please find the input script below.

Thank you in advance,

Chris

setup problem

units metal
boundary p p f
atom_style atomic

read_data …/data.mergedGrSiO2
lattice custom ${alat} a1 1 0 0 a2 0 1.73203 0 a3 0 0 10 basis 0 0 0 basis 0.5 0.16666666666666666 0 basis 0.5 0.5 0 basis 0 0.6666666666666666 0
group gr id <= 2304
group sub id > 2304

pair_style hybrid airebo 3.0 tersoff lj/cut 2.5
pair_coeff * * airebo /home/kpark39/programs/lammps-1Feb14/potentials/CH.airebo C NULL NULL
pair_coeff * * tersoff /home/kpark39/programs/lammps-1Feb14/potentials/SiO.tersoff NULL O Si
pair_coeff 1 2 lj/cut 0.006734 3.275
pair_coeff 1 3 lj/cut 0.008890 3.628

neighbor 2.0 bin
neigh_modify delay 10
timestep ${dt}

thermo_modify lost ignore
velocity all create $T 4928459 mom yes rot yes dist gaussian

dump 1 all cfg 10000 dump.config.*.cfg mass type xs ys zs vx vy vz x y z

1st equilibration run

#fix 1 all nvt temp $T $T 0.5
fix 1 all npt temp $T $T 0.5 x 0 0 0.5 y 0 0 0.5

thermo 1000
run 50000

velocity all scale $T

unfix 1

undump 1

thermal conductivity calculation

reset_timestep 0

compute myKE all ke/atom
compute myPE all pe/atom
compute myStress all stress/atom virial
compute flux all heat/flux myKE myPE myStress
variable Jx equal c_flux[1]/vol
variable Jy equal c_flux[2]/vol
variable Jz equal c_flux[3]/vol

dump 2 all cfg $d dump.config..cfg mass type xs ys zs vx vy vz x y z
fix 1 all nve
fix JJ all ave/correlate $s $p d & c_flux[1] c_flux[2] c_flux[3] type auto & file profile.heatflux ave running variable V2 equal 80*80*60 variable scale equal {convert}/${kB}
s*{dt}/$T/T/{V2}
variable k11 equal trap(f_JJ[3]){scale} variable k22 equal trap(f_JJ[4])*{scale}
variable k33 equal trap(f_JJ[5])
${scale}

thermo $d
thermo_style custom step temp v_Jx v_Jy v_Jz v_k11 v_k22 v_k33

run 2000000

Dear Lammps users,

Hi, I am calculating thermal conductivity of Graphene on amorphous SiO2
substrate using both EMD (Green-Kubo) and NEMD (MP) methods. I created an
initial structure set up with VMD topotools (nanotube & inorganic material
builder), which is contained in a data file. I am using airebo for C,
tersoff for Si and O, and vdW (lj/cut) for interactions between C-Si and
C-O.

For both EMD and NEMD simulations, problems occurred during the structure
relaxation. During NPT or NVT run, the graphene sheet on the amorphous SiO2
moved far away from the substrate, eventually reached the box top, and
broke. I have set the boundary condition to be (p p f) where graphene is on
xy plane. I first though it might have been the vdW repulsion-attraction
phenomena and checked the parameters for lj/cut interactions, but they
seemed to be correct (parameter values extracted from literature). I wonder
if anyone has any clue on why this is happening.

​most likely this is due to unequilibrated initial structures​ with high
potential energy and possibly other issues.

One more question. It does not seem to be very meaningful to relax
amorphous structure though, I just tried to relax both Graphene and SiO2.
When I did that, I observed some atoms were lost (or moved out of the box).
I set the simulation to ignore the lost atoms, so the simulation could run
nonetheless though, is there anyway to prevent the atoms to move that much?

​there are multiple issues with your overall approach.

- equilibration is *very* much needed, also for amorphous structure.
whatever typical structure building tools produce, is often not very
refined and rarely in keeping with the specific potential you choose.

- amorphous structures are often created by melting some starting
structures and then cooling it down and relaxing it with a given cooling
rate. please check the relevant literature

- for SiO2 you have make sure you have no dangling bonds, you may also need
to consider whether you have outward facing Si-OH on your surface instead
of Si-O.

- have you checked, whether the Si-O tersoff potential was validated for
modeling surfaces?

- you should only ignore lost atoms, when you *want* atoms to leave the
simulation box. this is clearly not the case here. thus lost atoms are a
sign of a bad choice of simulation parameters or a bad system setup.

- for slab systems as yours, you should immobilize the lowest layer(s) of
the slab to maintain a bulk-like structure. that will also avoid a drift of
the whole system (and thus lost atoms).

- it is probably best to first build and relax each part of your simulation
separately and verify that you get a correct simulation than can represent
your system well for each.

- before even trying advanced simulation techniques, that require practice
and experience like thermal conductivity studies or a multi-component
system, you should practice getting simpler simulations going and reproduce
published data for such simpler system. there is no point in trying
something more complex without knowing how to get it done right (and being
able to tell from having messed it up and corrected before).

it overall appears to me that you should find yourself a tutor/instructor
that can teach you on the fundamentals of MD simulations and best practices
in general and especially for the kinds of systems that you are interested
in. "flying blind" like you currently seem to be doing is extremely risky
and has a very high chance to produce meaningless data and waste a lot of
your time in the process.

axel.