I want to calculate enthalpy based on the paper referenced here.My potential is the same as the one used by the author, and I am using the LAMMPS code provided by the author. however, when I calculate the enthalpy difference for T=1000 and T=300 (i.e., h(1000)−h(300)), I do not obtain the same result as reported in the paper. (i subtract the average enthalpy by the number of atoms and then multiply the result by 96.49 to convert from eV to kJ/mol.)
here is my code:
# set periodic boundaries, units and thero-baro damping
boundary p p p
atom_style full
units metal
variable THERMO_DAMP equal 0.1
variable BARO_DAMP equal 0.5
# import UO2 single cell fluorite structure and assign atom labels
read_data UO2-single_cell.lmpstruct
variable O equal 1
variable U equal 2
# set charges
set type $O charge -1.1104
set type $U charge 2.2208
# define interatomic potential via coulombic and embed_UO2.fs tabulation
kspace_style pppm 1.0e-5
variable SR_CUTOFF equal 11.0
pair_style hybrid/overlay coul/long ${SR_CUTOFF} eam/alloy
pair_coeff * * coul/long
pair_coeff * * eam/alloy CeThUNpPuAmCmO.eam.alloy O U
#Perform minimization
fix boxrelax_fix all box/relax aniso 0.0 vmax 0.001
thermo 1
thermo_style custom step cpu press etotal lx cella cellb cellc cellalpha cellbeta cellgamma
min_style sd
min_modify line quadratic dmax 0.05
minimize 1.0e-25 1.0e-25 1000 10000
unfix boxrelax_fix
# Create 4x4x4 supercell containing 768 atoms
replicate 4 4 4
# create initial velocities and set timestep
velocity all create 300 160278
timestep 0.002
# Set-up thermostat
fix thermostat_fix all npt temp 300 300 ${THERMO_DAMP} aniso 0.0 0.0 ${BARO_DAMP}
thermo_style one
thermo 10
neigh_modify delay 0 every 1 check yes
# equilibrate for 50ps
run 25000
Here is the enthalpy per step from my simulation.