Hello,
I am very new to LAMMPS and molecular dynamics simulations, and I am currently trying to estimate the heat capacities cv and cp of W(CO)6 at 300K.
To estimate cv, I simulated a system containing 20,000 W(CO)₆ molecules under NVT conditions for 200,000 steps. I began averaging after 40,000 steps. At each step, I divided the total energy by the number of molecules. I then divided the dataset into 10 bins, computed the average energy within each bin, and took the average of those bin averages using (cv=<E^2> - ^2 /R*T^2) where kB is 0.0019872 kcal/molK. To estimate cp, I used the formula (cp - cv = R). However, the values I obtained for both
𝐶v and cp are far from the expected values. I would appreciate any guidance or suggestions. I have also attached my LAMMPS input script in case there are errors in the setup. Thank you
units real
atom_style charge
boundary p p p
region box block 0 935.1 0 935.1 0 935.1
create_box 3 box
pair_style reax/c NULL checkqeq no
pair_coeff * * ffield.WSe2 W C O
mass 1 183.84
mass 2 12.011
mass 3 15.999
molecule wco6 WCO6.mol
create_atoms 0 random 20000 12345 box mol wco6 45678
velocity all create 300.0 12345 dist gaussian
— Energy Minimization Step —
reset_timestep 0
min_style cg
minimize 1e-4 1e-6 1000 10000
— Time Integration —
timestep 0.25
fix 2 all nvt temp 300.0 300.0 100.0
— Computes —
compute pe1 all pe
compute ke1 all ke
variable Etotal equal c_pe1+c_ke1
variable Esquared equal (c_pe1+c_ke1)^2
— Averages for Cv —
Start averaging after step 40000
fix 3 all ave/time 1000 16 16000 v_Etotal v_Esquared file etotal_avg.dat start 40000
— Thermo Output —
thermo 100
thermo_style custom step temp press vol pe ke etotal
— Dump —
dump myDump all atom 100 dump1.lammpstrj
— Run Simulation —
run 200000
write_data WCO6_20000mol_equilibrated.lmp