X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=vasp_tools%2Fpc_calc.c;fp=vasp_tools%2Fpc_calc.c;h=db12f50c6e1a8f70cb758a6f17199a602e92d049;hp=0000000000000000000000000000000000000000;hb=6e6d7126ea9a845f11637d8e1b8eb2b570ac4dc9;hpb=97dc63eb6a519b8e1f4fbfaa9760dd94539436b0

diff --git a/vasp_tools/pc_calc.c b/vasp_tools/pc_calc.c
new file mode 100644
index 0000000..db12f50
--- /dev/null
+++ b/vasp_tools/pc_calc.c
@@ -0,0 +1,276 @@
+/*
+ * 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 get_max3_vals(char *buf,int amount,double *v1, double *v2,double *v3) {
+
+	int i,j,k;
+	char temp[128];
+
+	i=0;
+
+	for(k=0;k<amount;k++) {
+		while(buf[i]==' ')
+			i++;
+		j=0;
+		while((buf[i+j]!=' ')&&(buf[i+j]!=0))
+			j++;
+		memcpy(temp,buf+i,j);
+		temp[j]=0;
+		if(k==0) *v1=atof(temp);
+		if(k==1) *v2=atof(temp);
+		if(k==2) *v3=atof(temp);
+		i+=j;
+	}
+
+	return 0;
+}
+
+int get_2_ints(char *buf,int *i1,int *i2) {
+
+	int i,j,k;
+	char temp[128];
+
+	i=0;
+
+	for(k=0;k<2;k++) {
+		while(buf[i]==' ')
+			i++;
+		j=0;
+		while((buf[i+j]!=' ')&&(buf[i+j]!=0))
+			j++;
+		memcpy(temp,buf+i,j);
+		temp[j]=0;
+		if(k==0) *i1=atoi(temp);
+		if(k==1) *i2=atoi(temp);
+		i+=j;
+	}
+
+	return 0;
+}
+
+int main(int argc,char **argv) {
+
+	t_atom *atom;
+	double *aslot,*bslot,*cslot;
+	int cnt,count,fcnt,slots;
+	int fd;
+	char buf[256];
+	int i,j,k;
+	int nsi,nc,ntot;
+	double dx,dy,dz,dist,norm;
+	double dxt,dyt,dzt;
+	double X1,X2,X3,Y1,Y2,Y3,Z1,Z2,Z3;
+	double scale;
+
+	atom=NULL;
+
+	if(argc<1) {
+		printf("usage: %s file1 file2 ...\n",
+		       argv[0]);
+		return -1;
+	}
+
+	// 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);
+
+	// use all given files ...
+	printf("using files:\n");
+	for(fcnt=1;fcnt<argc;fcnt++)
+		printf(" %d: %s\n",fcnt,argv[fcnt]);
+
+	for(fcnt=1;fcnt<argc;fcnt++) {
+
+	fd=open(argv[fcnt],O_RDONLY);
+	if(fd<0) {
+		perror("open file\n");
+		return fd;
+	}
+
+	// first line
+	cnt=get_line(fd,buf,256);
+
+	// second line -> scale
+	cnt=get_line(fd,buf,256);
+	get_max3_vals(buf,1,&scale,NULL,NULL);
+	printf("scale: %f\n",scale);
+
+	// third line -> X
+	cnt=get_line(fd,buf,256);
+	get_max3_vals(buf,3,&X1,&X2,&X3);
+	printf("X: %f %f %f\n",X1,X2,X3);
+	
+	// 4th line -> Y
+	cnt=get_line(fd,buf,256);
+	get_max3_vals(buf,3,&Y1,&Y2,&Y3);
+	printf("Y: %f %f %f\n",Y1,Y2,Y3);
+	
+	// 5th line -> Z
+	cnt=get_line(fd,buf,256);
+	get_max3_vals(buf,3,&Z1,&Z2,&Z3);
+	printf("Z: %f %f %f\n",Z1,Z2,Z3);
+	
+	// 6th line -> nsi + nc = ntot
+	cnt=get_line(fd,buf,256);
+	get_2_ints(buf,&nsi,&nc);
+	ntot=nsi+nc;
+	printf("Si: %d - C: %d - tot: %d\n",nsi,nc,ntot);
+	
+	// 7th line
+	cnt=get_line(fd,buf,256);
+
+	// 8th line
+	cnt=get_line(fd,buf,256);
+
+	// read in atoms
+	count=0;
+	while(count<ntot) {
+		cnt=get_line(fd,buf,256);
+		if(cnt<=0) {
+			printf("Should never happen!\n");
+			return cnt;
+		}
+
+		atom=realloc(atom,(count+1)*sizeof(t_atom));
+		if(count<nsi)
+			atom[count].type='S';
+		else
+			atom[count].type='C';
+
+		get_max3_vals(buf,3,
+		              &(atom[count].x),
+		              &(atom[count].y),
+		              &(atom[count].z));
+
+		count+=1;
+	}
+
+	printf("i: read in %d atoms ...\n",count);
+
+	// calc pc
+	for(i=0;i<count;i++) {
+		for(j=0;j<i;j++) {
+			dx=atom[i].x-atom[j].x;
+			if(dx>0.5)
+				dx-=1;
+			else if(dx<-0.5)
+				dx+=1;
+
+			dy=atom[i].y-atom[j].y;
+			if(dy>0.5)
+				dy-=1;
+			else if(dy<-0.5)
+				dy+=1;
+
+			dz=atom[i].z-atom[j].z;
+			if(dz>0.5)
+				dz-=1;
+			else if(dz<-0.5)
+				dz+=1;
+
+			dxt=dx*X1+dy*Y1+dz*Z1;
+			dyt=dx*X2+dy*Y2+dz*Z2;
+			dzt=dx*X3+dy*Y3+dz*Z3;
+
+			dist=sqrt(dxt*dxt+dyt*dyt+dzt*dzt);
+			dist*=scale;
+
+			if(dist>=MAXR)
+				continue;
+			k=dist/DELTA;
+			if((atom[i].type=='S')&&(atom[j].type=='S')) {
+				aslot[k]+=1;
+			}
+			else if((atom[i].type=='C')&&(atom[j].type=='C')) {
+				bslot[k]+=1;
+			}
+			else {
+				cslot[k]+=1;
+			}
+		}
+	}
+
+	close(fd);
+
+	}
+
+	// normalization and output
+	for(i=1;i<slots;i++) {
+		norm=4*M_PI*(i*DELTA)*(i*DELTA)*DELTA;
+		norm*=(argc-1);
+		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);
+
+	return 0;
+
+}
+