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] 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.39.2