+int calculate_msd(t_moldyn *moldyn,double *msd) {
+
+ int i;
+ t_atom *atom;
+ t_3dvec dist;
+ t_3dvec final_r;
+ double d2;
+ int a_cnt;
+ int b_cnt;
+
+ atom=moldyn->atom;
+ msd[0]=0;
+ msd[1]=0;
+ msd[2]=0;
+ a_cnt=0;
+ b_cnt=0;
+
+ for(i=0;i<moldyn->count;i++) {
+
+ /* care for pb crossing */
+ if(atom[i].pbc[0]|atom[i].pbc[1]|atom[i].pbc[2]) {
+ printf("[moldyn] msd pb crossings for atom %d\n",i);
+ printf(" x: %d y: %d z: %d\n",
+ atom[i].pbc[0],atom[i].pbc[1],atom[i].pbc[2]);
+ }
+ final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x;
+ final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y;
+ final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z;
+ /* calculate distance */
+ v3_sub(&dist,&final_r,&(atom[i].r_0));
+ d2=v3_absolute_square(&dist);
+
+ if(atom[i].brand) {
+ b_cnt+=1;
+ msd[1]+=d2;
+ }
+ else {
+ a_cnt+=1;
+ msd[0]+=d2;
+ }
+
+ msd[2]+=d2;
+ }
+
+ msd[0]/=a_cnt;
+ msd[1]/=b_cnt;
+ msd[2]/=moldyn->count;
+
+ return 0;
+}
+