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 main(int argc,char **argv) {
  47
  48         t_atom *atom;
  49         double *aslot,*bslot,*cslot;
  50         int cnt,count,fcnt,slots;
  51         int fd;
  52         char buf[256],*wptr;
  53         int i,j,k;
  54         double dx2,dy2,dz2,dist,norm;
  55         double sx,sy,sz;
  56
  57         atom=NULL;
  58
  59         if(argc<5) {
  60                 printf("usage: %s sx sy sz file1 file2 ...\n",argv[0]);
  61                 return -1;
  62         }
  63
  64         sx=atof(argv[1]);
  65         sy=atof(argv[2]);
  66         sz=atof(argv[3]);
  67
  68         // prepare pc
  69
  70         slots=MAXR/DELTA;
  71         aslot=malloc(slots*sizeof(double));
  72         if(aslot==NULL) {
  73                 perror("slot a\n");
  74                 return -1;
  75         }
  76         memset(aslot,0,slots*sizeof(double));
  77         bslot=malloc(slots*sizeof(double));
  78         if(bslot==NULL) {
  79                 perror("slot a\n");
  80                 return -1;
  81         }
  82         memset(bslot,0,slots*sizeof(double));
  83         cslot=malloc(slots*sizeof(double));
  84         if(cslot==NULL) {
  85                 perror("slot a\n");
  86                 return -1;
  87         }
  88         memset(cslot,0,slots*sizeof(double));
  89
  90         printf("i: allocated 3 times %d slots ...\n",slots);
  91
  92         // use all given files ...
  93         printf("using files:\n");
  94         for(fcnt=4;fcnt<argc;fcnt++)
  95                 printf(" %d: %s\n",fcnt-4,argv[fcnt]);
  96
  97         for(fcnt=4;fcnt<argc;fcnt++) {
  98
  99         fd=open(argv[fcnt],O_RDONLY);
 100         if(fd<0) {
 101                 perror("open file\n");
 102                 return fd;
 103         }
 104
 105         // first line
 106         cnt=get_line(fd,buf,256);
 107
 108         // read in atoms
 109         count=0;
 110         while(1) {
 111                 cnt=get_line(fd,buf,256);
 112                 if(cnt<=0)
 113                         break;
 114
 115                 // there is (another) atom
 116
 117                 atom=realloc(atom,(count+1)*sizeof(t_atom));
 118                 wptr=strtok(buf," ");
 119                 atom[count].type=wptr[0];
 120
 121                 wptr=strtok(NULL," ");
 122                 atom[count].x=atof(wptr);
 123                 wptr=strtok(NULL," ");
 124                 atom[count].y=atof(wptr);
 125                 wptr=strtok(NULL," ");
 126                 atom[count].z=atof(wptr);
 127
 128                 count+=1;
 129         }
 130
 131         printf("i: read in %d atoms ...\n",count);
 132
 133         // calc pc
 134
 135         for(i=0;i<count;i++) {
 136                 for(j=0;j<i;j++) {
 137                         dx2=atom[i].x-atom[j].x;
 138                         if(dx2>sx/2.0)
 139                                 dx2-=sx;
 140                         else if(dx2<-sx/2.0)
 141                                 dx2+=sx;
 142                         dx2*=dx2;
 143                         dy2=atom[i].y-atom[j].y;
 144                         if(dy2>sy/2.0)
 145                                 dy2-=sy;
 146                         else if(dy2<-sy/2.0)
 147                                 dy2+=sy;
 148                         dy2*=dy2;
 149                         dz2=atom[i].z-atom[j].z;
 150                         if(dz2>sz/2.0)
 151                                 dz2-=sz;
 152                         else if(dz2<-sz/2.0)
 153                                 dz2+=sz;
 154                         dz2*=dz2;
 155                         dist=sqrt(dx2+dy2+dz2);
 156                         if(dist>=MAXR)
 157                                 continue;
 158                         k=dist/DELTA;
 159                         if((atom[i].type=='S')&&(atom[j].type=='S')) {
 160                                 aslot[k]+=1;
 161                         }
 162                         else if((atom[i].type=='C')&&(atom[j].type=='C')) {
 163                                 bslot[k]+=1;
 164                         }
 165                         else {
 166                                 cslot[k]+=1;
 167                         }
 168                 }
 169         }
 170
 171                 close(fd);
 172         }
 173
 174         // normalization and output
 175
 176         for(i=1;i<slots;i++) {
 177                 norm=4*M_PI*(i*DELTA)*(i*DELTA)*DELTA;
 178                 printf("pc: %f %f %f %f\n",i*DELTA,
 179                                            aslot[i]/norm,
 180                                            bslot[i]/norm,
 181                                            cslot[i]/norm);
 182         }
 183
 184         free(aslot);
 185         free(bslot);
 186         free(cslot);
 187
 188         free(atom);
 189
 190         return 0;
 191
 192 }
 193