X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=1aff5a8df207c2b5d22077c70ccbc31c128afdd6;hb=2b6a06cdf69b0d8f1dafc0e21b89c6bb4bfc192c;hp=a1bfcaea95c82580dc636f5655649a61d0e8cf8f;hpb=155e1cfea83209d09c2a06ae4fb7f5e1652fc00a;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index a1bfcae..1aff5a8 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -30,55 +30,6 @@
 #include "potentials/tersoff.h"
 #endif
 
-
-/*
- * global variables, pse and atom colors (only needed here)
- */ 
-
-static char *pse_name[]={
-	"*",
-	"H",
-	"He",
-	"Li",
-	"Be",
-	"B",
-	"C",
-	"N",
-	"O",
-	"F",
-	"Ne",
-	"Na",
-	"Mg",
-	"Al",
-	"Si",
-	"P",
-	"S",
-	"Cl",
-	"Ar",
-};
-
-static char *pse_col[]={
-	"*",
-	"White",
-	"He",
-	"Li",
-	"Be",
-	"B",
-	"Gray",
-	"N",
-	"Blue",
-	"F",
-	"Ne",
-	"Na",
-	"Mg",
-	"Al",
-	"Yellow",
-	"P",
-	"S",
-	"Cl",
-	"Ar",
-};
-
 /*
  * the moldyn functions
  */
@@ -131,24 +82,13 @@ int set_int_alg(t_moldyn *moldyn,u8 algo) {
 int set_cutoff(t_moldyn *moldyn,double cutoff) {
 
 	moldyn->cutoff=cutoff;
+	moldyn->cutoff_square=cutoff*cutoff;
 
 	printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff);
 
 	return 0;
 }
 
-int set_bondlen(t_moldyn *moldyn,double b0,double b1,double bm) {
-
-	moldyn->bondlen[0]=b0*b0;
-	moldyn->bondlen[1]=b1*b1;
-	if(bm<0)
-		moldyn->bondlen[2]=b0*b1;
-	else
-		moldyn->bondlen[2]=bm*bm;
-
-	return 0;
-}
-
 int set_temperature(t_moldyn *moldyn,double t_ref) {
 
 	moldyn->t_ref=t_ref;
@@ -167,6 +107,38 @@ int set_pressure(t_moldyn *moldyn,double p_ref) {
 	return 0;
 }
 
+int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) {
+
+	moldyn->pt_scale&=(~(P_SCALE_MASK));
+	moldyn->pt_scale|=ptype;
+	moldyn->p_tc=ptc;
+
+	printf("[moldyn] p scaling:\n");
+
+	printf("  p: %s",ptype?"yes":"no ");
+	if(ptype)
+		printf(" | type: %02x | factor: %f",ptype,ptc);
+	printf("\n");
+
+	return 0;
+}
+
+int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) {
+
+	moldyn->pt_scale&=(~(T_SCALE_MASK));
+	moldyn->pt_scale|=ttype;
+	moldyn->t_tc=ttc;
+
+	printf("[moldyn] t scaling:\n");
+
+	printf("  t: %s",ttype?"yes":"no ");
+	if(ttype)
+		printf(" | type: %02x | factor: %f",ttype,ttc);
+	printf("\n");
+
+	return 0;
+}
+
 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
 
 	moldyn->pt_scale=(ptype|ttype);
@@ -243,13 +215,14 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 
 	switch(type) {
 		case MOLDYN_POTENTIAL_TM:
-			moldyn->func1b=tersoff_mult_1bp;
-			moldyn->func3b_j1=tersoff_mult_3bp_j1;
-			moldyn->func3b_k1=tersoff_mult_3bp_k1;
-			moldyn->func3b_j2=tersoff_mult_3bp_j2;
-			moldyn->func3b_k2=tersoff_mult_3bp_k2;
-			// missing: check 2b bond func
+			moldyn->func_i0=tersoff_mult_1bp;
+			moldyn->func_j1=tersoff_mult_3bp_j1;
+			moldyn->func_j1_k0=tersoff_mult_3bp_k1;
+			moldyn->func_j1c=tersoff_mult_3bp_j2;
+			moldyn->func_j1_k1=tersoff_mult_3bp_k2;
+			moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
 			break;
+		/*
 		case MOLDYN_POTENTIAL_AM:
 			moldyn->func3b_j1=albe_mult_3bp_j1;
 			moldyn->func3b_k1=albe_mult_3bp_k1;
@@ -257,12 +230,23 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 			moldyn->func3b_k2=albe_mult_3bp_k2;
 			moldyn->check_2b_bond=albe_mult_check_2b_bond;
 			break;
+		*/
+		case MOLDYN_POTENTIAL_AM:
+			moldyn->func_i0=albe_mult_i0;
+			moldyn->func_j0=albe_mult_i0_j0;
+			moldyn->func_j0_k0=albe_mult_i0_j0_k0;
+			moldyn->func_j0e=albe_mult_i0_j1;
+			moldyn->func_j1=albe_mult_i0_j2;
+			moldyn->func_j1_k0=albe_mult_i0_j2_k0;
+			moldyn->func_j1c=albe_mult_i0_j3;
+			moldyn->check_2b_bond=albe_mult_check_2b_bond;
+			break;
 		case MOLDYN_POTENTIAL_HO:
-			moldyn->func2b=harmonic_oscillator;
+			moldyn->func_j0=harmonic_oscillator;
 			moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
 			break;
 		case MOLDYN_POTENTIAL_LJ:
-			moldyn->func2b=lennard_jones;
+			moldyn->func_j0=lennard_jones;
 			moldyn->check_2b_bond=lennard_jones_check_2b_bond;
 			break;
 		default:
@@ -507,6 +491,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	t_3dvec orig;
 	void *ptr;
 	t_atom *atom;
+	char name[16];
 
 	new=a*b*c;
 	count=moldyn->count;
@@ -541,18 +526,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 		case CUBIC:
 			set_nn_dist(moldyn,lc);
 			ret=cubic_init(a,b,c,lc,atom,&orig);
+			strcpy(name,"cubic");
 			break;
 		case FCC:
 			if(!origin)
 				v3_scale(&orig,&orig,0.5);
 			set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
 			ret=fcc_init(a,b,c,lc,atom,&orig);
+			strcpy(name,"fcc");
 			break;
 		case DIAMOND:
 			if(!origin)
 				v3_scale(&orig,&orig,0.25);
 			set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
 			ret=diamond_init(a,b,c,lc,atom,&orig);
+			strcpy(name,"diamond");
 			break;
 		default:
 			printf("unknown lattice type (%02x)\n",type);
@@ -569,7 +557,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	}
 
 	moldyn->count+=new;
-	printf("[moldyn] created lattice with %d atoms\n",new);
+	printf("[moldyn] created %s lattice with %d atoms\n",name,new);
 
 	for(ret=0;ret<new;ret++) {
 		atom[ret].element=element;
@@ -1581,7 +1569,6 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
 	/* sqaure of some variables */
 	moldyn->tau_square=moldyn->tau*moldyn->tau;
-	moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
 
 	/* get current time */
 	gettimeofday(&t1,NULL);
@@ -1604,8 +1591,9 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		printf("[moldyn] WARNING: forces too high / tau too small!\n");
 
 	/* zero absolute time */
-	moldyn->time=0.0;
-	moldyn->total_steps=0;
+	// should have right values!
+	//moldyn->time=0.0;
+	//moldyn->total_steps=0;
 
 	/* debugging, ignore */
 	moldyn->debug=0;
@@ -1637,7 +1625,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		temperature_calc(moldyn);
 		virial_sum(moldyn);
 		pressure_calc(moldyn);
-		//thermodynamic_pressure_calc(moldyn);
+		/*
+		thermodynamic_pressure_calc(moldyn);
+		printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
+		       moldyn->tp/BAR);
+		*/
 
 		/* calculate fluctuations + averages */
 		average_and_fluctuation_calc(moldyn);
@@ -1868,8 +1860,8 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 		/* single particle potential/force */
 		if(itom[i].attr&ATOM_ATTR_1BP)
-			if(moldyn->func1b)
-				moldyn->func1b(moldyn,&(itom[i]));
+			if(moldyn->func_i0)
+				moldyn->func_i0(moldyn,&(itom[i]));
 
 		if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
 			continue;
@@ -1884,8 +1876,14 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 		dnlc=lc->dnlc;
 
+#ifndef STATIC_LISTS
+		/* check for later 3 body interaction */
+		if(itom[i].attr&ATOM_ATTR_3BP)
+			memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
+#endif
+
 		/* first loop over atoms j */
-		if(moldyn->func2b) {
+		if(moldyn->func_j0) {
 			for(j=0;j<27;j++) {
 
 				bc_ij=(j<dnlc)?0:1;
@@ -1896,18 +1894,6 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 					jtom=&(atom[neighbour_i[j][p]]);
 					p++;
-
-					if(jtom==&(itom[i]))
-						continue;
-
-					if((jtom->attr&ATOM_ATTR_2BP)&
-					   (itom[i].attr&ATOM_ATTR_2BP)) {
-						moldyn->func2b(moldyn,
-						               &(itom[i]),
-						               jtom,
-					        	       bc_ij);
-					}
-				}
 #else
 				this=&(neighbour_i[j]);
 				list_reset_f(this);
@@ -1917,35 +1903,103 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 				do {
 					jtom=this->current->data;
+#endif
 
 					if(jtom==&(itom[i]))
 						continue;
 
+					/* reset 3bp run */
+					moldyn->run3bp=1;
+
 					if((jtom->attr&ATOM_ATTR_2BP)&
 					   (itom[i].attr&ATOM_ATTR_2BP)) {
-						moldyn->func2b(moldyn,
-						               &(itom[i]),
-						               jtom,
-					        	       bc_ij);
+						moldyn->func_j0(moldyn,
+						                &(itom[i]),
+						                jtom,
+					        	        bc_ij);
+					}
+
+					/* 3 body potential/force */
+
+					/* in j loop, 3bp run can be skipped */
+					if(!(moldyn->run3bp))
+						continue;
+
+					if(!(itom[i].attr&ATOM_ATTR_3BP))
+						continue;
+
+					if(!(jtom->attr&ATOM_ATTR_3BP))
+						continue;
+
+					if(moldyn->func_j0_k0==NULL)
+						continue;
+
+				/* first loop over atoms k in first j loop */
+				for(k=0;k<27;k++) {
+
+					bc_ik=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+					q=0;
+
+					while(neighbour_i[j][q]!=0) {
+
+						ktom=&(atom[neighbour_i[k][q]]);
+						q++;
+#else
+					that=&(neighbour_i2[k]);
+					list_reset_f(that);
+					
+					if(that->start==NULL)
+						continue;
+
+					do {
+						ktom=that->current->data;
+#endif
+
+						if(!(ktom->attr&ATOM_ATTR_3BP))
+							continue;
+
+						//if(ktom==jtom)
+						//	continue;
+
+						if(ktom==&(itom[i]))
+							continue;
+
+						moldyn->func_j0_k0(moldyn,
+						                   &(itom[i]),
+						                   jtom,
+						                   ktom,
+						                   bc_ik|bc_ij);
+#ifdef STATIC_LISTS
 					}
+#else
+					} while(list_next_f(that)!=\
+					        L_NO_NEXT_ELEMENT);
+#endif
+
+				}
+
+				/* finish of first j loop after first k loop */
+				if(moldyn->func_j0e)
+					moldyn->func_j0e(moldyn,
+					                 &(itom[i]),
+					                 jtom,
+					                 bc_ij);
+
+#ifdef STATIC_LISTS
+				}
+#else
 				} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
 
 			}
 		}
 
-		/* 3 body potential/force */
+		/* continued 3 body potential/force */
 
 		if(!(itom[i].attr&ATOM_ATTR_3BP))
 			continue;
 
-		/* copy the neighbour lists */
-#ifdef STATIC_LISTS
-		/* no copy needed for static lists */
-#else
-		memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
-#endif
-
 		/* second loop over atoms j */
 		for(j=0;j<27;j++) {
 
@@ -1978,18 +2032,18 @@ int potential_force_calc(t_moldyn *moldyn) {
 				/* reset 3bp run */
 				moldyn->run3bp=1;
 
-				if(moldyn->func3b_j1)
-					moldyn->func3b_j1(moldyn,
-					                  &(itom[i]),
-					                  jtom,
-					                  bc_ij);
+				if(moldyn->func_j1)
+					moldyn->func_j1(moldyn,
+					                &(itom[i]),
+					                jtom,
+					                bc_ij);
 
-				/* in first j loop, 3bp run can be skipped */
+				/* in j loop, 3bp run can be skipped */
 				if(!(moldyn->run3bp))
 					continue;
 			
-				/* first loop over atoms k */
-				if(moldyn->func3b_k1) {
+				/* first loop over atoms k in second j loop */
+				if(moldyn->func_j1_k0) {
 
 				for(k=0;k<27;k++) {
 
@@ -2015,17 +2069,17 @@ int potential_force_calc(t_moldyn *moldyn) {
 						if(!(ktom->attr&ATOM_ATTR_3BP))
 							continue;
 
-						if(ktom==jtom)
-							continue;
+						//if(ktom==jtom)
+						//	continue;
 
 						if(ktom==&(itom[i]))
 							continue;
 
-						moldyn->func3b_k1(moldyn,
-						                  &(itom[i]),
-						                  jtom,
-						                  ktom,
-						                  bc_ik|bc_ij);
+						moldyn->func_j1_k0(moldyn,
+						                   &(itom[i]),
+						                   jtom,
+						                   ktom,
+						                   bc_ik|bc_ij);
 #ifdef STATIC_LISTS
 					}
 #else
@@ -2037,14 +2091,15 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 				}
 
-				if(moldyn->func3b_j2)
-					moldyn->func3b_j2(moldyn,
-					                  &(itom[i]),
-					                  jtom,
-					                  bc_ij);
+				/* continued j loop after first k loop */
+				if(moldyn->func_j1c)
+					moldyn->func_j1c(moldyn,
+					                 &(itom[i]),
+					                 jtom,
+					                 bc_ij);
 
 				/* second loop over atoms k */
-				if(moldyn->func3b_k2) {
+				if(moldyn->func_j1_k1) {
 
 				for(k=0;k<27;k++) {
 
@@ -2070,17 +2125,17 @@ int potential_force_calc(t_moldyn *moldyn) {
 						if(!(ktom->attr&ATOM_ATTR_3BP))
 							continue;
 
-						if(ktom==jtom)
-							continue;
+						//if(ktom==jtom)
+						//	continue;
 
 						if(ktom==&(itom[i]))
 							continue;
 
-						moldyn->func3b_k2(moldyn,
-						                  &(itom[i]),
-						                  jtom,
-						                  ktom,
-						                  bc_ik|bc_ij);
+						moldyn->func_j1_k1(moldyn,
+						                   &(itom[i]),
+						                   jtom,
+						                   ktom,
+						                   bc_ik|bc_ij);
 
 #ifdef STATIC_LISTS
 					}
@@ -2093,11 +2148,11 @@ int potential_force_calc(t_moldyn *moldyn) {
 				
 				}
 
-				/* 2bp post function */
-				if(moldyn->func3b_j3) {
-					moldyn->func3b_j3(moldyn,
-					                  &(itom[i]),
-					                  jtom,bc_ij);
+				/* finish of j loop after second k loop */
+				if(moldyn->func_j1e) {
+					moldyn->func_j1e(moldyn,
+					                 &(itom[i]),
+					                 jtom,bc_ij);
 				}
 #ifdef STATIC_LISTS
 			}
@@ -2106,7 +2161,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 #endif
 		
 		}
-		
+
 #ifdef DEBUG
 	//printf("\n\n");
 #endif
@@ -2348,9 +2403,6 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
 	t_list *this;
 
 	lc=&(moldyn->lc);
-
-	link_cell_init(moldyn,VERBOSE);
-
 	itom=moldyn->atom;
 	
 	for(i=0;i<moldyn->count;i++) {
@@ -2414,7 +2466,8 @@ int get_line(int fd,char *line,int max) {
 		ret=read(fd,line+count,1);
 		if(ret<=0) return ret;
 		if(line[count]=='\n') {
-			line[count]='\0';
+			memset(line+count,0,max-count-1);
+			//line[count]='\0';
 			return count+1;
 		}
 		count+=1;
@@ -2609,8 +2662,6 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 	if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
 		return 0;
 
-	d=sqrt(d);
-
 	/* now count this bonding ... */
 	ba=data;
 
@@ -2705,39 +2756,48 @@ int visual_init(t_moldyn *moldyn,char *filebase) {
 	return 0;
 }
 
+int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+                         void *data,u8 bc) {
+
+	t_vb *vb;
+
+	vb=data;
+
+	if(itom->tag>=jtom->tag)
+		return 0;
+	
+	if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
+		return 0;
+
+	if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB))
+		dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n",
+		        itom->r.x,itom->r.y,itom->r.z,
+		        jtom->r.x,jtom->r.y,jtom->r.z);
+
+	return 0;
+}
+
 int visual_atoms(t_moldyn *moldyn) {
 
-	int i,j,fd;
+	int i;
  	char file[128+64];
 	t_3dvec dim;
 	double help;
 	t_visual *v;
 	t_atom *atom;
-	t_atom *btom;
-	t_linkcell *lc;
-#ifdef STATIC_LISTS
-	int *neighbour[27];
-	int p;
-#else
-	t_list neighbour[27];
-#endif
-	u8 bc;
-	t_3dvec dist;
-	double d2;
-	u8 brand;
+	t_vb vb;
 
 	v=&(moldyn->vis);
 	dim.x=v->dim.x;
 	dim.y=v->dim.y;
 	dim.z=v->dim.z;
 	atom=moldyn->atom;
-	lc=&(moldyn->lc);
 
 	help=(dim.x+dim.y);
 
 	sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
-	fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
-	if(fd<0) {
+	vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+	if(vb.fd<0) {
 		perror("open visual save file fd");
 		return -1;
 	}
@@ -2745,118 +2805,65 @@ int visual_atoms(t_moldyn *moldyn) {
 	/* write the actual data file */
 
 	// povray header
-	dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n",
+	dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n",
 	        moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
 
 	// atomic configuration
-	for(i=0;i<moldyn->count;i++) {
+	for(i=0;i<moldyn->count;i++)
 		// atom type, positions, color and kinetic energy
-		dprintf(fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
-		                                 atom[i].r.x,
-		                                 atom[i].r.y,
-		                                 atom[i].r.z,
-		                                 pse_col[atom[i].element],
-		                                 atom[i].ekin);
-
-		/*
-		 * bond detection should usually be done by potential
-		 * functions. brrrrr! EVIL!
-		 * 
-		 * todo: potentials need to export a 'find_bonds' function!
-		 */
-
-		// bonds between atoms
-		if(!(atom[i].attr&ATOM_ATTR_VB))
-			continue;
-		link_cell_neighbour_index(moldyn,
-		                          (atom[i].r.x+moldyn->dim.x/2)/lc->x,
-		                          (atom[i].r.y+moldyn->dim.y/2)/lc->y,
-		                          (atom[i].r.z+moldyn->dim.z/2)/lc->z,
-		                          neighbour);
-		for(j=0;j<27;j++) {
-			bc=j<lc->dnlc?0:1;
-#ifdef STATIC_LISTS
-			p=0;
-			while(neighbour[j][p]!=0) {
-				btom=&(atom[neighbour[j][p]]);
-				p++;
-#else
-			list_reset_f(&neighbour[j]);
-			if(neighbour[j].start==NULL)
-				continue;
-			do {
-				btom=neighbour[j].current->data;
-#endif
-				if(btom==&atom[i])	// skip identical atoms
-					continue;
-				//if(btom<&atom[i])	// skip half of them
-				//	continue;
-				v3_sub(&dist,&(atom[i].r),&(btom->r));
-				if(bc) check_per_bound(moldyn,&dist);
-				d2=v3_absolute_square(&dist);
-				brand=atom[i].brand;
-				if(brand==btom->brand) {
-					if(d2>moldyn->bondlen[brand])
-						continue;
-				}
-				else {
-					if(d2>moldyn->bondlen[2])
-						continue;
-				}
-				dprintf(fd,"# [B] %f %f %f %f %f %f\n",
-				        atom[i].r.x,atom[i].r.y,atom[i].r.z,
-		       		        btom->r.x,btom->r.y,btom->r.z);
-#ifdef STATIC_LISTS
-			}
-#else
-			} while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT);
-#endif
-		}
-	}
-
+		dprintf(vb.fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
+		                                    atom[i].r.x,
+		                                    atom[i].r.y,
+		                                    atom[i].r.z,
+		                                    pse_col[atom[i].element],
+		                                    atom[i].ekin);
+	
+	// bonds between atoms
+	process_2b_bonds(moldyn,&vb,visual_bonds_process);
+	
 	// boundaries
 	if(dim.x) {
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,-dim.y/2,-dim.z/2,
 		        dim.x/2,-dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,-dim.y/2,-dim.z/2,
 		        -dim.x/2,dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        dim.x/2,dim.y/2,-dim.z/2,
 		        dim.x/2,-dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,dim.y/2,-dim.z/2,
 		        dim.x/2,dim.y/2,-dim.z/2);
 
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,-dim.y/2,dim.z/2,
 		        dim.x/2,-dim.y/2,dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,-dim.y/2,dim.z/2,
 		        -dim.x/2,dim.y/2,dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        dim.x/2,dim.y/2,dim.z/2,
 		        dim.x/2,-dim.y/2,dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,dim.y/2,dim.z/2,
 		        dim.x/2,dim.y/2,dim.z/2);
 
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,-dim.y/2,dim.z/2,
 		        -dim.x/2,-dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,dim.y/2,dim.z/2,
 		        -dim.x/2,dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        dim.x/2,-dim.y/2,dim.z/2,
 		        dim.x/2,-dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        dim.x/2,dim.y/2,dim.z/2,
 		        dim.x/2,dim.y/2,-dim.z/2);
 	}
 
-	close(fd);
+	close(vb.fd);
 
 	return 0;
 }