From 0cc57bba82ec0e4208f48caf793726050a30b998 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Sun, 13 Sep 2009 13:20:49 +0200
Subject: [PATCH] added a pc calc for ascci files produced by outcar2moldyn ...

---
 vasp_tools/pc_calc.c | 190 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 190 insertions(+)
 create mode 100644 vasp_tools/pc_calc.c

diff --git a/vasp_tools/pc_calc.c b/vasp_tools/pc_calc.c
new file mode 100644
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;
+
+}
+
-- 
2.39.2