From 30fafe0cb366808315649bd07bfaa7133ca0ecd2 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Mon, 18 Feb 2008 19:42:42 +0100
Subject: [PATCH] added diffusion code

---
 Makefile         |  6 ++--
 diff_calc_script | 74 ++++++++++++++++++++++++++++++++++++++++++++++++
 diffusion_calc.c | 58 +++++++++++++++++++++++++++++++++++++
 3 files changed, 136 insertions(+), 2 deletions(-)
 create mode 100755 diff_calc_script
 create mode 100644 diffusion_calc.c

diff --git a/Makefile b/Makefile
index b25487e..1b875fe 100644
--- a/Makefile
+++ b/Makefile
@@ -4,7 +4,7 @@ CFLAGS = -Wall
 #CFLAGS += -Wextra -pedantic
 CFLAGS += -O3
 CFLAGS += -g
-#CFLAGS += -ffloat-store
+CFLAGS += -ffloat-store
 
 CFLAGS += -DALBE
 #CFLAGS += -DTERSOFF_ORIG
@@ -22,7 +22,7 @@ DEPS = moldyn.o random/random.o list/list.o
 DEPS += potentials/lennard_jones.o potentials/harmonic_oscillator.o
 DEPS += potentials/tersoff.o potentials/albe.o
 
-all: posic sic fluctuation_calc postproc pair_correlation_calc
+all: posic sic fluctuation_calc postproc pair_correlation_calc diffusion_calc
 
 posic: $(DEPS)
 
@@ -32,6 +32,8 @@ postproc: $(DEPS)
 
 pair_correlation_calc: $(DEPS)
 
+diffusion_calc: $(DEPS)
+
 .PHONY:clean
 clean:
 	rm -vf sic postproc fluctuation_calc pair_correlation_calc *.o */*.o
diff --git a/diff_calc_script b/diff_calc_script
new file mode 100755
index 0000000..272c7bd
--- /dev/null
+++ b/diff_calc_script
@@ -0,0 +1,74 @@
+#!/bin/sh
+
+#
+# calculate diffusion coefficient within a certain time period
+# frank.zirkelbach@physik.uni-augsburg.de
+#
+
+do_it() {
+	echo "$3 `./diffusion_calc $1 | grep coefficients | \
+	      awk '{ print $3 " " $4 " " $5 }'`" >> $2
+}
+
+if [ -d $1 ]; then
+	rm -f $1/diff_coeff.txt
+	echo "processing $1 ..."
+	for file in $1/*.save; do
+		time="`basename $file | awk -F- '{ print $2 }' | \
+		       sed 's/\.save//'`"
+		if [ ! -z "$3" ]; then
+			[ "$time" -gt "$4" -o "$time" -lt "$3" ] && continue
+		fi
+		do_it $file $1/diff_coeff.txt $time
+	done
+else
+	echo "not a valid directory -> $1"
+	exit
+fi
+
+# gnuplot
+
+ddir=$1
+dfile=$ddir/diff_coeff.scr
+
+cat > $dfile <<-EOF
+set autoscale 
+unset log 
+unset label 
+set xtic auto 
+set ytic auto 
+set title 'Diffusion coefficients' 
+set xlabel 't [fs]' 
+set ylabel 'D [cm cm / s]'
+set terminal postscript eps enhanced color solid lw 1 'Helvetica' 14 
+set output '$ddir/diff_coeff.eps' 
+EOF
+
+echo -en "plot " >> $dfile
+
+komma=0
+
+if [ ! -z `echo $2 | grep a` ]; then
+	[ "$komma" = "1" ] &&
+		echo -en ", " >> $dfile
+	echo -en "\"$1/diff_coeff.txt\" u 1:2 w l t \"a\"" >> $dfile
+	komma=1
+fi
+
+if [ ! -z `echo $2 | grep b` ]; then
+	[ "$komma" = "1" ] &&
+		echo -en ", " >> $dfile
+	echo -en "\"$1/diff_coeff.txt\" u 1:3 w l t \"b\"" >> $dfile
+	komma=1
+fi
+
+if [ ! -z `echo $2 | grep t` ]; then
+	[ "$komma" = "1" ] &&
+		echo -en ", " >> $dfile
+	echo -en "\"$1/diff_coeff.txt\" u 1:4 w l t \"t\"" >> $dfile
+fi
+
+echo -en "\n" >> $dfile
+
+gnuplot $dfile
+
diff --git a/diffusion_calc.c b/diffusion_calc.c
new file mode 100644
index 0000000..576064a
--- /dev/null
+++ b/diffusion_calc.c
@@ -0,0 +1,58 @@
+/*
+ * calculation of diffusion coefficient
+ *
+ * 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"
+
+int usage(char *prog) {
+
+	printf("\nusage:\n");
+	printf("  %s <save file>\n\n",prog);
+
+	return -1;
+}
+
+int main(int argc,char **argv) {
+
+	t_moldyn moldyn;
+	int ret;
+	double dc[3];
+
+	if(argc!=2) {
+		usage(argv[0]);
+		return -1;
+	}
+
+	memset(&moldyn,0,sizeof(t_moldyn));
+
+	printf("[diffusion calc] reading save file ...\n");
+	ret=moldyn_read_save_file(&moldyn,argv[1]);
+	if(ret) {
+		printf("[diffusion calc] exit!\n");
+		return ret;
+	}
+
+	calculate_diffusion_coefficient(&moldyn,dc);
+	
+	dc[0]*=(100*100*SECOND/(METER*METER));
+	dc[1]*=(100*100*SECOND/(METER*METER));
+	dc[2]*=(100*100*SECOND/(METER*METER));
+
+	printf("diffusion coefficients: %.10f %.10f %.10f [cm^2/s]\n",
+	       dc[0],dc[1],dc[2]);
+
+	return 0;
+}
+
-- 
2.39.2