From 66e48ac1e6b133a607b5fbe2d4d8b7f523659021 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Mon, 6 Apr 2009 18:29:01 +0200
Subject: [PATCH] added another diff calc code

---
 diffusion_calc_ver2.c | 113 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 113 insertions(+)
 create mode 100644 diffusion_calc_ver2.c

diff --git a/diffusion_calc_ver2.c b/diffusion_calc_ver2.c
new file mode 100644
index 0000000..ef071a0
--- /dev/null
+++ b/diffusion_calc_ver2.c
@@ -0,0 +1,113 @@
+/*
+ * 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;
+}
+
-- 
2.20.1