Dear Bo,
I have not observed this behavior myself. Did you start off with the correct stacking for the two graphite planes? Here is a link to the correct structure for graphite: http://cst-www.nrl.navy.mil/lattice/struk/a9.html If you started with all the carbons in each plane on top of one another, one layer should shift with respect to the other such that the proper structure is obtained.
Good luck,
Dave Schall
Dear Dave,
Thanks for the link! I created the graphite structure from the unit cell looks like below:
two layers in the unit cell in c direction, one at fractional 0.25, the other at 0.75
one layer looks like:
* *
Hi Bo,
Your diagram looks correct. I just double-checked the Stuart Airebo paper, looks like he’s using 3 sigma to reproduce g® data for molecule systems. I imagine they used that number for fitting the graphite data as well. When you do your calculation double check that you get the correct interplanar spacing of 3.354 Angstroms. I’m guessing your number will be slightly different, because you said you were only doing 2 layers.
If you are doing two layers with free surfaces top and bottom, I imagine you would need ppf for your boundaries or at least you’d want to have plenty of room (rule of thumb: >3 sigma) above and below the surface. I haven’t done this calculation with Lammps, but I was a grad student with Brenner and postdoc with Harrison, two of the coauthors on the airebo paper and used their research code extensively. It’s been awhile but I am almost certain I got the correct structure. Way back I wrote a paper on sliding and rolling of nanotubes on graphite where I had multiple layers of graphite.
Let me know what you find.
Dave
Dear Dave,
I've been doing various tests, however, everyone fails:
1. ppp, graphite, airebo 4.0 1 0, slide of layers
2. ppf, graphene 2 ly, airebo 1.0, atoms lost within short time
3. ppf, graphene 2 ly, airebo 4.0, slide of layers
4. ppp, graphite, airebo 3.0, tstep reduced from 0.0005 to 0.0001,slide of layers
To your question, I did NPT at the beginning, but only relaxes the xy dimensions while z dimension is fixed, so yes, the interlayer separation is about 3.34 A, and the measured inplane bonds are about 1.42 A. Just don't know why the sliding happens, so frustrating.
below is the input file I have for test 4, attached is the xyz trajectory from test 4:
units metal
dimension 3
boundary p p p
atom_style atomic
pair_style airebo 3.0
read_data lammps.txt
pair_coeff * * CH.airebo C
velocity all create 300 102486 mom yes rot yes dist gaussian
#compute 1 all temp
#velocity all create 300 873586443 temp 1
timestep 0.0001
thermo 1000
# ------------- Equilibration and thermalisation ----------------
dump 1 all xyz 100000 aa.xyz
fix NPT all npt temp 300 300 0.5 x 0.0 0.0 10.0 y 0.0 0.0 1.5 drag 0.2
run 1000000
unfix NPT
# --------------- Equilibration in nve -----------------
#dump 1 all xyz 100 aa.xyz
fix NVE all nve
run 21000000
Please help me take a look, thanks a lot!
Best,
Bo
Dave Schall wrote:
aa.xyz (935 KB)
Also attached is the input geometry I have for graphite. Thanks!
Bo
Bo Qiu wrote:
lammps.txt (11.9 KB)
Bo,
I am not seeing any sliding. I changed only one thing. Starting from case 4, I changed your zhi zlo in your coord file to a number greater than 3 sigma. Look in the airebo paper for the values of sigma if you need to. Anytime your system size is smaller than your cutoff for the interatomic (or in this case intermolecular) potential atoms in your system can see themselves across the period boundaries. Obviously that is a big problem and will usually lead to your system exploding or a stacking dump or some other nastiness. I suspect that was your problem.
Also, is there any reason why you need to run for 21e6 timesteps? That’s a really long time. There is a remote chance that in that amount of time one layer could diffuse relative to the other. Probably not very likely. Anyway, I ran for 10K steps in NPT and 100K steps in NVE (with the changes above) and saw nothing to indicate any sliding. I also output xyz files every 100 steps an set thermo = 10. Looked like energy was nicely conserved in NVE. Seemed good to me.
Give it a try and let me know how it goes.
Dave
Dear Dave,
Your suggestions are super helpful, thanks a lot! Shame I should have thought about it earlier.
Though it should work for graphite with large z of the box, I'm wondering how that can help remove the sliding problem in double graphene sheets case, since the boundary is ppf, the z box size should not matter.
I'm trying to calculate some correlation function which generally requires several ns data, corresponding to 2000K data length with 0.5fs timestep. I set 21e6 steps in the input file only to see whether the system is stable in sufficiently long time.
I'll surely give it a try and keep you posted.
Thanks again for your great help!
Bo
Dave Schall wrote:
Dear Bo and Dave
i hope its alright if i ask something ontop of your questions 
i am a new comer to lammps and MD in general.
when dave said that “energy was nicely conserved in NVE” how much of a loss in energy is acceptable
i am asking this since i ran the simulation script Bo attached to see and it seems that the energy fluctuates around and interval of give or take 4.30 i know that there is an error while doing numerical integration with verlet algorithm but is this acceptable ?
i am asking since i have also ran a very simple script of a harmonic chain and there also the energy fluctuates strongly but it should be a constant ,any and all explanations will be helpful of what am i doing wrong
i attach the input script and data as well since i have another problem that when i try to run it with MPI i get an error of bonds missing on prcs .
thanks a lot in advance !! and i apologize for piggyback riding on your posts
Tal
the simulation:
#simulation for harmonic chain
units lj
dimension 2
boundary p p p
atom_style bond
read_data bonds.txt
mass 1 1.0
velocity all create 0.1 4928459 rot no dist gaussian
bond_style harmonic
bond_coeff 1 1 16
neighbor 0.3 nsq
neigh_modify delay 5
atom_modify sort 1 1
fix NVE all nve
fix Z all enforce2d
fix 1Dy all lineforce 1 0 0
dump 1 all xyz 100 data.xyz
thermo 250
thermo_style custom pe ke emol etotal
run 10000
Bo Qiu wrote:
bonds.txt (2.66 KB)
You still need a large z dimension with ppf. In your initial set up, the z-coordinate of the atoms in the graphene sheet were close to zlo. This is why you were loosing atoms towards the beginning of your run. Even with large zlo and zhi, eventually your system may drift past zlo and be lost, especially during long run times. To keep this from happening you would have to subtract the COM from the motion of the atoms adding a little additional computational cost. The best solution is to keep the system in ppp and set zhi zlo large enough for the periodic image to not see itself (> 3 sigma) otherwise there is no reason why the atoms can’t escape from the box to be lost forever.
Dear Dave,
I'm back with the results:
the system momentum and angular momentum are all zeroed.
test #1:
ppp graphite with 3 unit cells in z direction (1200 atoms total),
0.0 25.29298 xlo xhi
0.0 21.90436 ylo yhi
0.0 20.112 zlo zhi
in input file:
pair_style airebo 2.8
timestep 0.0005
so the cutoff is not larger than half of the box size in z direction.
results: shortly after about 500K steps, all layers slide randomly in x-y plane, while their z coordinates are not changed much. Attached is some snapshots of the trajectory.
test #2:
ppf 2-layer graphene (400 atoms total) with box size in z enlarged:
-6.0 12.704 zlo zhi
so the cutoff is not larger than half of the box size in z direction.
results: after a while the two layers slide randomly in x-y plane, while their z coordinates are not changed much.
By the way, I tried 200K steps on 2-layer graphene with large box and ppp, as you have tested previously. Indeed, I didn't see any sliding. However, when I have longer simulation, the sliding still shows up.
I still feel there are something not right somewhere, but just cannot tell...
Thanks as always,
Bo
Dave Schall wrote:
aa.xyz (359 KB)