Dear Aidan,
As you may remember, I'm attermpting to carry out simulations of
crystalline ZrO2 using the reaxFF potential. I've been unable to
equilibrate the system at finite temperature, so I checked the
zero-T energy per atom as a function of lattice parameter, which
gives a measure of the interatomic potential in the crystalline
environment. This turns out to have a double minimum (plot column
14 as a function of column 13), which is unphysical. I can think
of several possible causes of this error:
(1) My LAMMPS script (below) contains an error. It seems unlikely,
because it's very simple, but could you please have a look at it?
(2) The results aren't coverged. However, I've checked that they're
independent of the parameters governing hbcut and the charge
equilibration precision, and the double minimum is robust to
variations of the lower and upper limits to fix deform.
(3) The potential I'm using (published as Supplementary Information
to J. Chem Phys 112, 3133 (2008), and attached) contains an error. I've written
to the authors but haven't yet received a reply.
(4) LAMMPS and the reaxFF libraries weren't properly compiled /
linked on my system. Could a problem here cause the observed
behaviour? If so, would it be possible for you to run my script on
your system and compare the results?
Thanks for your help,
Philip
############### Parameters ########################################
######## parameters for the potential ############################
# cutoff distance for hydrogen bonds in Armstrongs
variable hbcut equal 1
# relative accuracy of precision for charge equilibration
variable charge_precision equal 1.0e-6
######## parameters for structure ################################
# number of unit cells in x-, y- and z-directions:
variable size equal 6
# initial lattice constants:
variable ax_i equal 4.8
variable ay_i equal 4.8
variable az_i equal 4.8
# final lattice constants:
variable a_f equal 5.4
# number of steps for deformation:
variable deform_steps equal 41
# final size of simulation cell
variable L_f equal \{a\_f\}\*{size}
############## Initialize ########################################
boundary p p p
units real
atom_style charge
variable relca equal \{az\_i\}/{ax_i}
variable relba equal \{ay\_i\}/{ax_i}
lattice custom \{ax\_i\} a1 1\.0 0\.0 0\.0 a2 0\.0 {relba} 0.0 a3 0.0 0.0 ${relca}&
origin 0.0 0.0 0.0 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.5 0.0 0.5 basis 0.0 0.5 0.5&
basis 0.25 0.25 0.25 basis 0.75 0.25 0.25 basis 0.25 0.75 0.25 basis 0.75 0.75 0.25 basis 0.25 0.25 0.75 basis 0.75 0.25 0.75 basis 0.25 0.75 0.75 basis 0.75 0.75 0.75
region mybox block 0 \{size\} 0 {size} 0 ${size}
create_box 2 mybox
create_atoms 2 box basis 1 2 basis 2 2 basis 3 2 basis 4 2&
basis 5 1 basis 6 1 basis 7 1 basis 8 1&
basis 9 1 basis 10 1 basis 11 1 basis 12 1
# atom 1 is Oxygen, atom 2 is Zr
group O type 1
group Zr type 2
set group O charge -2
set group Zr charge 4
mass 1 16
mass 2 91.22
neigh_modify one 20000 page 200000
############### Definition of potential #########################
pair_style reax \{hbcut\} {charge_precision}
pair_coeff * * ffield.reax 3 5
log Murnaghan_hb${hbcut}_p${charge_precision}.dat
############### Deformation of unit cell #######################
variable ax equal lx/v_size
variable ay equal ly/v_size
variable az equal lz/v_size
variable vol_per_ZrO2 equal vol/(\{size\}\*{size}*\{size\}\*4\)
\# Shift the origin of the energy scale to \(approximately\) the
\# energy of the ground\-state monoclinic structure
variable E\_per\_ZrO2 equal etotal/\({size}*\{size\}\*{size}*4)+565.76
thermo_style custom step press pxx pyy pzz etotal pe evdwl ecoul v_ax v_ay v_az v_E_per_ZrO2 v_vol_per_ZrO2
thermo_modify flush yes
thermo 1
fix DEF all deform 1 x final 0 \{L\_f\} y final 0 {L_f} z final 0 \{L\_f\} units box
run {deform_steps}
unfix DEF
ffield.reax (17.5 KB)
