Cv and Cp calculation

LAMMPS LAMMPS Beginners
1

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

2

This is more a question about how to do your research than it is about the specifics of using LAMMPS. That makes it mostly off-topic for the LAMMPS forum categories. You should discuss this with your adviser or tutor or competent collaborators, i.e. people that know and care about your research.

Since this is not my area of research, I have nothing to recommend myself. Out of curiosity, I have started a search on perplexity.ai for available methodologies to compute heat capacities from MD.
Check out: https://www.perplexity.ai/search/how-do-i-compute-heat-capacity-FwkPD5GjTI6J_so4hIt90g
The response looks comprehensive and on-topic (but with AI based searches one always has to be careful), so I suggest to have a look at the references it uses (those tend to be more reliable, if from a credible source).