--- /dev/null
+/*
+ * calculation of diffusion coefficient (version 2)
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+
+#include "moldyn.h"
+
+#define CM 1.0e8
+
+int usage(char *prog) {
+
+ printf("\nusage:\n");
+ printf(" %s <save file 1> <save file 2>\n\n",prog);
+
+ return -1;
+}
+
+int main(int argc,char **argv) {
+
+ t_moldyn moldyn_0,moldyn_1;
+ int ret,i,a_cnt,b_cnt;
+ double dx,dy,dz,dc[3];
+
+ if(argc!=3) {
+ usage(argv[0]);
+ return -1;
+ }
+
+ memset(&moldyn_0,0,sizeof(t_moldyn));
+ memset(&moldyn_1,0,sizeof(t_moldyn));
+
+ a_cnt=0;
+ b_cnt=0;
+
+ dc[0]=0.0;
+ dc[1]=0.0;
+ dc[2]=0.0;
+
+ printf("[diffusion calc] reading save files ...\n");
+ ret=moldyn_read_save_file(&moldyn_0,argv[1]);
+ if(ret) {
+ printf("[diffusion calc] exit!\n");
+ return ret;
+ }
+ ret=moldyn_read_save_file(&moldyn_1,argv[2]);
+ if(ret) {
+ printf("[diffusion calc] exit!\n");
+ return ret;
+ }
+
+ if(moldyn_0.count!=moldyn_1.count) {
+ printf("[diffusion calc] atom count mismatch\n");
+ return -1;
+ }
+
+ for(i=0;i<moldyn_0.count;i++) {
+ dx=moldyn_0.atom[i].r.x/moldyn_0.dim.x;
+ dy=moldyn_0.atom[i].r.y/moldyn_0.dim.y;
+ dz=moldyn_0.atom[i].r.z/moldyn_0.dim.z;
+ dx-=(moldyn_1.atom[i].r.x/moldyn_1.dim.x);
+ dy-=(moldyn_1.atom[i].r.y/moldyn_1.dim.y);
+ dz-=(moldyn_1.atom[i].r.z/moldyn_1.dim.z);
+ if(dx>0.5)
+ dx-=0.5;
+ if(dx<-0.5)
+ dx+=0.5;
+ if(dy>0.5)
+ dy-=0.5;
+ if(dy<-0.5)
+ dy+=0.5;
+ if(dz>0.5)
+ dz-=0.5;
+ if(dz<-0.5)
+ dz+=0.5;
+ dx*=moldyn_1.dim.x;
+ dy*=moldyn_1.dim.y;
+ dz*=moldyn_1.dim.z;
+ if(moldyn_0.atom[i].brand) {
+ b_cnt+=1;
+ dc[1]+=dx*dx+dy*dy+dz*dz;
+ }
+ else {
+ a_cnt+=1;
+ dc[0]+=dx*dx+dy*dy+dz*dz;
+ }
+ dc[2]+=dx*dx+dy*dy+dz*dz;
+ }
+
+ dc[0]*=(1.0/(6.0*(moldyn_1.time-moldyn_0.time)*a_cnt));
+ dc[1]*=(1.0/(6.0*(moldyn_1.time-moldyn_0.time)*b_cnt));
+ dc[2]*=(1.0/(6.0*(moldyn_1.time-moldyn_0.time)*moldyn_0.count));
+
+ dc[0]*=(SECOND/(CM*CM));
+ dc[1]*=(SECOND/(CM*CM));
+ dc[2]*=(SECOND/(CM*CM));
+
+ printf("diffusion coefficients: %.10f %.10f %.10f [cm^2/s]\n",
+ dc[0],dc[1],dc[2]);
+
+ return 0;
+}
+