- for(i=0;i<moldyn->count;i++)
- virial+=v3_norm(&(atom[i].virial));
-
-printf("%f\n",virial);
- /* get pressure from virial */
- moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial;
- moldyn->p/=moldyn->volume;
-printf("%f\n",moldyn->p/(ATM));
-
- /* scale factor */
- if(moldyn->pt_scale&P_SCALE_BERENDSEN)
- scale=3*sqrt(1-(moldyn->p_ref-moldyn->p)/moldyn->p_tc);
- else
- /* should actually never be used */
- scale=pow(moldyn->p/moldyn->p_ref,1.0/3.0);
-
-printf("scale = %f\n",scale);
- /* actual scaling */
- dim->x*=scale;
- dim->y*=scale;
- dim->z*=scale;
- if(vdim->x) vdim->x=dim->x;
- if(vdim->y) vdim->y=dim->y;
- if(vdim->z) vdim->z=dim->z;
- moldyn->volume*=(scale*scale*scale);
-
- /* check whether we need a new linkcell init */
- if((dim->x/moldyn->cutoff!=lc->nx)||
- (dim->y/moldyn->cutoff!=lc->ny)||
- (dim->z/moldyn->cutoff!=lc->nx)) {
- link_cell_shutdown(moldyn);
- link_cell_init(moldyn);
+ return p;
+}
+
+double virial_sum(t_moldyn *moldyn) {
+
+ int i;
+ t_virial *virial;
+
+ /* virial (sum over atom virials) */
+ moldyn->virial=0.0;
+ moldyn->vir.xx=0.0;
+ moldyn->vir.yy=0.0;
+ moldyn->vir.zz=0.0;
+ moldyn->vir.xy=0.0;
+ moldyn->vir.xz=0.0;
+ moldyn->vir.yz=0.0;
+ for(i=0;i<moldyn->count;i++) {
+ virial=&(moldyn->atom[i].virial);
+ moldyn->virial+=(virial->xx+virial->yy+virial->zz);
+ moldyn->vir.xx+=virial->xx;
+ moldyn->vir.yy+=virial->yy;
+ moldyn->vir.zz+=virial->zz;
+ moldyn->vir.xy+=virial->xy;
+ moldyn->vir.xz+=virial->xz;
+ moldyn->vir.yz+=virial->yz;
+ }
+
+ /* global virial (absolute coordinates) */
+ //virial=&(moldyn->gvir);
+ //moldyn->gv=virial->xx+virial->yy+virial->zz;
+
+ return moldyn->virial;
+}
+
+double pressure_calc(t_moldyn *moldyn) {
+
+ /*
+ * PV = NkT + <W>
+ * with W = 1/3 sum_i f_i r_i (- skipped!)
+ * virial = sum_i f_i r_i
+ *
+ * => P = (2 Ekin + virial) / (3V)
+ */
+
+ /* assume up to date virial & up to date kinetic energy */
+
+ /* pressure (atom virials) */
+ moldyn->p=2.0*moldyn->ekin+moldyn->virial;
+ moldyn->p/=(3.0*moldyn->volume);
+
+ //moldyn->px=2.0*moldyn->ekinx+moldyn->vir.xx;
+ //moldyn->px/=moldyn->volume;
+ //moldyn->py=2.0*moldyn->ekiny+moldyn->vir.yy;
+ //moldyn->py/=moldyn->volume;
+ //moldyn->pz=2.0*moldyn->ekinz+moldyn->vir.zz;
+ //moldyn->pz/=moldyn->volume;
+
+ /* pressure (absolute coordinates) */
+ //moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
+ //moldyn->gp/=(3.0*moldyn->volume);
+
+ return moldyn->p;
+}
+
+int average_reset(t_moldyn *moldyn) {
+
+ printf("[moldyn] average reset\n");
+
+ /* update skip value */
+ moldyn->avg_skip=moldyn->total_steps;
+
+ /* kinetic energy */
+ moldyn->k_sum=0.0;
+ moldyn->k2_sum=0.0;
+
+ /* potential energy */
+ moldyn->v_sum=0.0;
+ moldyn->v2_sum=0.0;
+
+ /* temperature */
+ moldyn->t_sum=0.0;
+
+ /* virial */
+ moldyn->virial_sum=0.0;
+ //moldyn->gv_sum=0.0;
+
+ /* pressure */
+ moldyn->p_sum=0.0;
+ //moldyn->gp_sum=0.0;
+ moldyn->tp_sum=0.0;
+
+ return 0;
+}
+
+int average_and_fluctuation_calc(t_moldyn *moldyn) {
+
+ int denom;
+
+ if(moldyn->total_steps<moldyn->avg_skip)
+ return 0;
+
+ denom=moldyn->total_steps+1-moldyn->avg_skip;
+
+ /* assume up to date energies, temperature, pressure etc */
+
+ /* kinetic energy */
+ moldyn->k_sum+=moldyn->ekin;
+ moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
+ moldyn->k_avg=moldyn->k_sum/denom;
+ moldyn->k2_avg=moldyn->k2_sum/denom;
+ moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg);
+
+ /* potential energy */
+ moldyn->v_sum+=moldyn->energy;
+ moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
+ moldyn->v_avg=moldyn->v_sum/denom;
+ moldyn->v2_avg=moldyn->v2_sum/denom;
+ moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg);
+
+ /* temperature */
+ moldyn->t_sum+=moldyn->t;
+ moldyn->t_avg=moldyn->t_sum/denom;
+
+ /* virial */
+ moldyn->virial_sum+=moldyn->virial;
+ moldyn->virial_avg=moldyn->virial_sum/denom;
+ //moldyn->gv_sum+=moldyn->gv;
+ //moldyn->gv_avg=moldyn->gv_sum/denom;
+
+ /* pressure */
+ moldyn->p_sum+=moldyn->p;
+ moldyn->p_avg=moldyn->p_sum/denom;
+ //moldyn->gp_sum+=moldyn->gp;
+ //moldyn->gp_avg=moldyn->gp_sum/denom;
+ moldyn->tp_sum+=moldyn->tp;
+ moldyn->tp_avg=moldyn->tp_sum/denom;
+
+ return 0;
+}
+
+int get_heat_capacity(t_moldyn *moldyn) {
+
+ double temp2,ighc;
+
+ /* averages needed for heat capacity calc */
+ if(moldyn->total_steps<moldyn->avg_skip)
+ return 0;
+
+ /* (temperature average)^2 */
+ temp2=moldyn->t_avg*moldyn->t_avg;
+ printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
+ moldyn->t_avg);
+
+ /* ideal gas contribution */
+ ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
+ printf(" ideal gas contribution: %f\n",
+ ighc/moldyn->mass*KILOGRAM/JOULE);
+
+ /* specific heat for nvt ensemble */
+ moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc;
+ moldyn->c_v_nvt/=moldyn->mass;
+
+ /* specific heat for nve ensemble */
+ moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2)));
+ moldyn->c_v_nve/=moldyn->mass;
+
+ printf(" NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
+ printf(" NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
+printf(" --> <dV2> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->count*K_B2*moldyn->t_avg*moldyn->t_avg*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM)));
+
+ return 0;
+}
+
+double thermodynamic_pressure_calc(t_moldyn *moldyn) {
+
+ t_3dvec dim;
+ //t_3dvec *tp;
+ double h,dv;
+ double y0,y1;
+ double su,sd;
+ t_atom *store;
+
+ /*
+ * dU = - p dV
+ *
+ * => p = - dU/dV
+ *
+ */
+
+ /* store atomic configuration + dimension */
+ store=malloc(moldyn->count*sizeof(t_atom));
+ if(store==NULL) {
+ printf("[moldyn] allocating store mem failed\n");
+ return -1;