X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=vasp_tools%2Fpc_calc.c;h=db12f50c6e1a8f70cb758a6f17199a602e92d049;hb=facb6c3b0ab8d6a6dc1c9e15298fb2ec5e0317f5;hp=10275d2d5d8b957765d79657b3b74e5cbb7d645b;hpb=0cc57bba82ec0e4208f48caf793726050a30b998;p=physik%2Fposic.git

diff --git a/vasp_tools/pc_calc.c b/vasp_tools/pc_calc.c
index 10275d2..db12f50 100644
--- a/vasp_tools/pc_calc.c
+++ b/vasp_tools/pc_calc.c
@@ -43,62 +43,76 @@ int get_line(int fd,char *line,int max) {
         }
 }
 
-int main(int argc,char **argv) {
+int get_max3_vals(char *buf,int amount,double *v1, double *v2,double *v3) {
 
-	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;
+	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;
 	}
 
-	fd=open(argv[1],O_RDONLY);
-	if(fd<0) {
-		perror("open file\n");
-		return fd;
-	}
+	return 0;
+}
 
-	sx=atof(argv[2]);
-	sy=atof(argv[3]);
-	sz=atof(argv[4]);
+int get_2_ints(char *buf,int *i1,int *i2) {
 
-	// first line
-	cnt=get_line(fd,buf,256);
+	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;
+	}
 
-	// read in atoms
-	count=0;
-	while(1) {
-		cnt=get_line(fd,buf,256);
-		if(cnt<=0)
-			break;
+	return 0;
+}
 
-		// there is (another) atom	
+int main(int argc,char **argv) {
 
-		atom=realloc(atom,(count+1)*sizeof(t_atom));
-		wptr=strtok(buf," ");
-		atom[count].type=wptr[0];
+	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;
 
-		wptr=strtok(NULL," ");
-		atom[count].x=atof(wptr);
-		wptr=strtok(NULL," ");
-		atom[count].y=atof(wptr);
-		wptr=strtok(NULL," ");
-		atom[count].z=atof(wptr);
+	atom=NULL;
 
-		count+=1;
+	if(argc<1) {
+		printf("usage: %s file1 file2 ...\n",
+		       argv[0]);
+		return -1;
 	}
 
-	printf("i: read in %d atoms ...\n",count);
-
 	// prepare pc
-
 	slots=MAXR/DELTA;
 	aslot=malloc(slots*sizeof(double));
 	if(aslot==NULL) {
@@ -118,58 +132,132 @@ int main(int argc,char **argv) {
 		return -1;
 	}
 	memset(cslot,0,slots*sizeof(double));
-
 	printf("i: allocated 3 times %d slots ...\n",slots);
 
-	acnt=0;
-	bcnt=0;
-	ccnt=0;
+	// use all given files ...
+	printf("using files:\n");
+	for(fcnt=1;fcnt<argc;fcnt++)
+		printf(" %d: %s\n",fcnt,argv[fcnt]);
 
-	// calc pc
+	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++) {
-			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);
+			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;
-				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
+	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,
@@ -182,8 +270,6 @@ int main(int argc,char **argv) {
 
 	free(atom);
 
-	close(fd);
-
 	return 0;
 
 }