]> hackdaworld.org Git - physik/posic.git/blob - vasp_tools/tXYZp.c
count spins ...
[physik/posic.git] / 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