vasp_tools/tXYZp.c

   1 /*
   2  * tXYZp.c
   3  *
   4  * author: frank.zirkelbach@physik.uni-augsburg.de
   5  *
   6  */
   7
   8 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <sys/types.h>
  12 #include <sys/stat.h>
  13 #include <fcntl.h>
  14 #include <string.h>
  15 #include <unistd.h>
  16 #include <math.h>
  17
  18 int get_line(int fd,char *line,int max) {
  19
  20         int count,ret;
  21
  22         count=0;
  23
  24         while(1) {
  25                 if(count==max) return count;
  26                 ret=read(fd,line+count,1);
  27                 if(ret<=0) return ret;
  28                 if(line[count]=='\n') {
  29                         memset(line+count,0,max-count-1);
  30                         //line[count]='\0';
  31                         return count+1;
  32                 }
  33                 count+=1;
  34         }
  35 }
  36
  37 void usage(void) {
  38         printf("usage: tXYZp -X,Y,Z <angle> -X,Y,Z <angle> ...\n");
  39 }
  40
  41 int b_transform(char rtype,double angle,double cosangle,double sinangle,
  42                 double o[3][3],double t[3][3]) {
  43
  44         int i;
  45
  46         switch(rtype) {
  47                 case 'X':
  48                         for(i=0;i<3;i++) {
  49                                 t[i][0]=o[i][0];
  50                                 t[i][1]=cosangle*o[i][1]+sinangle*o[i][2];
  51                                 t[i][2]=cosangle*o[i][2]-sinangle*o[i][1];
  52                         }
  53                         break;
  54                 case 'Y':
  55                         for(i=0;i<3;i++) {
  56                                 t[i][0]=cosangle*o[i][0]+sinangle*o[i][2];
  57                                 t[i][1]=o[i][1];
  58                                 t[i][2]=cosangle*o[i][2]-sinangle*o[i][0];
  59                         }
  60                         break;
  61                 case 'Z':
  62                         for(i=0;i<3;i++) {
  63                                 t[i][0]=cosangle*o[i][0]+sinangle*o[i][1];
  64                                 t[i][1]=cosangle*o[i][1]-sinangle*o[i][0];
  65                                 t[i][2]=o[i][2];
  66                         }
  67                         break;
  68                 default:
  69                         break;
  70         }
  71
  72         printf("Transformed (%c - %f) basis:\n",rtype,angle);
  73         printf(" b1: (%f, %f, %f)\n",t[0][0],t[0][1],t[0][2]);
  74         printf(" b2: (%f, %f, %f)\n",t[1][0],t[1][1],t[1][2]);
  75         printf(" b3: (%f, %f, %f)\n",t[2][0],t[2][1],t[2][2]);
  76
  77         return 1;
  78 }
  79
  80 int a_transform(char rtype,double angle,double cosangle,double sinangle,
  81                 double o[3],double t[3]) {
  82
  83         switch(rtype) {
  84                 case 'X':
  85                         t[0]=o[0];
  86                         t[1]=cosangle*o[1]-sinangle*o[2];
  87                         t[2]=cosangle*o[2]+sinangle*o[1];
  88                         break;
  89                 case 'Y':
  90                         t[0]=cosangle*o[0]-sinangle*o[2];
  91                         t[1]=o[1];
  92                         t[2]=cosangle*o[2]+sinangle*o[0];
  93                         break;
  94                 case 'Z':
  95                         t[0]=cosangle*o[0]-sinangle*o[1];
  96                         t[1]=cosangle*o[1]+sinangle*o[0];
  97                         t[2]=o[2];
  98                         break;
  99                 default:
 100                         break;
 101         }
 102
 103         return 1;
 104 }
 105
 106 int main(int argc,char **argv) {
 107
 108         int i,j;
 109         double angle[3];
 110         char rtype[3];
 111         double cosangle[3];
 112         double sinangle[3];
 113         int posr,poso;
 114         char file[128];
 115         int cnt,count;
 116         char buf[256];
 117         char *wptr;
 118
 119         /* basis in cartesian coordinates */
 120         double b[3][3];
 121
 122         /* transformed basis */
 123         double tb[3][3][3];
 124
 125         double x[3],X[3];
 126         char t1,t2,t3;
 127
 128         memset(angle,0,sizeof(double)*3);
 129         memset(rtype,0,sizeof(char)*3);
 130
 131         count=0;
 132         for(i=1;i<argc;i++) {
 133                 if(argv[i][0]=='-') {
 134                         switch(argv[i][1]) {
 135                                 case 'X':
 136                                 case 'Y':
 137                                 case 'Z':
 138                                         rtype[count]=argv[i][1];
 139                                         angle[count]=atof(argv[++i])/180.0*M_PI;
 140                                         count+=1;
 141                                         break;
 142                                 default:
 143                                         usage();
 144                                         return -1;
 145                         }
 146                 }
 147                 else {
 148                         usage();
 149                         return -1;
 150                 }
 151         }
 152
 153         printf("\nOperations: ");
 154         for(i=0;i<3;i++) {
 155                 if(rtype[i]) {
 156                         cosangle[i]=cos(angle[i]);
 157                         sinangle[i]=sin(angle[i]);
 158                         printf("%c %f (%d) | ",rtype[i],angle[i]/M_PI*180.0,i);
 159                 }
 160         }
 161         printf("\n\n");
 162
 163         posr=open("POSCAR",O_RDONLY);
 164         if(posr<0) {
 165                 perror("open POSCAR (read) file\n");
 166                 return posr;
 167         }
 168
 169         snprintf(file,255,"POSCAR.");
 170         for(i=0;i<3;i++) {
 171                 if(rtype[i])
 172                         snprintf(file,255,"%s%c%f",file,rtype[i],angle[i]);
 173         }
 174         poso=open(file,O_WRONLY|O_CREAT);
 175         if(poso<0) {
 176                 perror("open POSCAR (write) file\n");
 177                 return poso;
 178         }
 179
 180         // first line
 181         cnt=get_line(posr,buf,256);
 182         buf[cnt-1]='\n';
 183         write(poso,buf,cnt);
 184
 185         // second line
 186         cnt=get_line(posr,buf,256);
 187         buf[cnt-1]='\n';
 188         write(poso,buf,cnt);
 189
 190         // basis in cartesian coordinates
 191         for(j=0;j<3;j++) {
 192                 cnt=get_line(posr,buf,256);
 193                 for(i=0;i<3;i++) {
 194                         if(i==0)
 195                                 wptr=strtok(buf," ");
 196                         else
 197                                 wptr=strtok(NULL," ");
 198                         b[j][i]=atof(wptr);
 199                 }
 200         }
 201
 202         printf("Basis:\n");
 203         printf(" b1: (%f, %f, %f)\n",b[0][0],b[0][1],b[0][2]);
 204         printf(" b2: (%f, %f, %f)\n",b[1][0],b[1][1],b[1][2]);
 205         printf(" b3: (%f, %f, %f)\n",b[2][0],b[2][1],b[2][2]);
 206
 207         /* transformation */
 208         // remember:
 209         // first rotation A
 210         // second rotation B
 211         // B needs to get transformed into B' = ABA^-1
 212         // => total: ABA^-1A = AB
 213         // => reverse sequence of transformations!
 214         for(i=count-1;i>-1;i--) {
 215                 if(rtype[i]) {
 216                         if(i==count-1)
 217                                 b_transform(rtype[i],angle[i],
 218                                             cosangle[i],sinangle[i],
 219                                             b,tb[i]);
 220                         else
 221                                 b_transform(rtype[i],angle[i],
 222                                             cosangle[i],sinangle[i],
 223                                             tb[i+1],tb[i]);
 224                 }
 225         }
 226
 227         for(i=0;i<3;i++)
 228                 dprintf(poso," %f %f %f\n",
 229                         tb[0][i][0],tb[0][i][1],tb[0][i][2]);
 230
 231         // 6th line
 232         cnt=get_line(posr,buf,256);
 233         buf[cnt-1]='\n';
 234         write(poso,buf,cnt);
 235
 236         // 7th line
 237         cnt=get_line(posr,buf,256);
 238         buf[cnt-1]='\n';
 239         write(poso,buf,cnt);
 240
 241         // 8th line
 242         cnt=get_line(posr,buf,256);
 243         buf[cnt-1]='\n';
 244         write(poso,buf,cnt);
 245
 246         while(1) {
 247                 cnt=get_line(posr,buf,256);
 248                 if(cnt<=0)
 249                         break;
 250                 wptr=strtok(buf," ");
 251                 x[0]=atof(wptr);
 252                 wptr=strtok(NULL," ");
 253                 x[1]=atof(wptr);
 254                 wptr=strtok(NULL," ");
 255                 x[2]=atof(wptr);
 256
 257                 wptr=strtok(NULL," ");
 258                 t1=wptr[0];
 259                 wptr=strtok(NULL," ");
 260                 t2=wptr[0];
 261                 wptr=strtok(NULL," ");
 262                 t3=wptr[0];
 263
 264                 for(i=0;i<3;i++)  {
 265                         if(rtype[i]) {
 266                                 X[0]=x[0];
 267                                 X[1]=x[1];
 268                                 X[2]=x[2];
 269                                 a_transform(rtype[i],angle[i],
 270                                             cosangle[i],sinangle[i],
 271                                             X,x);
 272                         }
 273                 }
 274
 275                 dprintf(poso," %f %f %f %c %c %c\n",x[0],x[1],x[2],t1,t2,t3);
 276         }
 277
 278         close(posr);
 279         close(poso);
 280
 281         printf("done!\n");
 282
 283         return 0;
 284 }
 285