a more generic trafo tool

diff --git a/vasp_tools/tXYZp.c b/vasp_tools/tXYZp.c
new file mode 100644 (file)
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;
+}
+