[lammps-users] Measuring the Stacking Fault Energy

Hi everybody,

I want to measure the instable stacking fault energy for some potentials. For this I want to use the method described in Vitek (Intrinsic Stacking Faults in Body-centred Cubic Crystals; 1968). The chosen potential is Mishin1 copper from PRB63:224106 (other potentials deliver analogous results). Shortly the method I want to use is:
(1) Shift the crystal in one plane (i.e. [-211] with 0<=dx<1)
(1a) The particles mustn’t move in lateral direction
(2) Calculate E_total(dx)/A_surface as a function of dx
(3) Plot ( E_total(dx) - E_total(dx=0) )/A_surface

The problem is now, that:
(1) The function E_total(dx)/A_surface is not smooth in dx
(2) The concrete values given in Mishins paper cannot be reproduced
(3) There is a strong dependency of the results on the number of units cells

This is basically what I have done:
******* SCRIPT *******
boundary p s p

units metal
atom_style atomic

lattice fcc 3.615 orient x 1 1 -2 orient y 1 1 1 orient z 1 -1 0

variable normal equal 5
variable lateral equal 10
variable llo equal 0
variable lhi equal {normal} variable dy equal div(1,3) variable mlo equal add({lhi},{dy}) variable mhi equal add(add({mlo},{normal}),{dy})
region lower block 0 {lateral} {llo} {lhi} 0 {lateral} units lattice
region middle block 0 {lateral} {mlo} {mhi} 0 {lateral} units lattice
region box union 2 lower middle
create_box 2 box

create_atoms 1 region lower

variable dx equal .5
variable dz equal 0
log log-x${dx}-z${dz}

lattice fcc 3.615 orient x 1 1 -2 orient y 1 1 1 orient z 1 -1 0 origin {dx} .0 {dz}
create_atoms 2 region middle

mass * 63.55
pair_style eam/alloy/opt
pair_coeff * * Cu_Mishin1.eam Cu Cu

thermo 1
variable s0 equal 0
variable sunit equal 16021.765
variable s equal sub(mult(div(etotal,mult(lx,lz)),{sunit}),{s0})
compute s all variable s
compute x all variable dx
compute z all variable dz
thermo_style custom c_x c_z lx ly lz etotal c_s temp step atoms pyy
thermo_modify format float %15.15g

fix notlateralforce all setforce 0 NULL 0

dump out all custom 100 dump.*.bin x y z

min_style cg
minimize 1e-20 1e6 1e6
***** SCRIPT *****

What am I missing? :confused:

Best regards:

I don't see where you're shifting the system or applying a dx != 0.0
All you're doing is relaxing a perfect crystal with freedom in
y, but no movement in x,z ?