From 7acc7a495744e10db14c3d85c9829ceeddfea61c Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Mon, 17 Aug 2009 15:58:39 +0200
Subject: [PATCH 01/16] new defaults (+ alternatives) for the Makefile

---
 Makefile | 15 ++++++++++++---
 1 file changed, 12 insertions(+), 3 deletions(-)

diff --git a/Makefile b/Makefile
index 224fe72..fc5705a 100644
--- a/Makefile
+++ b/Makefile
@@ -1,29 +1,38 @@
 CC = gcc-4.3
+#CC = gcc-3.4
 
 CFLAGS = -Wall -Winline
 #CFLAGS += -Wextra -pedantic
+
 CFLAGS += -O3 -march=native -msse2 -mfpmath=sse
+#CFLAGS += -O3 -march=athlon64
+
 CFLAGS += -g
 #CFLAGS += -pg
 #CFLAGS += -ffloat-store
 
 #CFLAGS += -DPARALLEL -fopenmp
 #CFLAGS += -DPTHREADS -lpthread
-CFLAGS += -DVISUAL_THREAD -lpthread
+#CFLAGS += -DVISUAL_THREAD -lpthread
 #CFLAGS += -DPTHREADS -DVISUAL_THREAD -lpthread
 
-CFLAGS += -DALBE
+#CFLAGS += -DALBE
 CFLAGS += -DALBE_FAST
 #CFLAGS += -DTERSOFF_ORIG
 
 #CFLAGS += -DSTATIC_LISTS
 CFLAGS += -DLOWMEM_LISTS
 
+#CFLAGS += -DPDEBUG
 #CFLAGS += -DDEBUG
-#CFLAGS += -DDSTART=-1 -DDEND=40 -DDATOM=0
+#CFLAGS += -DDSTART=50 -DDEND=60 -DDATOM=0
 #CFLAGS += -DVDEBUG
 
+#CFLAGS += -DCONSTRAINT_110_5832
+#CFLAGS += -DQUENCH
+
 LDFLAGS = -lm
+
 #LDFLAGS += -lc_p
 #LDFLAGS += -lefence
 
-- 
2.20.1


From dc980b4d12e3ad4f2421e09db6acd4ce83132dec Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Mon, 17 Aug 2009 15:59:13 +0200
Subject: [PATCH 02/16] a modified posselt variation + the right one commented

---
 calc_delta_e | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/calc_delta_e b/calc_delta_e
index 6a97cc5..2b4130e 100755
--- a/calc_delta_e
+++ b/calc_delta_e
@@ -20,5 +20,6 @@ echo "-------------------------------------------------------------------------"
 echo "$si_cnt $c_cnt $e1 $e0" | \
 	awk '{ print "   Gao: "($3-$4)*($1+$2) }'
 echo "$si_cnt $c_cnt $musi $muc $e1 $music" | \
-	awk '{ print "   Posselt: "$5*($1+$2)-0.5*($1+$2)*$6-0.5*($1-$2)*($3-$4) }'
+	#awk '{ print "   Posselt: "$5*($1+$2)-0.5*($1+$2)*$6-0.5*($1-$2)*($3-$4) }'
+	awk '{ print "   Posselt: "$5*($1+$2)-$1*$3-$2*$4) }'
 
-- 
2.20.1


From f6f9984d115aeabdbbcbfec80a2aba2a35e601fb Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Mon, 17 Aug 2009 16:00:00 +0200
Subject: [PATCH 03/16] init of a and b values

---
 random/random.c | 1 +
 1 file changed, 1 insertion(+)

diff --git a/random/random.c b/random/random.c
index 902cdef..292bc69 100644
--- a/random/random.c
+++ b/random/random.c
@@ -109,6 +109,7 @@ double rand_get_gauss(t_random *random) {
 		return random->gauss;
 	}
 
+	a=0; b=0;
 	w=0;
 	while((w>=1.0)||(w==0.0)) {
 		a=-2.0*rand_get_double(random)+1.0;
-- 
2.20.1


From 5d016e95f7763347385ea545fdc5faadc80af25e Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Mon, 17 Aug 2009 16:02:07 +0200
Subject: [PATCH 04/16] put logs to logs directory + -p for the mkdir commands

---
 runmd | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/runmd b/runmd
index d50d9ae..541734b 100755
--- a/runmd
+++ b/runmd
@@ -8,13 +8,14 @@ if [ ! -f ./config ]; then
 	exit
 fi
 
-[ ! -d $1 ] && mkdir $1
+[ ! -d $1 ] && mkdir -p $1
 
 ./clean $1
 
 cp -v config $1/config
 
-time ./mdrun -c ./config -s $1
+mkdir -p logs
+time ./mdrun -c ./config -s $1 | tee logs/run_`basename $1`.log
 
 if [ "$?" == "0" ]; then
 	#./perms
@@ -24,7 +25,7 @@ if [ "$?" == "0" ]; then
 
 		# center unit cell
 		./visualize -w 640 -h 480 -d $1 \
-		            -nll -0.56 -0.56 -0.56 -fur 0.56 0.56 0.56 \
+		            -nll -0.56 -0.56 -0.76 -fur 0.56 0.56 0.56 \
 			    -b -0.5 -0.5 -0.5 0.5 0.5 0.5 \
 			    -c -0.2 -2.0 0.6 -L 0 0 -0.1 \
 			    -r 0.6 -B 0.1
-- 
2.20.1


From 2998bed05e84541809a34de696eeb1a9267ffe4d Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Mon, 17 Aug 2009 16:02:57 +0200
Subject: [PATCH 05/16] no verbose printfs for total v attrr changes

---
 mdrun.c | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/mdrun.c b/mdrun.c
index 31e5f0a..f5acaa1 100644
--- a/mdrun.c
+++ b/mdrun.c
@@ -1108,7 +1108,9 @@ int chaatr(t_moldyn *moldyn,t_mdrun *mdrun) {
 			if(cap->z1<atom->r.z)
 				continue;
 		}
-		printf("  changing attributes of atom %d (0x%x)\n",i,cap->attr);
+		if(!(cap->type&CHAATTR_TOTALV))
+			printf("  changing attributes of atom %d (0x%x)\n",
+			       i,cap->attr);
 		atom->attr=cap->attr;
 	}
 
-- 
2.20.1


From fbf49170d213d948e9523d727a76b5337de90a05 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Mon, 17 Aug 2009 16:03:35 +0200
Subject: [PATCH 06/16] QUENCH feature (well, or sth like that ...)

---
 moldyn.c | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/moldyn.c b/moldyn.c
index d3d63eb..056047c 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -2283,7 +2283,9 @@ int velocity_verlet(t_moldyn *moldyn) {
 			delta.y=-delta.x;
 		}
 #endif
+#ifndef QUENCH
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
+#endif
 		v3_scale(&delta,&(atom[i].f),h*tau_square);
 #ifdef CONSTRAINT_110_5832
 		if(i==5832) {
-- 
2.20.1


From f0e4c6422ec1aff0cf86597fef919335bba75c1b Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Mon, 17 Aug 2009 17:16:20 +0200
Subject: [PATCH 07/16] added offset feature

---
 mdrun.c  | 17 ++++++++++++++---
 moldyn.c | 11 ++++++++---
 moldyn.h |  3 ++-
 3 files changed, 24 insertions(+), 7 deletions(-)

diff --git a/mdrun.c b/mdrun.c
index f5acaa1..062224b 100644
--- a/mdrun.c
+++ b/mdrun.c
@@ -393,6 +393,14 @@ int mdrun_parse_config(t_mdrun *mdrun) {
 					continue;
 
 				}
+				// offset
+				if(word[i][0]=='o') {
+					fp.o_params.o.x=atof(word[++i])*fp.lc;
+					fp.o_params.o.y=atof(word[++i])*fp.lc;
+					fp.o_params.o.z=atof(word[++i])*fp.lc;
+					fp.o_params.use=1;
+					continue;
+				}
 				i+=1;
 			}
 			add_stage(mdrun,STAGE_FILL,&fp);
@@ -1317,7 +1325,8 @@ int mdrun_hook(void *ptr1,void *ptr2) {
 					               fp->lx,fp->ly,fp->lz,
 					               &o,
 					               &(fp->p_params),
-					               &(fp->d_params));
+					               &(fp->d_params),
+					               &(fp->o_params));
 					o.x+=0.25*fp->lc;
 					o.y=o.x;
 					o.z=o.x;
@@ -1328,7 +1337,8 @@ int mdrun_hook(void *ptr1,void *ptr2) {
 					               fp->lx,fp->ly,fp->lz,
 					               &o,
 					               &(fp->p_params),
-					               &(fp->d_params));
+					               &(fp->d_params),
+					               &(fp->o_params));
 					break;
 
 					default:
@@ -1341,7 +1351,8 @@ int mdrun_hook(void *ptr1,void *ptr2) {
 					               fp->lx,fp->ly,fp->lz,
 					               NULL,
 					               &(fp->p_params),
-					               &(fp->d_params));
+					               &(fp->d_params),
+					               &(fp->o_params));
 					break;
 				}
 				moldyn_bc_check(moldyn);
diff --git a/moldyn.c b/moldyn.c
index 056047c..15dcede 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -520,7 +520,8 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin,
-                   t_part_params *p_params,t_defect_params *d_params) {
+                   t_part_params *p_params,t_defect_params *d_params,
+                   t_offset_params *o_params) {
 
 	int new,count;
 	int ret;
@@ -592,6 +593,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
 
 	switch(type) {
 		case CUBIC:
+			if(o_params->use)
+				v3_add(&orig,&orig,&(o_params->o));
 			set_nn_dist(moldyn,lc);
 			ret=cubic_init(a,b,c,lc,atom,&orig,p_params,d_params);
 			strcpy(name,"cubic");
@@ -599,6 +602,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
 		case FCC:
 			if(!origin)
 				v3_scale(&orig,&orig,0.5);
+			if(o_params->use)
+				v3_add(&orig,&orig,&(o_params->o));
 			set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
 			ret=fcc_init(a,b,c,lc,atom,&orig,p_params,d_params);
 			strcpy(name,"fcc");
@@ -606,6 +611,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
 		case DIAMOND:
 			if(!origin)
 				v3_scale(&orig,&orig,0.25);
+			if(o_params->use)
+				v3_add(&orig,&orig,&(o_params->o));
 			set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
 			ret=diamond_init(a,b,c,lc,atom,&orig,p_params,d_params);
 			strcpy(name,"diamond");
@@ -806,8 +813,6 @@ int del_atom(t_moldyn *moldyn,int tag) {
 			case DEFECT_STYPE_DB_Z:\
 				d_o.z=d_params->od;\
 				d_d.z=d_params->dd;\
-d_d.x=0.9;\
-d_d.y=0.9;\
 				break;\
 			case DEFECT_STYPE_DB_R:\
 				break;\
diff --git a/moldyn.h b/moldyn.h
index 6310416..bdb3758 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -444,7 +444,8 @@ int moldyn_log_shutdown(t_moldyn *moldyn);
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin,
-                   t_part_params *p_params,t_defect_params *d_params);
+                   t_part_params *p_params,t_defect_params *d_params,
+                   t_offset_params *o_params);
 int add_atom(t_moldyn *moldyn,int element,u8 brand,u8 attr,
              t_3dvec *r,t_3dvec *v);
 int del_atom(t_moldyn *moldyn,int tag);
-- 
2.20.1


From be1d19844c16b2d49a376c8afc574fe468cbb2fb Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Thu, 20 Aug 2009 11:14:32 +0200
Subject: [PATCH 08/16] added trafo c code

---
 vasp_tools/tXp.c | 214 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 214 insertions(+)
 create mode 100644 vasp_tools/tXp.c

diff --git a/vasp_tools/tXp.c b/vasp_tools/tXp.c
new file mode 100644
index 0000000..54fe5e2
--- /dev/null
+++ b/vasp_tools/tXp.c
@@ -0,0 +1,214 @@
+/*
+ * tXp.c
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <string.h>
+#include <unistd.h>
+#include <math.h>
+
+int get_line(int fd,char *line,int max) {
+
+        int count,ret;
+
+        count=0;
+
+        while(1) {
+                if(count==max) return count;
+                ret=read(fd,line+count,1);
+                if(ret<=0) return ret;
+                if(line[count]=='\n') {
+                        memset(line+count,0,max-count-1);
+                        //line[count]='\0';
+                        return count+1;
+                }
+                count+=1;
+        }
+}
+
+
+int main(int argc,char **argv) {
+
+	double x1,x2,x3,y1,y2,y3,z1,z2,z3;
+	double x_1,x_2,x_3,y_1,y_2,y_3,z_1,z_2,z_3;
+	double X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3;
+	double theta;
+	double costheta,sintheta;
+	double normx,normy,normz;
+	int posr,poso;
+	char file[128],buf[256];
+	char *wptr;
+	char t1,t2,t3;
+	int cnt;
+
+	posr=open("POSCAR",O_RDONLY);
+	if(posr<0) {
+		perror("open POSCAR (read) file\n");
+		return posr;
+	}
+
+	theta=(atof(argv[1])/180.0*M_PI);
+	costheta=cos(theta);
+	sintheta=sin(theta);
+
+	snprintf(file,127,"POSCAR.X%f",theta);
+	poso=open(file,O_WRONLY|O_CREAT);
+	if(poso<0) {
+		perror("open POSCAR (write) file\n");
+		return poso;
+	}
+	
+	// first line
+	cnt=get_line(posr,buf,256);
+	buf[cnt-1]='\n';
+	write(poso,buf,cnt);
+
+	// second line
+	cnt=get_line(posr,buf,256);
+	buf[cnt-1]='\n';
+	write(poso,buf,cnt);
+
+	// basis
+	cnt=get_line(posr,buf,256);
+	wptr=strtok(buf," ");
+	x1=atof(wptr);
+	wptr=strtok(NULL," ");
+	x2=atof(wptr);
+	wptr=strtok(NULL," ");
+	x3=atof(wptr);
+
+	cnt=get_line(posr,buf,256);
+	wptr=strtok(buf," ");
+	y1=atof(wptr);
+	wptr=strtok(NULL," ");
+	y2=atof(wptr);
+	wptr=strtok(NULL," ");
+	y3=atof(wptr);
+
+	cnt=get_line(posr,buf,256);
+	wptr=strtok(buf," ");
+	z1=atof(wptr);
+	wptr=strtok(NULL," ");
+	z2=atof(wptr);
+	wptr=strtok(NULL," ");
+	z3=atof(wptr);
+
+	/* norm */
+	normx=sqrt(x1*x1+x2*x2+x3*x3);
+	normy=sqrt(y1*y1+y2*y2+y3*y3);
+	normz=sqrt(z1*z1+z2*z2+z3*z3);
+
+	/* basis in given basis */
+	X1=(x1*x1+x2*x2+x3*x3)/normx;
+	X2=(y1*x1+y2*x2+y3*x3)/normy;
+	X3=(z1*x1+z2*x2+z3*x3)/normz;
+
+	Y1=(x1*y1+x2*y2+x3*y3)/normx;
+	Y2=(y1*y1+y2*y2+y3*y3)/normy;
+	Y3=(z1*y1+z2*y2+z3*y3)/normz;
+
+	Z1=(x1*z1+x2*z2+x3*z3)/normx;
+	Z2=(y1*z1+y2*z2+y3*z3)/normy;
+	Z3=(z1*z1+z2*z2+z3*z3)/normz;
+
+	printf("Basis expressed by itself:\n");
+	printf("     %f     %f     %f\n",X1,Y1,Z1);
+	printf(" x = %f y = %f z = %f\n",X2,Y2,Z2);
+	printf("     %f     %f     %f\n",X3,Y3,Z3);
+
+	/* transformed basis */
+	x_1=X1;
+	x_2=costheta*X2-sintheta*X3;
+	x_3=sintheta*X2+costheta*X3;
+
+	y_1=Y1;
+	y_2=costheta*Y2-sintheta*Y3;
+	y_3=sintheta*Y2+costheta*Y3;
+
+	z_1=Z1;
+	z_2=costheta*Z2-sintheta*Z3;
+	z_3=sintheta*Z2+costheta*Z3;
+
+	printf("Transformed basis in the given basis:\n");
+	printf("     %f     %f     %f\n",x_1,y_1,z_1);
+	printf(" x = %f y = %f z = %f\n",x_2,y_2,z_2);
+	printf("     %f     %f     %f\n",x_3,y_3,z_3);
+
+	/* transformed basis in cartesian coordinates */
+	X1=(x1/normx*x_1+y1/normy*x_2+z1/normz*x_3);
+	X2=(x2/normx*x_1+y2/normy*x_2+z2/normz*x_3);
+	X3=(x3/normx*x_1+y3/normy*x_2+z3/normz*x_3);
+
+	Y1=(x1/normx*y_1+y1/normy*y_2+z1/normz*y_3);
+	Y2=(x2/normx*y_1+y2/normy*y_2+z2/normz*y_3);
+	Y3=(x3/normx*y_1+y3/normy*y_2+z3/normz*y_3);
+
+	Z1=(x1/normx*z_1+y1/normy*z_2+z1/normz*z_3);
+	Z2=(x2/normx*z_1+y2/normy*z_2+z2/normz*z_3);
+	Z3=(x3/normx*z_1+y3/normy*z_2+z3/normz*z_3);
+
+	printf("Transformed basis in cartesian coordinates:\n");
+	printf("     %f     %f     %f\n",X1,Y1,Z1);
+	printf(" x = %f y = %f z = %f\n",X2,Y2,Z2);
+	printf("     %f     %f     %f\n",X3,Y3,Z3);
+
+	dprintf(poso," %f %f %f\n",X1,X2,X3);
+	dprintf(poso," %f %f %f\n",Y1,Y2,Y3);
+	dprintf(poso," %f %f %f\n",Z1,Z2,Z3);
+
+	// 6th line
+	cnt=get_line(posr,buf,256);
+	buf[cnt-1]='\n';
+	write(poso,buf,cnt);
+
+	// 7th line
+	cnt=get_line(posr,buf,256);
+	buf[cnt-1]='\n';
+	write(poso,buf,cnt);
+
+	// 8th line
+	cnt=get_line(posr,buf,256);
+	buf[cnt-1]='\n';
+	write(poso,buf,cnt);
+
+	while(1) {
+		cnt=get_line(posr,buf,256);
+		if(cnt<=0)
+			break;
+		wptr=strtok(buf," ");
+		x1=atof(wptr);
+		wptr=strtok(NULL," ");
+		x2=atof(wptr);
+		wptr=strtok(NULL," ");
+		x3=atof(wptr);
+
+		wptr=strtok(NULL," ");
+		t1=wptr[0];
+		wptr=strtok(NULL," ");
+		t2=wptr[0];
+		wptr=strtok(NULL," ");
+		t3=wptr[0];
+
+		X1=x1;
+		X2=costheta*x2+sintheta*x3;
+		X3=costheta*x3-sintheta*x2;
+
+		dprintf(poso," %f %f %f %c %c %c\n",X1,X2,X3,t1,t2,t3);
+	}
+
+	close(posr);
+	close(poso);
+
+	printf("done!\n");
+
+	return 0;
+}
+
-- 
2.20.1


From 764be9105e36d59a0b9c910503b68818e849a982 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Thu, 20 Aug 2009 11:15:43 +0200
Subject: [PATCH 09/16] more trafo (still sth is wrong!)

---
 vasp_tools/trafoXposcar | 101 +++++++++++++++++++++++++++++++++++++---
 1 file changed, 95 insertions(+), 6 deletions(-)

diff --git a/vasp_tools/trafoXposcar b/vasp_tools/trafoXposcar
index 1aebc4f..1cd5286 100755
--- a/vasp_tools/trafoXposcar
+++ b/vasp_tools/trafoXposcar
@@ -2,10 +2,10 @@
 
 echo "parsing POSCAR file ..."
 
-trg=POSCAR.transformed
-
 theta=$1
 
+trg=POSCAR.Xtrafo$theta
+
 sicnt=`sed -n 6p POSCAR | awk '{ print $1 }'`
 ccnt=`sed -n 6p POSCAR | awk '{ print $2 }'`
 
@@ -35,27 +35,116 @@ echo
 echo "  -----"
 echo "  Trafo: x axis, $theta degree"
 echo "  -----"
+echo
+
+# determine trafo of cartesian coordinates to this one
+normx=`echo $x1 $x2 $x3 | awk '{ print sqrt($1*$1+$2*$2+$3*$3) }'`
+normy=`echo $y1 $y2 $y3 | awk '{ print sqrt($1*$1+$2*$2+$3*$3) }'`
+normz=`echo $z1 $z2 $z3 | awk '{ print sqrt($1*$1+$2*$2+$3*$3) }'`
+echo $normx $normz $normy
+t11=`echo $x1 $normx | awk '{ print $1/$2 }'`
+t12=`echo $x2 $normx | awk '{ print $1/$2 }'`
+t13=`echo $x3 $normx | awk '{ print $1/$2 }'`
+t21=`echo $y1 $normy | awk '{ print $1/$2 }'`
+t22=`echo $y2 $normy | awk '{ print $1/$2 }'`
+t23=`echo $y3 $normy | awk '{ print $1/$2 }'`
+t31=`echo $z1 $normz | awk '{ print $1/$2 }'`
+t32=`echo $z2 $normz | awk '{ print $1/$2 }'`
+t33=`echo $z3 $normz | awk '{ print $1/$2 }'`
+echo "  Matrix from cartesian to used coordinates:"
+echo "  | $t11 $t21 $t31 |"
+echo "  | $t12 $t22 $t32 |"
+echo "  | $t13 $t23 $t33 |"
+echo
+i11=$t11
+i12=$t21
+i13=$t31
+i21=$t12
+i22=$t22
+i23=$t32
+i31=$t13
+i32=$t23
+i33=$t33
+echo "  Inverse matrix (i=t^T):"
+echo "  | $i11 $i21 $i31 |"
+echo "  | $i12 $i22 $i32 |"
+echo "  | $i13 $i23 $i33 |"
+echo
+# i * basis
+X1=`echo $i11 $i21 $i31 $x1 $x2 $x3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+X2=`echo $i12 $i22 $i32 $x1 $x2 $x3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+X3=`echo $i13 $i23 $i33 $x1 $x2 $x3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+
+Y1=`echo $i11 $i21 $i31 $y1 $y2 $y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+Y2=`echo $i12 $i22 $i32 $y1 $y2 $y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+Y3=`echo $i13 $i23 $i33 $y1 $y2 $y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+
+Z1=`echo $i11 $i21 $i31 $z1 $z2 $z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+Z2=`echo $i12 $i22 $i32 $z1 $z2 $z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+Z3=`echo $i13 $i23 $i33 $z1 $z2 $z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+echo "  Basis (in rotation system):"
+echo "      $X1      $Y1      $Z1"
+echo "  x = $X2  y = $Y2  z = $Z2"
+echo "      $X3      $Y3      $Z3"
+echo
 
 costheta=`echo $theta | awk '{ print cos($1*3.1415927/180.0) }'`
 sintheta=`echo $theta | awk '{ print sin($1*3.1415927/180.0) }'`
 
+x1=$X1
+x2=$X2
+x3=$X3
 X1=$x1
 X2=`echo $x2 $x3 $costheta $sintheta | awk '{ print $3*$1-$4*$2 }'`
 X3=`echo $x2 $x3 $costheta $sintheta | awk '{ print $4*$1+$3*$2 }'`
 
+y1=$Y1
+y2=$Y2
+y3=$Y3
 Y1=$y1
 Y2=`echo $y2 $y3 $costheta $sintheta | awk '{ print $3*$1-$4*$2 }'`
 Y3=`echo $y2 $y3 $costheta $sintheta | awk '{ print $4*$1+$3*$2 }'`
 
+z1=$Z1
+z2=$Z2
+z3=$Z3
 Z1=$z1
 Z2=`echo $z2 $z3 $costheta $sintheta | awk '{ print $3*$1-$4*$2 }'`
 Z3=`echo $z2 $z3 $costheta $sintheta | awk '{ print $4*$1+$3*$2 }'`
 
+echo "  Transformed basis (in rotation system):"
+echo "       $X1       $Y1       $Z1"
+echo "  x' = $X2  y' = $Y2  z' = $Z2"
+echo "       $X3       $Y3       $Z3"
 echo
-echo "  Transformed basis:"
-echo "      $X1      $Y1      $Z1"
-echo "  x = $X2  y = $Y2  z = $Z2"
-echo "      $X3      $Y3      $Z3"
+
+# t * basis
+x1=`echo $t11 $t21 $t31 $X1 $X2 $X3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+x2=`echo $t12 $t22 $t32 $X1 $X2 $X3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+x3=`echo $t13 $t23 $t33 $X1 $X2 $X3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+X1=$x1
+X2=$x2
+X3=$x3
+
+y1=`echo $t11 $t21 $t31 $Y1 $Y2 $Y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+y2=`echo $t12 $t22 $t32 $Y1 $Y2 $Y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+y3=`echo $t13 $t23 $t33 $Y1 $Y2 $Y3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+Y1=$y1
+Y2=$y2
+Y3=$y3
+
+z1=`echo $t11 $t21 $t31 $Z1 $Z2 $Z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+z2=`echo $t12 $t22 $t32 $Z1 $Z2 $Z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+z3=`echo $t13 $t23 $t33 $Z1 $Z2 $Z3 | awk '{ print $1*$4+$2*$5+$3*$6 }'`
+Z1=$z1
+Z2=$z2
+Z3=$z3
+
+
+echo "  Transformed basis (cartesian coordinates):"
+echo "       $X1       $Y1       $Z1"
+echo "  x' = $X2  y' = $Y2  z' = $Z2"
+echo "       $X3       $Y3       $Z3"
 echo
 
 cnt=0
-- 
2.20.1


From 7b97b630ab3bee5cc967142d1a4bfedd29005179 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Thu, 20 Aug 2009 13:02:22 +0200
Subject: [PATCH 10/16] more easy ...

---
 vasp_tools/poscar2moldyn | 15 ++++++---------
 1 file changed, 6 insertions(+), 9 deletions(-)

diff --git a/vasp_tools/poscar2moldyn b/vasp_tools/poscar2moldyn
index 29a527e..b0ad584 100755
--- a/vasp_tools/poscar2moldyn
+++ b/vasp_tools/poscar2moldyn
@@ -50,15 +50,12 @@ tail -${total} $file | \
 			type="C"
 			color="Gray"
 		fi
-		((X=0))
-		((Y=0))
-		((Z=0))
-		X=`echo $X $x $y $z $x1 $y1 $z1 | \
-			awk '{ print $1+($2*$5)+($3*$6)+($4*$7) }'`
-		Y=`echo $Y $x $y $z $x2 $y2 $z2 | \
-			awk '{ print $1+($2*$5)+($3*$6)+($4*$7) }'`
-		Z=`echo $Z $x $y $z $x3 $y3 $z3 | \
-			awk '{ print $1+($2*$5)+($3*$6)+($4*$7) }'`
+		X=`echo $x $y $z $x1 $y1 $z1 | \
+			awk '{ print $1*$4+$2*$5+$3*$6 }'`
+		Y=`echo $x $y $z $x2 $y2 $z2 | \
+			awk '{ print $1*$4+$2*$5+$3*$6 }'`
+		Z=`echo $x $y $z $x3 $y3 $z3 | \
+			awk '{ print $1*$4+$2*$5+$3*$6 }'`
 		X=`echo $lc $X | awk '{ print $1*$2 }'`
 		Y=`echo $lc $Y | awk '{ print $1*$2 }'`
 		Z=`echo $lc $Z | awk '{ print $1*$2 }'`
-- 
2.20.1


From 0b7c1bce98bbad88ec20df66abb513348e37adbe Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Thu, 20 Aug 2009 13:02:54 +0200
Subject: [PATCH 11/16] added a temp check (still sth is wrong!)

---
 vasp_tools/tXp.c | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/vasp_tools/tXp.c b/vasp_tools/tXp.c
index 54fe5e2..6b6e6a1 100644
--- a/vasp_tools/tXp.c
+++ b/vasp_tools/tXp.c
@@ -201,6 +201,14 @@ int main(int argc,char **argv) {
 		X2=costheta*x2+sintheta*x3;
 		X3=costheta*x3-sintheta*x2;
 
+		/* check! */
+		double tmp1,tmp2,tmp3;
+		tmp1=(X1*x_1+X2*y_1+X3*z_1)/normx;
+		tmp2=(X1*x_2+X2*y_2+X3*z_2)/normy;
+		tmp3=(X1*x_3+X2*y_3+X3*z_3)/normz;
+		printf("%f %f %f - %f %f %f | %f %f %f\n",
+		       x1,x2,x3,tmp1,tmp2,tmp3,x1-tmp1,x2-tmp2,x3-tmp3);
+
 		dprintf(poso," %f %f %f %c %c %c\n",X1,X2,X3,t1,t2,t3);
 	}
 
-- 
2.20.1


From 705f50a9951d521843ecf30b71ec24703200f25f Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Wed, 26 Aug 2009 10:41:41 +0200
Subject: [PATCH 12/16] a more generic trafo tool

---
 vasp_tools/tXYZp.c | 285 +++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 285 insertions(+)
 create mode 100644 vasp_tools/tXYZp.c

diff --git a/vasp_tools/tXYZp.c b/vasp_tools/tXYZp.c
new file mode 100644
index 0000000..cd96040
--- /dev/null
+++ b/vasp_tools/tXYZp.c
@@ -0,0 +1,285 @@
+/*
+ * tXYZp.c
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <string.h>
+#include <unistd.h>
+#include <math.h>
+
+int get_line(int fd,char *line,int max) {
+
+        int count,ret;
+
+        count=0;
+
+        while(1) {
+                if(count==max) return count;
+                ret=read(fd,line+count,1);
+                if(ret<=0) return ret;
+                if(line[count]=='\n') {
+                        memset(line+count,0,max-count-1);
+                        //line[count]='\0';
+                        return count+1;
+                }
+                count+=1;
+        }
+}
+
+void usage(void) {
+	printf("usage: tXYZp -X,Y,Z <angle> -X,Y,Z <angle> ...\n");
+}
+
+int b_transform(char rtype,double angle,double cosangle,double sinangle,
+                double o[3][3],double t[3][3]) {
+
+	int i;
+
+	switch(rtype) {
+		case 'X':
+			for(i=0;i<3;i++) {
+				t[i][0]=o[i][0];
+				t[i][1]=cosangle*o[i][1]+sinangle*o[i][2];
+				t[i][2]=cosangle*o[i][2]-sinangle*o[i][1];
+			}
+			break;
+		case 'Y':
+			for(i=0;i<3;i++) {
+				t[i][0]=cosangle*o[i][0]+sinangle*o[i][2];
+				t[i][1]=o[i][1];
+				t[i][2]=cosangle*o[i][2]-sinangle*o[i][0];
+			}
+			break;
+		case 'Z':
+			for(i=0;i<3;i++) {
+				t[i][0]=cosangle*o[i][0]+sinangle*o[i][1];
+				t[i][1]=cosangle*o[i][1]-sinangle*o[i][0];
+				t[i][2]=o[i][2];
+			}
+			break;
+		default:
+			break;
+	}
+
+	printf("Transformed (%c - %f) basis:\n",rtype,angle);
+	printf(" b1: (%f, %f, %f)\n",t[0][0],t[0][1],t[0][2]);
+	printf(" b2: (%f, %f, %f)\n",t[1][0],t[1][1],t[1][2]);
+	printf(" b3: (%f, %f, %f)\n",t[2][0],t[2][1],t[2][2]);
+
+	return 1;
+}
+
+int a_transform(char rtype,double angle,double cosangle,double sinangle,
+                double o[3],double t[3]) {
+
+	switch(rtype) {
+		case 'X':
+			t[0]=o[0];
+			t[1]=cosangle*o[1]-sinangle*o[2];
+			t[2]=cosangle*o[2]+sinangle*o[1];
+			break;
+		case 'Y':
+			t[0]=cosangle*o[0]-sinangle*o[2];
+			t[1]=o[1];
+			t[2]=cosangle*o[2]+sinangle*o[0];
+			break;
+		case 'Z':
+			t[0]=cosangle*o[0]-sinangle*o[1];
+			t[1]=cosangle*o[1]+sinangle*o[0];
+			t[2]=o[2];
+			break;
+		default:
+			break;
+	}
+
+	return 1;
+}
+
+int main(int argc,char **argv) {
+
+	int i,j;
+	double angle[3];
+	char rtype[3];
+	double cosangle[3];
+	double sinangle[3];
+	int posr,poso;
+	char file[128];
+	int cnt,count;
+	char buf[256];
+	char *wptr;
+
+	/* basis in cartesian coordinates */
+	double b[3][3];
+
+	/* transformed basis */	
+	double tb[3][3][3];
+
+	double x[3],X[3];
+	char t1,t2,t3;
+
+	memset(angle,0,sizeof(double)*3);
+	memset(rtype,0,sizeof(char)*3);
+
+	count=0;
+	for(i=1;i<argc;i++) {
+		if(argv[i][0]=='-') {
+			switch(argv[i][1]) {
+				case 'X':
+				case 'Y':
+				case 'Z':
+					rtype[count]=argv[i][1];
+					angle[count]=atof(argv[++i])/180.0*M_PI;
+					count+=1;
+					break;
+				default:
+					usage();
+					return -1;
+			}
+		}
+		else {
+			usage();
+			return -1;
+		}
+	}
+
+	printf("\nOperations: ");
+	for(i=0;i<3;i++) {
+		if(rtype[i]) {
+			cosangle[i]=cos(angle[i]);
+			sinangle[i]=sin(angle[i]);
+			printf("%c %f (%d) | ",rtype[i],angle[i]/M_PI*180.0,i);
+		}
+	}
+	printf("\n\n");
+
+	posr=open("POSCAR",O_RDONLY);
+	if(posr<0) {
+		perror("open POSCAR (read) file\n");
+		return posr;
+	}
+
+	snprintf(file,255,"POSCAR.");
+	for(i=0;i<3;i++) {
+		if(rtype[i]) 
+			snprintf(file,255,"%s%c%f",file,rtype[i],angle[i]);
+	}
+	poso=open(file,O_WRONLY|O_CREAT);
+	if(poso<0) {
+		perror("open POSCAR (write) file\n");
+		return poso;
+	}
+	
+	// first line
+	cnt=get_line(posr,buf,256);
+	buf[cnt-1]='\n';
+	write(poso,buf,cnt);
+
+	// second line
+	cnt=get_line(posr,buf,256);
+	buf[cnt-1]='\n';
+	write(poso,buf,cnt);
+
+	// basis in cartesian coordinates
+	for(j=0;j<3;j++) {
+		cnt=get_line(posr,buf,256);
+		for(i=0;i<3;i++) {
+			if(i==0)
+				wptr=strtok(buf," ");
+			else
+				wptr=strtok(NULL," ");
+			b[j][i]=atof(wptr);
+		}
+	}
+
+	printf("Basis:\n");
+	printf(" b1: (%f, %f, %f)\n",b[0][0],b[0][1],b[0][2]);
+	printf(" b2: (%f, %f, %f)\n",b[1][0],b[1][1],b[1][2]);
+	printf(" b3: (%f, %f, %f)\n",b[2][0],b[2][1],b[2][2]);
+
+	/* transformation */
+	// remember:
+	// first rotation A
+	// second rotation B
+	// B needs to get transformed into B' = ABA^-1
+	// => total: ABA^-1A = AB
+	// => reverse sequence of transformations!
+	for(i=count-1;i>-1;i--) {
+		if(rtype[i]) {
+			if(i==count-1)
+				b_transform(rtype[i],angle[i],
+				            cosangle[i],sinangle[i],
+				            b,tb[i]);
+			else
+				b_transform(rtype[i],angle[i],
+				            cosangle[i],sinangle[i],
+				            tb[i+1],tb[i]);
+		}
+	}
+
+	for(i=0;i<3;i++)
+		dprintf(poso," %f %f %f\n",
+		        tb[0][i][0],tb[0][i][1],tb[0][i][2]);
+			
+	// 6th line
+	cnt=get_line(posr,buf,256);
+	buf[cnt-1]='\n';
+	write(poso,buf,cnt);
+
+	// 7th line
+	cnt=get_line(posr,buf,256);
+	buf[cnt-1]='\n';
+	write(poso,buf,cnt);
+
+	// 8th line
+	cnt=get_line(posr,buf,256);
+	buf[cnt-1]='\n';
+	write(poso,buf,cnt);
+
+	while(1) {
+		cnt=get_line(posr,buf,256);
+		if(cnt<=0)
+			break;
+		wptr=strtok(buf," ");
+		x[0]=atof(wptr);
+		wptr=strtok(NULL," ");
+		x[1]=atof(wptr);
+		wptr=strtok(NULL," ");
+		x[2]=atof(wptr);
+
+		wptr=strtok(NULL," ");
+		t1=wptr[0];
+		wptr=strtok(NULL," ");
+		t2=wptr[0];
+		wptr=strtok(NULL," ");
+		t3=wptr[0];
+
+		for(i=0;i<3;i++)  {
+			if(rtype[i]) {
+				X[0]=x[0];
+				X[1]=x[1];
+				X[2]=x[2];
+				a_transform(rtype[i],angle[i],
+				            cosangle[i],sinangle[i],
+				            X,x);
+			}
+		}
+
+		dprintf(poso," %f %f %f %c %c %c\n",x[0],x[1],x[2],t1,t2,t3);
+	}
+
+	close(posr);
+	close(poso);
+
+	printf("done!\n");
+
+	return 0;
+}
+
-- 
2.20.1


From 31a82cd2a10be5c6c7095c784588bb43c48a4c3d Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Wed, 26 Aug 2009 16:52:36 +0200
Subject: [PATCH 13/16] care for appendices

---
 vasp_tools/poscar2moldyn | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/vasp_tools/poscar2moldyn b/vasp_tools/poscar2moldyn
index b0ad584..a4f2562 100755
--- a/vasp_tools/poscar2moldyn
+++ b/vasp_tools/poscar2moldyn
@@ -2,8 +2,13 @@
 
 mkdir -p video
 
-file=$1
-[ -z $1 ] && file=POSCAR
+if [ -z $1 ]; then
+	file=POSCAR
+	app=""
+else
+	file=$1
+	app=".`echo $file | awk -F. '{ print $2 }'`"
+fi
 
 echo "parsing $file file ..."
 
@@ -25,6 +30,7 @@ z2=`sed -n 5p $file | awk '{ print $2 }'`
 z3=`sed -n 5p $file | awk '{ print $3 }'`
 
 ((total=sicnt+ccnt))
+((eoa=total+8))
 
 echo "  Si: $sicnt, C: $ccnt, total: $total"
 echo "  Lattice constant: $lc A"
@@ -40,9 +46,9 @@ cx=1.0
 cy=1.0
 cz=0.8
 
-echo "# P $total init <$cx,$cy,$cz>" > video/atomic_conf_init.xyz
+echo "# P $total init <$cx,$cy,$cz>" > video/atomic_conf_init${app}.xyz
 
-tail -${total} $file | \
+sed -n 9,${eoa}p $file | \
 	while read x y z fx fy fz; do
 		type="Si"
 		color="Yellow"
@@ -59,7 +65,7 @@ tail -${total} $file | \
 		X=`echo $lc $X | awk '{ print $1*$2 }'`
 		Y=`echo $lc $Y | awk '{ print $1*$2 }'`
 		Z=`echo $lc $Z | awk '{ print $1*$2 }'`
-		echo "$type $X $Y $Z $color 0.0" >> video/atomic_conf_init.xyz
+		echo "$type $X $Y $Z $color 0.0" >> video/atomic_conf_init${app}.xyz
 		((cnt+=1))
 	done
 
-- 
2.20.1


From 2d50cbdcdbeaa83ee5a538ab5e6114b643679174 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Thu, 10 Sep 2009 08:34:25 +0200
Subject: [PATCH 14/16] total energies printed out too

---
 vasp_tools/e_coh | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/vasp_tools/e_coh b/vasp_tools/e_coh
index 6d8a107..44adfd9 100755
--- a/vasp_tools/e_coh
+++ b/vasp_tools/e_coh
@@ -65,9 +65,9 @@ echo "  Si correction: $free_si eV, C correction: $free_c eV"
 echo $energy $total | \
 	awk '{ print "  total e: " $1 " eV, per atom: " $1/$2 " eV"}'
 echo "$energy $sicnt $ccnt $free_si $free_c" | \
-	awk '{ print "  cohesive energy (Si and C): " ($1-$2*$4-$3*$5)/($2+$3) " eV" }'
+	awk '{ print "  cohesive energy (Si and C): " ($1-$2*$4-$3*$5)/($2+$3) " eV, " $1-$2*$4-$3*$5 " eV" }'
 echo "$energy $sicnt $free_c" | \
-	awk '{ print "  cohesive energy (C only): " ($1-$2*$3)/$2 " eV" }'
+	awk '{ print "  cohesive energy (C only): " ($1-$2*$3)/$2 " eV, " $1-$2*$3 " eV" }'
 #echo "$energy $total $free_si_250" | \
 #	awk '{ print "  cohesive energy (Si only 250): " ($1-$2*$3)/$2 " eV" }'
 #echo "$energy $sicnt $free_c_650" | \
-- 
2.20.1


From cba0d347201db1c99dc23736658ed78ed12ad5bd Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Fri, 11 Sep 2009 18:31:39 +0200
Subject: [PATCH 15/16] added msd calc + 11x constraints as ifdef

---
 moldyn.c   | 151 +++++++++++++++++++++++++++++++++++++++++++++++++++++
 moldyn.h   |   1 +
 msd_calc.c |  54 +++++++++++++++++++
 3 files changed, 206 insertions(+)
 create mode 100644 msd_calc.c

diff --git a/moldyn.c b/moldyn.c
index 15dcede..bd29afc 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -81,6 +81,10 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 #ifdef CONSTRAINT_110_5832
 	printf("\n\n\nWARNING! WARNING! WARNING!\n\n\n");
 	printf("\n\n\n!! -- constraints enabled -- !!\n\n\n");
+#endif
+#ifdef CONSTRAINT_11X_5832
+	printf("\n\n\nWARNING! WARNING! WARNING!\n\n\n");
+	printf("\n\n\n!! -- constraints enabled -- !!\n\n\n");
 #endif
 	return 0;
 }
@@ -2264,6 +2268,10 @@ int velocity_verlet(t_moldyn *moldyn) {
 	double tau,tau_square,h;
 	t_3dvec delta;
 	t_atom *atom;
+#ifdef CONSTRAINT_11X_5832
+	double xt,yt,zt;
+	double xtt,ytt,ztt;
+#endif
 
 	atom=moldyn->atom;
 	count=moldyn->count;
@@ -2275,6 +2283,28 @@ int velocity_verlet(t_moldyn *moldyn) {
 		atom[5832].f.x=0.5*(atom[5832].f.x-atom[5832].f.y);
 		atom[5832].f.y=-atom[5832].f.x;
 	}
+#endif
+#ifdef CONSTRAINT_11X_5832
+	if(count==5833) {
+		// second trafo
+		xt=atom[5832].f.x;
+		yt=atom[5832].f.y*cos(-0.16935129)-atom[5832].f.z*sin(-0.16935129);
+		zt=atom[5832].f.y*sin(-0.16935129)+atom[5832].f.z*cos(-0.16935129);
+		// first trafo
+		xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
+		ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
+		ztt=zt;
+		// apply constraints
+		ytt=0.0;
+		// first trafo backwards
+		xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
+		yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
+		zt=ztt;
+		// second trafo backwards
+		atom[5832].f.x=xt;
+		atom[5832].f.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
+		atom[5832].f.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
+	}
 #endif
 	for(i=0;i<count;i++) {
 		/* check whether fixed atom */
@@ -2288,6 +2318,30 @@ int velocity_verlet(t_moldyn *moldyn) {
 			delta.y=-delta.x;
 		}
 #endif
+#ifdef JHDVHJDV
+#ifdef CONSTRAINT_11X_5832
+	if(i==5833) {
+		// second trafo
+		xt=delta.x;
+		yt=delta.y*cos(-0.16935129)-delta.z*sin(-0.16935129);
+		zt=delta.y*sin(-0.16935129)+delta.z*cos(-0.16935129);
+		// first trafo
+		xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
+		ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
+		ztt=zt;
+		// apply constraints
+		ytt=0.0;
+		// first trafo backwards
+		xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
+		yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
+		zt=ztt;
+		// second trafo backwards
+		delta.x=xt;
+		delta.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
+		delta.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
+	}
+#endif
+#endif
 #ifndef QUENCH
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
 #endif
@@ -2296,6 +2350,30 @@ int velocity_verlet(t_moldyn *moldyn) {
 		if(i==5832) {
 			delta.y=-delta.x;
 		}
+#endif
+#ifdef GHDJHBSJBSD
+#ifdef CONSTRAINT_11X_5832
+	if(i==5833) {
+		// second trafo
+		xt=delta.x;
+		yt=delta.y*cos(-0.16935129)-delta.z*sin(-0.16935129);
+		zt=delta.y*sin(-0.16935129)+delta.z*cos(-0.16935129);
+		// first trafo
+		xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
+		ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
+		ztt=zt;
+		// apply constraints
+		ytt=0.0;
+		// first trafo backwards
+		xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
+		yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
+		zt=ztt;
+		// second trafo backwards
+		delta.x=xt;
+		delta.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
+		delta.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
+	}
+#endif
 #endif
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
 		//check_per_bound_and_care_for_pbc(moldyn,&(atom[i]));
@@ -2328,6 +2406,28 @@ int velocity_verlet(t_moldyn *moldyn) {
 		atom[5832].f.x=0.5*(atom[5832].f.x-atom[5832].f.y);
 		atom[5832].f.y=-atom[5832].f.x;
 	}
+#endif
+#ifdef CONSTRAINT_11X_5832
+	if(count==5833) {
+		// second trafo
+		xt=atom[5832].f.x;
+		yt=atom[5832].f.y*cos(-0.16935129)-atom[5832].f.z*sin(-0.16935129);
+		zt=atom[5832].f.y*sin(-0.16935129)+atom[5832].f.z*cos(-0.16935129);
+		// first trafo
+		xtt=xt*cos(M_PI/4.0)-yt*sin(M_PI/4.0);
+		ytt=xt*sin(M_PI/4.0)+yt*sin(M_PI/4.0);
+		ztt=zt;
+		// apply constraints
+		ytt=0.0;
+		// first trafo backwards
+		xt=xtt*cos(M_PI/4.0)+ytt*sin(M_PI/4.0);
+		yt=-xtt*sin(M_PI/4.0)+ytt*sin(M_PI/4.0);
+		zt=ztt;
+		// second trafo backwards
+		atom[5832].f.x=xt;
+		atom[5832].f.y=yt*cos(-0.16935129)+zt*sin(-0.16935129);
+		atom[5832].f.z=-yt*sin(-0.16935129)+zt*cos(-0.16935129);
+	}
 #endif
 	for(i=0;i<count;i++) {
 		/* check whether fixed atom */
@@ -3197,6 +3297,57 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 	return 0;
 }
 
+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;
+}
+
 int bonding_analyze(t_moldyn *moldyn,double *cnt) {
 
 	return 0;
diff --git a/moldyn.h b/moldyn.h
index bdb3758..a1866a9 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -523,6 +523,7 @@ int get_line(int fd,char *line,int max);
 
 int pair_correlation_init(t_moldyn *moldyn,double dr);
 int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc);
+int calculate_msd(t_moldyn *moldyn,double *msd);
 int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
                                        t_atom *jtom,void *data,u8 bc);
 int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr);
diff --git a/msd_calc.c b/msd_calc.c
new file mode 100644
index 0000000..80388d8
--- /dev/null
+++ b/msd_calc.c
@@ -0,0 +1,54 @@
+/*
+ * calculation of mean square displacement
+ *
+ * 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 msd[3];
+
+	if(argc!=2) {
+		usage(argv[0]);
+		return -1;
+	}
+
+	memset(&moldyn,0,sizeof(t_moldyn));
+
+	printf("[msd calc] reading save file ...\n");
+	ret=moldyn_read_save_file(&moldyn,argv[1]);
+	if(ret) {
+		printf("[msd calc] exit!\n");
+		return ret;
+	}
+
+	calculate_msd(&moldyn,msd);
+	
+	printf("MSD - %f ps: %.10f %.10f %.10f\n",moldyn.time,
+	       msd[0],msd[1],msd[2]);
+
+	return 0;
+}
+
-- 
2.20.1


From 0cc57bba82ec0e4208f48caf793726050a30b998 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Sun, 13 Sep 2009 13:20:49 +0200
Subject: [PATCH 16/16] added a pc calc for ascci files produced by
 outcar2moldyn ...

---
 vasp_tools/pc_calc.c | 190 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 190 insertions(+)
 create mode 100644 vasp_tools/pc_calc.c

diff --git a/vasp_tools/pc_calc.c b/vasp_tools/pc_calc.c
new file mode 100644
index 0000000..10275d2
--- /dev/null
+++ b/vasp_tools/pc_calc.c
@@ -0,0 +1,190 @@
+/*
+ * pc_calc.c - calculate the pc
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+
+#include <math.h>
+
+
+#define MAXR	6.0
+#define DELTA	0.01
+
+typedef struct s_atom {
+	double x,y,z;
+	char type;
+} t_atom;
+
+int get_line(int fd,char *line,int max) {
+
+        int count,ret;
+
+        count=0;
+
+        while(1) {
+                if(count==max) return count;
+                ret=read(fd,line+count,1);
+                if(ret<=0) return ret;
+                if(line[count]=='\n') {
+                        memset(line+count,0,max-count-1);
+                        //line[count]='\0';
+                        return count+1;
+                }
+                count+=1;
+        }
+}
+
+int main(int argc,char **argv) {
+
+	t_atom *atom;
+	double *aslot,*bslot,*cslot;
+	int acnt,bcnt,ccnt,cnt,count,slots;
+	int fd;
+	char buf[256],*wptr;
+	int i,j,k;
+	double dx2,dy2,dz2,dist,norm;
+	double sx,sy,sz;
+
+	if(argc!=5) {
+		printf("usage: %s file sx sy sz\n",argv[0]);
+		return -1;
+	}
+
+	fd=open(argv[1],O_RDONLY);
+	if(fd<0) {
+		perror("open file\n");
+		return fd;
+	}
+
+	sx=atof(argv[2]);
+	sy=atof(argv[3]);
+	sz=atof(argv[4]);
+
+	// first line
+	cnt=get_line(fd,buf,256);
+
+	// read in atoms
+	count=0;
+	while(1) {
+		cnt=get_line(fd,buf,256);
+		if(cnt<=0)
+			break;
+
+		// there is (another) atom	
+
+		atom=realloc(atom,(count+1)*sizeof(t_atom));
+		wptr=strtok(buf," ");
+		atom[count].type=wptr[0];
+
+		wptr=strtok(NULL," ");
+		atom[count].x=atof(wptr);
+		wptr=strtok(NULL," ");
+		atom[count].y=atof(wptr);
+		wptr=strtok(NULL," ");
+		atom[count].z=atof(wptr);
+
+		count+=1;
+	}
+
+	printf("i: read in %d atoms ...\n",count);
+
+	// prepare pc
+
+	slots=MAXR/DELTA;
+	aslot=malloc(slots*sizeof(double));
+	if(aslot==NULL) {
+		perror("slot a\n");
+		return -1;
+	}
+	memset(aslot,0,slots*sizeof(double));
+	bslot=malloc(slots*sizeof(double));
+	if(bslot==NULL) {
+		perror("slot a\n");
+		return -1;
+	}
+	memset(bslot,0,slots*sizeof(double));
+	cslot=malloc(slots*sizeof(double));
+	if(cslot==NULL) {
+		perror("slot a\n");
+		return -1;
+	}
+	memset(cslot,0,slots*sizeof(double));
+
+	printf("i: allocated 3 times %d slots ...\n",slots);
+
+	acnt=0;
+	bcnt=0;
+	ccnt=0;
+
+	// calc pc
+
+	for(i=0;i<count;i++) {
+		for(j=0;j<i;j++) {
+			dx2=atom[i].x-atom[j].x;
+			if(dx2>sx/2.0)
+				dx2-=sx;
+			else if(dx2<-sx/2.0)
+				dx2+=sx;
+			dx2*=dx2;
+			dy2=atom[i].y-atom[j].y;
+			if(dy2>sy/2.0)
+				dy2-=sy;
+			else if(dy2<-sy/2.0)
+				dy2+=sy;
+			dy2*=dy2;
+			dz2=atom[i].z-atom[j].z;
+			if(dz2>sz/2.0)
+				dz2-=sz;
+			else if(dz2<-sz/2.0)
+				dz2+=sz;
+			dz2*=dz2;
+			dist=sqrt(dx2+dy2+dz2);
+			if(dist>=MAXR)
+				continue;
+			k=dist/DELTA;
+			if((atom[i].type=='S')&&(atom[j].type=='S')) {
+				aslot[k]+=1;
+				acnt+=1;
+			}
+			else if((atom[i].type=='C')&&(atom[j].type=='C')) {
+				bslot[k]+=1;
+				bcnt+=1;
+			}
+			else {
+				cslot[k]+=1;
+				ccnt+=1;
+			}
+		}
+	}
+
+	// normalization and output
+
+	for(i=1;i<slots;i++) {
+		norm=4*M_PI*(i*DELTA)*(i*DELTA)*DELTA;
+		printf("pc: %f %f %f %f\n",i*DELTA,
+		                           aslot[i]/norm,
+		                           bslot[i]/norm,
+		                           cslot[i]/norm);
+	}
+
+	free(aslot);
+	free(bslot);
+	free(cslot);
+
+	free(atom);
+
+	close(fd);
+
+	return 0;
+
+}
+
-- 
2.20.1