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