+
+int get_heat_capacity(t_moldyn *moldyn) {
+
+ double temp2,mass,ighc;
+ int i;
+
+ /* (temperature average)^2 */
+ temp2=2.0*moldyn->k_m*EV/(3.0*K_BOLTZMANN);
+ printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",temp2);
+ temp2*=temp2;
+
+ /* total mass */
+ mass=0.0;
+ for(i=0;i<moldyn->count;i++)
+ mass+=moldyn->atom[i].mass;
+
+ /* ideal gas contribution */
+ ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
+ printf(" ideal gas contribution: %f\n",ighc/mass*KILOGRAM/JOULE);
+
+ moldyn->c_v_nvt=moldyn->dp2_m*moldyn->count*moldyn->count*EV/(K_BOLTZMANN*temp2)+ighc;
+ moldyn->c_v_nvt/=mass;
+ moldyn->c_v_nve=ighc/(1.0-(moldyn->dp2_m*moldyn->count*moldyn->count*EV/(ighc*K_BOLTZMANN*temp2)));
+ moldyn->c_v_nve/=mass;
+
+ printf(" NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
+ printf(" NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
+
+ return 0;
+}