added a pc calc for ascci files produced by outcar2moldyn ...

diff --git a/vasp_tools/pc_calc.c b/vasp_tools/pc_calc.c
new file mode 100644 (file)
index 0000000..10275d2
--- /dev/null
+++ b/vasp_tools/pc_calc.c
@@ -0,0 +1,190 @@
+/*
+ * pc_calc.c - calculate the pc
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+
+#include <math.h>
+
+
+#define MAXR   6.0
+#define DELTA  0.01
+
+typedef struct s_atom {
+       double x,y,z;
+       char type;
+} t_atom;
+
+int get_line(int fd,char *line,int max) {
+
+        int count,ret;
+
+        count=0;
+
+        while(1) {
+                if(count==max) return count;
+                ret=read(fd,line+count,1);
+                if(ret<=0) return ret;
+                if(line[count]=='\n') {
+                        memset(line+count,0,max-count-1);
+                        //line[count]='\0';
+                        return count+1;
+                }
+                count+=1;
+        }
+}
+
+int main(int argc,char **argv) {
+
+       t_atom *atom;
+       double *aslot,*bslot,*cslot;
+       int acnt,bcnt,ccnt,cnt,count,slots;
+       int fd;
+       char buf[256],*wptr;
+       int i,j,k;
+       double dx2,dy2,dz2,dist,norm;
+       double sx,sy,sz;
+
+       if(argc!=5) {
+               printf("usage: %s file sx sy sz\n",argv[0]);
+               return -1;
+       }
+
+       fd=open(argv[1],O_RDONLY);
+       if(fd<0) {
+               perror("open file\n");
+               return fd;
+       }
+
+       sx=atof(argv[2]);
+       sy=atof(argv[3]);
+       sz=atof(argv[4]);
+
+       // first line
+       cnt=get_line(fd,buf,256);
+
+       // read in atoms
+       count=0;
+       while(1) {
+               cnt=get_line(fd,buf,256);
+               if(cnt<=0)
+                       break;
+
+               // there is (another) atom      
+
+               atom=realloc(atom,(count+1)*sizeof(t_atom));
+               wptr=strtok(buf," ");
+               atom[count].type=wptr[0];
+
+               wptr=strtok(NULL," ");
+               atom[count].x=atof(wptr);
+               wptr=strtok(NULL," ");
+               atom[count].y=atof(wptr);
+               wptr=strtok(NULL," ");
+               atom[count].z=atof(wptr);
+
+               count+=1;
+       }
+
+       printf("i: read in %d atoms ...\n",count);
+
+       // prepare pc
+
+       slots=MAXR/DELTA;
+       aslot=malloc(slots*sizeof(double));
+       if(aslot==NULL) {
+               perror("slot a\n");
+               return -1;
+       }
+       memset(aslot,0,slots*sizeof(double));
+       bslot=malloc(slots*sizeof(double));
+       if(bslot==NULL) {
+               perror("slot a\n");
+               return -1;
+       }
+       memset(bslot,0,slots*sizeof(double));
+       cslot=malloc(slots*sizeof(double));
+       if(cslot==NULL) {
+               perror("slot a\n");
+               return -1;
+       }
+       memset(cslot,0,slots*sizeof(double));
+
+       printf("i: allocated 3 times %d slots ...\n",slots);
+
+       acnt=0;
+       bcnt=0;
+       ccnt=0;
+
+       // calc pc
+
+       for(i=0;i<count;i++) {
+               for(j=0;j<i;j++) {
+                       dx2=atom[i].x-atom[j].x;
+                       if(dx2>sx/2.0)
+                               dx2-=sx;
+                       else if(dx2<-sx/2.0)
+                               dx2+=sx;
+                       dx2*=dx2;
+                       dy2=atom[i].y-atom[j].y;
+                       if(dy2>sy/2.0)
+                               dy2-=sy;
+                       else if(dy2<-sy/2.0)
+                               dy2+=sy;
+                       dy2*=dy2;
+                       dz2=atom[i].z-atom[j].z;
+                       if(dz2>sz/2.0)
+                               dz2-=sz;
+                       else if(dz2<-sz/2.0)
+                               dz2+=sz;
+                       dz2*=dz2;
+                       dist=sqrt(dx2+dy2+dz2);
+                       if(dist>=MAXR)
+                               continue;
+                       k=dist/DELTA;
+                       if((atom[i].type=='S')&&(atom[j].type=='S')) {
+                               aslot[k]+=1;
+                               acnt+=1;
+                       }
+                       else if((atom[i].type=='C')&&(atom[j].type=='C')) {
+                               bslot[k]+=1;
+                               bcnt+=1;
+                       }
+                       else {
+                               cslot[k]+=1;
+                               ccnt+=1;
+                       }
+               }
+       }
+
+       // normalization and output
+
+       for(i=1;i<slots;i++) {
+               norm=4*M_PI*(i*DELTA)*(i*DELTA)*DELTA;
+               printf("pc: %f %f %f %f\n",i*DELTA,
+                                          aslot[i]/norm,
+                                          bslot[i]/norm,
+                                          cslot[i]/norm);
+       }
+
+       free(aslot);
+       free(bslot);
+       free(cslot);
+
+       free(atom);
+
+       close(fd);
+
+       return 0;
+
+}
+