Dear LAMMPS Users&Developers,
I would like to ask your help in the following problem.
I would like to calculate the surface energy of copper (1,1,1).
For this purpose I created 2 LAMMPS scripts to determine the energy of bulk copper and copper slab.
Here is the input file for bulk:
units real
dimension 3
boundary p p p
atom_style charge
lattice fcc 3.615
region total_box block -22 -6 -22 -6 -4 4 units lattice
create_box 1 total_box
create_atoms 1 region total_box
mass 1 63.546
pair_style reaxff lmp_control
pair_coeff * * ffield.reax.cuo Cu
fix 1 all qeq/reax 1 0.0 10.0 1e-6 param.qeq
neighbor 2.0 bin
neigh_modify every 10 delay 0 check yes
thermo 100
thermo_style custom step pe lx ly lz etotal
min_style cg
minimize 1.0e-4 1.0e-6 8000 10000
write_data data.bulk_minimized_v02
And the file for copper slab:
units real
dimension 3
boundary p p p
atom_style charge
lattice fcc 3.615 orient x 1 -1 0 orient y 1 1 -2 orient z 1 1 1
region total_box block -22 -6 -22 -6 -34 34 units lattice
create_box 1 total_box
mass 1 63.546
region slab1 block -22 -6 -22 -6 -4 4 units lattice
create_atoms 1 region slab1
pair_style reaxff lmp_control
pair_coeff * * ffield.reax.cuo Cu
fix 1 all qeq/reax 1 0.0 10.0 1e-6 param.qeq
neighbor 2.0 bin
neigh_modify every 10 delay 0 check yes
thermo 100
thermo_style custom step pe lx ly lz etotal
min_style cg
minimize 1.0e-4 1.0e-6 8000 10000
write_data data.slab_minimized_v02
In case of slab I track lx ly lz to determine the surface.
After that I use gamma=(E_slab-N_slabE_bulk)/(2A) correlation.
I got the result 0.124 eV/A^2 but the real value is far from it.
What did I do wrong? Is there any issue in the slab simulation or already in the concept?
Thank you in advance for your help!
Patrik
Have you discussed this with your adviser/tutor already?
What is a “real value”?
What is a “real value”?
Yes, sorry for wrong designation. I mean the calculated value in the literatures 0.08-0.09 eV/A^2.
My result(s) is (are) far from that, around 0.15 eV/A^2.
Hi! Following up on this as I am looking at something similar? Did you ever figure it out? Your accuracy seems reasonable when comparing to experiment in my opinion.
Also, why do you multiply N_slab*E_bulk. I would think you would first need E_bulk as per atom before multiplying? I think my script is having an error at this stage
So my code looks like below, I simplified to just a Na case to see what was going on. There must be an issue with my method as the minimized energies are always the same. I take a periodic slab first, minimize and calculate pe. I think open this up to vacuum on both sides, mninimze again, and get pe again. But the number is the same?
# Initialization -------------------------------
units real
dimension 3
processors * * *
boundary p p p
atom_style full
box tilt large
atom_modify map array
read_data crystalNa.data
neighbor 2 bin
neigh_modify every 1 check yes
# reaxff force field
pair_style reaxff NULL checkqeq yes
pair_coeff * * albitereax.ff.txt Na
fix 21 all qeq/reax 1 0.0 10.0 1.0e-6 reaxff
timestep 0.5
#Data variable setup
variable output_data string "output.txt" #Text file
# ---- dimension variables ----
# box length in x and y
variable xlo equal "xlo"
variable xhi equal "xhi"
variable ylo equal "ylo"
variable yhi equal "yhi"
# height of layer
variable zlo equal "zlo"
variable zhi equal "zhi"
variable albite_temp equal 1 #kelvin
variable mass_Al equal 26.98
variable mass_Si equal 28.09
variable mass_Na equal 22.99
variable mass_O equal 16
mass 1 ${mass_Na}
# BULK -------------------------------
thermo_style custom atoms step dt time temp etotal pe
min_style cg
minimize 1.0e-8 1.0e-10 8000 10000
write_data bulkNa.data
#Print to output text file.
variable bulkte equal "etotal"
variable bulkpe equal "pe"
variable N_bulk equal count(all)
print ${bulkpe}
shell echo bulk pe is ${bulkpe}. >> ${output_data}
variable E_bulk_per_atom equal ${bulkpe}/${N_bulk}
# slab -------------------------------
#make surface
change_box all z final -20 100
change_box all boundary p p f
min_style cg
minimize 1.0e-8 1.0e-10 8000 10000
write_date surfaceNa.data
#Print to output text file.
variable slabte equal "etotal"
variable slabpe equal "pe"
variable N_surf equal "atoms"
variable E_bulk_for_slab equal ${E_bulk_per_atom}*${N_surf}
shell echo slab pe is ${slabpe}. >> ${output_data}
# ---- SFE ----
# Calculate surface energy per unit area (A = surface area of slab)
variable A equal "lx*ly" # Area of the surface
print "Areas is ${A} slab pe is ${slabpe} and bulk pe is print ${bulkpe}"
variable gamma equal ((0.043*(${slabpe}-${E_bulk_for_slab}))/2*${A})
print ${gamma}
shell echo SFE is ${gamma}. >> ${output_data}