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