Equilibrium lattice parameter

Hello all,

I am trying to find out equilibrium lattice parameter of Cu by applying NPT thermodynamics. However, it keeps on changing in output. With each given input of lattice constant, it is coming out to be different How to find accurately the Equilibrium lattice constant?

Following is the program:

---------- Initialize Simulation ---------------------

units metal
dimension 3
boundary p p p
atom_style atomic
atom_modify map array

---------- Create Atoms ---------------------

lattice fcc 3.631618387
region box block 0 1 0 1 0 1 units lattice
create_box 1 box
lattice fcc 3.631618387 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
create_atoms 1 box
replicate 10 10 10

---------- Define Interatomic Potential ---------------------

pair_style eam
pair_coeff * * Cu_u3.eam
neighbor 2.0 bin
neigh_modify delay 10 check yes

reset_timestep 0

compute new all temp
compute eng all pe/atom
compute eatoms all reduce sum c_eng

velocity all create 300 2223333 temp new
fix 1 all npt temp 300 300 0.2 iso 0.0 0.0 1
###fix 2 all box/relax iso 0.0 vmax 0.001

#------ output the lattice constant -----#

thermo 100
#thermo_style custom step temp lx ly lz
thermo_style custom step temp press pxx pyy pzz lx ly lz c_eatoms
run 200000
variable natoms equal “count(all)”
variable teng equal “c_eatoms”
variable a equal “lx/10”
variable ecoh equal “v_teng/v_natoms”

print “Lattice constant(Angstroms) = {a};" print "No. of atoms = {natoms};”
print “Total energy (eV) = {teng};" print "Cohesive energy (eV) = {ecoh};”

variable length equal “lx/10”
print “Lattice constant (Angstoms) = ${length};”
print “All done!”

Did you try plotting lx, ly and lz as a function of time and seeing if it has converged. If these curves are oscillating about an average value, take the average value and divide it by ten to get the lattice constant.

The oscillations are due to thermal vibrations at 300 K.

Best Regards

Thank you Mr. Warrier.