Dear everyone,
Recently I am working on how to calculate Cp and Cv from MD simulation from energy fluctuation. I found in the book ‘Computer Simulation of Liquids’.
Based on the above equation, I wrote my own LAMMPS commands to calculate Cp of water at 300K&1atm. However, I just cannot get the right results. I guess the main reason is the unit conversion process during calculation. Here is my related LAMMPS commands, (my system has 900 water molecules)
variable h equal abs(enthalpy)*4184/900 # J/mol
variable h2 equal v_h^2
fix 2 all ave/time 10 1000 10000 v_h
fix 3 all ave/time 10 1000 10000 v_h2
variable cp equal (f_3-f_2^2)/NA/kB/T/T/18 # kJ/kgK (I removed dollar sign for simplicity)
And the output is like,
Step Temp Press KinEng PotEng TotEng Density v_cp
210000 304.97165 720.11485 2453.5599 -9122.7096 -6669.1497 1.0418024 0.16568487
220000 298.43008 -1222.8428 2400.9316 -9077.2395 -6676.3079 1.0347138 0.1921633
230000 297.43112 -558.0904 2392.8947 -9137.544 -6744.6493 1.0431676 0.21482302
240000 296.86561 1863.72 2388.3451 -9093.3951 -6705.05 1.0395024 0.15437286
250000 293.22998 -97.355582 2359.0957 -9170.5971 -6811.5014 1.0494803 0.25051753
260000 293.19078 1523.8496 2358.7803 -9200.3628 -6841.5826 1.0608247 0.23747013
270000 303.63927 -101.09751 2442.8406 -9129.2868 -6686.4461 1.0344554 0.13832415
280000 294.13219 889.73517 2366.3541 -9251.0854 -6884.7312 1.0536651 0.23688666
290000 303.75677 -438.42439 2443.7859 -9037.7087 -6593.9228 1.0266251 0.27586393
300000 307.04574 -994.36334 2470.2463 -9120.3572 -6650.1109 1.0352707 0.18916265
The value of Cp should be 4.184kJ/kgK for water at 300K&1atm. But what I got is about 0.2 or so.
Could anyone help me fix it right? I think the problem should be in the formula of variable cp.
Thanks a lot and may you have a nice day!