]> www.hackdaworld.org Git - physik/posic.git/blob - vasp_tools/pc_calc.c
added fractional output (search bonds) + video{.*} feature (stdvis)
[physik/posic.git] / vasp_tools / pc_calc.c
1 /*
2  * pc_calc.c - calculate the pc
3  *
4  * author: frank.zirkelbach@physik.uni-augsburg.de
5  *
6  */
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <string.h>
11 #include <sys/types.h>
12 #include <sys/stat.h>
13 #include <fcntl.h>
14 #include <unistd.h>
15
16 #include <math.h>
17
18
19 #define MAXR    6.0
20 #define DELTA   0.01
21
22 typedef struct s_atom {
23         double x,y,z;
24         char type;
25 } t_atom;
26
27 int get_line(int fd,char *line,int max) {
28
29         int count,ret;
30
31         count=0;
32
33         while(1) {
34                 if(count==max) return count;
35                 ret=read(fd,line+count,1);
36                 if(ret<=0) return ret;
37                 if(line[count]=='\n') {
38                         memset(line+count,0,max-count-1);
39                         //line[count]='\0';
40                         return count+1;
41                 }
42                 count+=1;
43         }
44 }
45
46 int get_max3_vals(char *buf,int amount,double *v1, double *v2,double *v3) {
47
48         int i,j,k;
49         char temp[128];
50
51         i=0;
52
53         for(k=0;k<amount;k++) {
54                 while(buf[i]==' ')
55                         i++;
56                 j=0;
57                 while((buf[i+j]!=' ')&&(buf[i+j]!=0))
58                         j++;
59                 memcpy(temp,buf+i,j);
60                 temp[j]=0;
61                 if(k==0) *v1=atof(temp);
62                 if(k==1) *v2=atof(temp);
63                 if(k==2) *v3=atof(temp);
64                 i+=j;
65         }
66
67         return 0;
68 }
69
70 int get_2_ints(char *buf,int *i1,int *i2) {
71
72         int i,j,k;
73         char temp[128];
74
75         i=0;
76
77         for(k=0;k<2;k++) {
78                 while(buf[i]==' ')
79                         i++;
80                 j=0;
81                 while((buf[i+j]!=' ')&&(buf[i+j]!=0))
82                         j++;
83                 memcpy(temp,buf+i,j);
84                 temp[j]=0;
85                 if(k==0) *i1=atoi(temp);
86                 if(k==1) *i2=atoi(temp);
87                 i+=j;
88         }
89
90         return 0;
91 }
92
93 int main(int argc,char **argv) {
94
95         t_atom *atom;
96         double *aslot,*bslot,*cslot;
97         int cnt,count,fcnt,slots;
98         int fd;
99         char buf[256];
100         int i,j,k;
101         int nsi,nc,ntot;
102         double dx,dy,dz,dist,norm;
103         double dxt,dyt,dzt;
104         double X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3;
105         double scale;
106
107         atom=NULL;
108
109         if(argc<1) {
110                 printf("usage: %s file1 file2 ...\n",
111                        argv[0]);
112                 return -1;
113         }
114
115         // prepare pc
116         slots=MAXR/DELTA;
117         aslot=malloc(slots*sizeof(double));
118         if(aslot==NULL) {
119                 perror("slot a\n");
120                 return -1;
121         }
122         memset(aslot,0,slots*sizeof(double));
123         bslot=malloc(slots*sizeof(double));
124         if(bslot==NULL) {
125                 perror("slot a\n");
126                 return -1;
127         }
128         memset(bslot,0,slots*sizeof(double));
129         cslot=malloc(slots*sizeof(double));
130         if(cslot==NULL) {
131                 perror("slot a\n");
132                 return -1;
133         }
134         memset(cslot,0,slots*sizeof(double));
135         printf("i: allocated 3 times %d slots ...\n",slots);
136
137         // use all given files ...
138         printf("using files:\n");
139         for(fcnt=1;fcnt<argc;fcnt++)
140                 printf(" %d: %s\n",fcnt,argv[fcnt]);
141
142         for(fcnt=1;fcnt<argc;fcnt++) {
143
144         fd=open(argv[fcnt],O_RDONLY);
145         if(fd<0) {
146                 perror("open file\n");
147                 return fd;
148         }
149
150         // first line
151         cnt=get_line(fd,buf,256);
152
153         // second line -> scale
154         cnt=get_line(fd,buf,256);
155         get_max3_vals(buf,1,&scale,NULL,NULL);
156         printf("scale: %f\n",scale);
157
158         // third line -> X
159         cnt=get_line(fd,buf,256);
160         get_max3_vals(buf,3,&X1,&X2,&X3);
161         printf("X: %f %f %f\n",X1,X2,X3);
162         
163         // 4th line -> Y
164         cnt=get_line(fd,buf,256);
165         get_max3_vals(buf,3,&Y1,&Y2,&Y3);
166         printf("Y: %f %f %f\n",Y1,Y2,Y3);
167         
168         // 5th line -> Z
169         cnt=get_line(fd,buf,256);
170         get_max3_vals(buf,3,&Z1,&Z2,&Z3);
171         printf("Z: %f %f %f\n",Z1,Z2,Z3);
172         
173         // 6th line -> nsi + nc = ntot
174         cnt=get_line(fd,buf,256);
175         get_2_ints(buf,&nsi,&nc);
176         ntot=nsi+nc;
177         printf("Si: %d - C: %d - tot: %d\n",nsi,nc,ntot);
178         
179         // 7th line
180         cnt=get_line(fd,buf,256);
181
182         // 8th line
183         cnt=get_line(fd,buf,256);
184
185         // read in atoms
186         count=0;
187         while(count<ntot) {
188                 cnt=get_line(fd,buf,256);
189                 if(cnt<=0) {
190                         printf("Should never happen!\n");
191                         return cnt;
192                 }
193
194                 atom=realloc(atom,(count+1)*sizeof(t_atom));
195                 if(count<nsi)
196                         atom[count].type='S';
197                 else
198                         atom[count].type='C';
199
200                 get_max3_vals(buf,3,
201                               &(atom[count].x),
202                               &(atom[count].y),
203                               &(atom[count].z));
204
205                 count+=1;
206         }
207
208         printf("i: read in %d atoms ...\n",count);
209
210         // calc pc
211         for(i=0;i<count;i++) {
212                 for(j=0;j<i;j++) {
213                         dx=atom[i].x-atom[j].x;
214                         if(dx>0.5)
215                                 dx-=1;
216                         else if(dx<-0.5)
217                                 dx+=1;
218
219                         dy=atom[i].y-atom[j].y;
220                         if(dy>0.5)
221                                 dy-=1;
222                         else if(dy<-0.5)
223                                 dy+=1;
224
225                         dz=atom[i].z-atom[j].z;
226                         if(dz>0.5)
227                                 dz-=1;
228                         else if(dz<-0.5)
229                                 dz+=1;
230
231                         dxt=dx*X1+dy*Y1+dz*Z1;
232                         dyt=dx*X2+dy*Y2+dz*Z2;
233                         dzt=dx*X3+dy*Y3+dz*Z3;
234
235                         dist=sqrt(dxt*dxt+dyt*dyt+dzt*dzt);
236                         dist*=scale;
237
238                         if(dist>=MAXR)
239                                 continue;
240                         k=dist/DELTA;
241                         if((atom[i].type=='S')&&(atom[j].type=='S')) {
242                                 aslot[k]+=1;
243                         }
244                         else if((atom[i].type=='C')&&(atom[j].type=='C')) {
245                                 bslot[k]+=1;
246                         }
247                         else {
248                                 cslot[k]+=1;
249                         }
250                 }
251         }
252
253         close(fd);
254
255         }
256
257         // normalization and output
258         for(i=1;i<slots;i++) {
259                 norm=4*M_PI*(i*DELTA)*(i*DELTA)*DELTA;
260                 norm*=(argc-1);
261                 printf("pc: %f %f %f %f\n",i*DELTA,
262                                            aslot[i]/norm,
263                                            bslot[i]/norm,
264                                            cslot[i]/norm);
265         }
266
267         free(aslot);
268         free(bslot);
269         free(cslot);
270
271         free(atom);
272
273         return 0;
274
275 }
276