added diffusion code

diff --git a/Makefile b/Makefile
index b25487e..1b875fe 100644 (file)
--- 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 (executable)
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 (file)
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;
+}
+