From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Wed, 17 Dec 2008 16:18:12 +0000 (+0100)
Subject: cleaned partial lattice create + added defect feature, adabted md code +
X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dfb47bde6bf974d1b231048ec8bcda712404ccc1;p=physik%2Fposic.git

cleaned partial lattice create + added defect feature, adabted md code +
some minor changes: ask diff! :) - xmas checkin ...
---

diff --git a/TODO.txt b/TODO.txt
index 570ee54..8712c9d 100644
--- a/TODO.txt
+++ b/TODO.txt
@@ -1,19 +1,23 @@
 interstitial simulations
 ************************
 
-  DOUBLECHECK with PARCAS !!!!
-
-  - visualize md.movie (converter tool - maybe good!)
-  - more checking ...
+-  more DOUBLECHECKs with PARCAS !!!!
 
 implement
 *********
 
  - improve ins_m_atoms (merge to general ins_atoms function maybe)
+
+   * modified first: try rules to better distribute Si and C	***
+   * general cleanup
+   * type 1 square, type 2 radius!
+
  - clean up mdrun
  - config sanity checks
  - introduce/improve fill command (multiple fills)
 
+ - anneal -> from current to T with rate R
+
  - check virial calc, where does the - come from?
 
  - angular distribution
@@ -25,8 +29,18 @@ implement
 
  - improve diff calc
 
- - make it parallel! (mpi/openmp) <- email to ralfu (asap)
-   * openmp: doch auf verlet listen (pointer problem!)
+ - pthreads:
+
+   * threads nur einmal oeffnen
+   * verteilung auf laufende threads
+
+ - potentials:
+   
+   * tersoff fast			***
+   * modifief tersoff/albe		***
+   * stillinger weber			***
+   * eam
+   * edip
 
  - optimize code!
 
@@ -44,6 +58,9 @@ implement
    * -> 3.50 ^ + static lists
    * -> 3.37 source c, albe f, arch opts, c,d,h,gamma
              + lowmem lists (test for bigger + atoms!)
+   * -> 3.31(36) ^ + inline v_calc
+
+   * -> 4.44 orig albe, lowmem, source c, arch opts
 
    todo:
     - listen ! estimate time
@@ -61,7 +78,7 @@ simulation runs
 - tctrl only in outer regions
 - only 1 atom per timestep
 - EXTENDED C-C cutoff
-- different sized SiC prec in Si (4:5! diff temperatures)
+- different sized SiC prec in Si (4:5! diff temperatures)	***
 - melting exps (both, anneal + interface method)
 - interstitials:
   - more interstitial combinations
diff --git a/calc_delta_e b/calc_delta_e
index 4ff46fb..455ce57 100755
--- a/calc_delta_e
+++ b/calc_delta_e
@@ -2,7 +2,7 @@
 
 file=`ls $1/atomic_conf_* | tail -1`
 atom_cnt=`grep '# \[P\]' $file | awk '{ print $3 }'`
-e0=`grep ^0 $1/energy | awk '{ print $4 }'`
+e0=`awk 'NR==2' $1/energy | awk '{ print $4 }'`
 e1=`tail -n 1 $1/energy | awk '{ print $4 }'`
 
 echo "$e0 $e1 $atom_cnt" | \
diff --git a/mdrun.c b/mdrun.c
index 6effe1a..d063693 100644
--- a/mdrun.c
+++ b/mdrun.c
@@ -277,66 +277,114 @@ int mdrun_parse_config(t_mdrun *mdrun) {
 				mdrun->lattice=DIAMOND;
 			if(!strncmp(word[1],"none",4))
 				mdrun->lattice=NONE;
+			if(wcnt==3)
+				mdrun->lc=atof(word[2]);
 		}
 		else if(!strncmp(word[0],"element1",8)) {
 			mdrun->element1=atoi(word[1]);
-			mdrun->m1=pse_mass[mdrun->element1];
 		}
 		else if(!strncmp(word[0],"element2",8)) {
 			mdrun->element2=atoi(word[1]);
-			mdrun->m2=pse_mass[mdrun->element2];
 		}
 		else if(!strncmp(word[0],"fill",6)) {
-			fp.lx=atoi(word[2]);
-			fp.ly=atoi(word[3]);
-			fp.lz=atoi(word[4]);
-			fp.lc=atof(word[5]);
-			mdrun->lc=fp.lc;
-			if(!strncmp(word[1],"lc",2)) {
-				if(wcnt==8) {
-					fp.fill_element=atoi(word[6]);
-					fp.fill_brand=atoi(word[7]);
+			// default values
+			fp.fill_element=mdrun->element1;
+			fp.fill_brand=0;
+			fp.lattice=mdrun->lattice;
+			fp.p_params.type=0;
+			fp.d_params.type=0;
+			// parse fill command
+			i=1;
+			while(i<wcnt) {
+				if(!strncmp(word[i],"lc",2)) {
+					fp.lx=atoi(word[++i]);
+					fp.ly=atoi(word[++i]);
+					fp.lz=atoi(word[++i]);
+					fp.lc=atof(word[++i]);
+					mdrun->lc=fp.lc;
 				}
-				else {
-					fp.fill_element=mdrun->element1;
-					fp.fill_brand=0;
+				if(!strncmp(word[i],"eb",2)) {
+					fp.fill_element=atoi(word[++i]);
+					fp.fill_brand=atoi(word[++i]);
 				}
-			}
-			else {
-				switch(word[6][0]) {
+				if(word[i][0]=='p') {
+					i+=1;
+					switch(word[i][0]) {
 					case 'i':
-						if(word[6][1]=='r')
-							fp.p_type=PART_INSIDE_R;
+						if(word[i][1]=='r')
+						fp.p_params.type=PART_INSIDE_R;
 						else
-							fp.p_type=PART_INSIDE_D;
+						fp.p_params.type=PART_INSIDE_D;
 						break;
 					case 'o':
-						if(word[6][1]=='r')
-						fp.p_type=PART_OUTSIDE_R;
+						if(word[i][1]=='r')
+						fp.p_params.type=PART_OUTSIDE_R;
 						else
-						fp.p_type=PART_OUTSIDE_D;
+						fp.p_params.type=PART_OUTSIDE_D;
 						break;
 					default:
 						break;
+					}
+					if((fp.p_params.type==PART_INSIDE_R)||
+                         		  (fp.p_params.type==PART_OUTSIDE_R)) {
+						fp.p_params.r=atof(word[++i]);	
+						fp.p_params.p.x=atof(word[++i]);
+						fp.p_params.p.y=atof(word[++i]);
+						fp.p_params.p.z=atof(word[++i]);
+					}
+					if((fp.p_params.type==PART_INSIDE_D)||
+					   (fp.p_params.type==PART_OUTSIDE_D)) {
+						fp.p_params.p.x=atof(word[++i]);
+						fp.p_params.p.y=atof(word[++i]);
+						fp.p_params.p.z=atof(word[++i]);
+						fp.p_params.d.x=atof(word[++i]);
+						fp.p_params.d.y=atof(word[++i]);
+						fp.p_params.d.z=atof(word[++i]);
+					}
 				}
+				if(word[i][0]=='d') {
+					switch(word[++i][0]) {
+						case '0':
+
+				fp.d_params.type=DEFECT_TYPE_0D;
+				if(!strncmp(word[i+1],"dbx",3)) {
+					fp.d_params.stype=DEFECT_STYPE_DB_X;
+				}
+				if(!strncmp(word[i+1],"dby",3)) {
+					fp.d_params.stype=DEFECT_STYPE_DB_Y;
+				}
+				if(!strncmp(word[i+1],"dbz",3)) {
+					fp.d_params.stype=DEFECT_STYPE_DB_Z;
+				}
+				if(!strncmp(word[i+1],"dbr",3)) {
+					fp.d_params.stype=DEFECT_STYPE_DB_R;
+				}
+				i+=1;
+				fp.d_params.od=atof(word[++i]);
+				fp.d_params.dd=atof(word[++i]);
+				fp.d_params.element=atoi(word[++i]);
+				fp.d_params.brand=atoi(word[++i]);
+				// parsed in future
+		fp.d_params.attr=ATOM_ATTR_HB|ATOM_ATTR_VA;
+		fp.d_params.attr|=ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP;
+				break;
+
+						case '1':
+				fp.d_params.type=DEFECT_TYPE_1D;
+				break;
+						case '2':
+				fp.d_params.type=DEFECT_TYPE_2D;
+				break;
+						case '3':
+				fp.d_params.type=DEFECT_TYPE_3D;
+				break;
+						default:
+							break;
+					}
+
+				}
+				i+=1;
 			}
-			if((fp.p_type==PART_INSIDE_R)||
-                           (fp.p_type==PART_OUTSIDE_R)) {
-				fp.p_vals.r=atof(word[7]);	
-				fp.p_vals.p.x=atof(word[8]);	
-				fp.p_vals.p.y=atof(word[9]);	
-				fp.p_vals.p.z=atof(word[10]);
-			}
-			if((fp.p_type==PART_INSIDE_D)||
-                           (fp.p_type==PART_OUTSIDE_D)) {
-				fp.p_vals.p.x=atof(word[7]);	
-				fp.p_vals.p.y=atof(word[8]);	
-				fp.p_vals.p.z=atof(word[9]);
-				fp.p_vals.d.x=atof(word[10]);	
-				fp.p_vals.d.y=atof(word[11]);	
-				fp.p_vals.d.z=atof(word[12]);
-			}
-			fp.lattice=mdrun->lattice;
 			add_stage(mdrun,STAGE_FILL,&fp);
 		}
 		else if(!strncmp(word[0],"thermal_init",12)) {
@@ -777,7 +825,7 @@ int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
 				}
 			}
 		}
-		add_atom(moldyn,iap->element,pse_mass[iap->element],
+		add_atom(moldyn,iap->element,
 		         iap->brand,iap->attr,&r,&v);
 		printf("%s atom inserted (%d/%d): %f %f %f\n",
 		       ME,(iap->cnt_steps+1)*iap->ins_atoms,
@@ -846,7 +894,7 @@ int insert_mixed_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
 					if(dmin>cmax)
 						retry=1;
 			}
-			add_atom(moldyn,imp->element1,pse_mass[imp->element1],
+			add_atom(moldyn,imp->element1,
                                  imp->brand1,imp->attr1,&r,&v);
 			printf("%s (mixed) atom inserted (%d): %f %f %f\n",
                                ME,imp->amount1,r.x,r.y,r.z);
@@ -880,7 +928,7 @@ int insert_mixed_atoms(t_moldyn *moldyn,t_mdrun *mdrun) {
 					if(dmin>cmax)
 						retry=1;
 			}
-			add_atom(moldyn,imp->element2,pse_mass[imp->element2],
+			add_atom(moldyn,imp->element2,
                                  imp->brand2,imp->attr2,&r,&v);
 			printf("%s (mixed) atom inserted (%d): %f %f %f\n",
                                ME,imp->amount2,r.x,r.y,r.z);
@@ -1103,39 +1151,44 @@ int mdrun_hook(void *ptr1,void *ptr2) {
 			case STAGE_FILL:
 				stage_print("  -> fill lattice\n\n");
 				fp=stage->params;
-				if(fp->lattice!=ZINCBLENDE) {
-					create_lattice(moldyn,
-					               fp->lattice,fp->lc,
-					               mdrun->element1,
-					               mdrun->m1,
-					               DEFAULT_ATOM_ATTR,0,
-					               fp->lx,fp->ly,fp->lz,
-					               NULL,fp->p_type,
-					               &(fp->p_vals));
-				}
-				else {
+				switch(fp->lattice) {
+					case ZINCBLENDE:
+
 					o.x=0.5*0.25*fp->lc;
 					o.y=o.x;
 					o.z=o.x;
 					create_lattice(moldyn,
 					               FCC,fp->lc,
 					               mdrun->element1,
-					               mdrun->m1,
 					               DEFAULT_ATOM_ATTR,0,
 					               fp->lx,fp->ly,fp->lz,
-					               &o,fp->p_type,
-					               &(fp->p_vals));
+					               &o,
+					               &(fp->p_params),
+					               &(fp->d_params));
 					o.x+=0.25*fp->lc;
 					o.y=o.x;
 					o.z=o.x;
 					create_lattice(moldyn,
 					               FCC,fp->lc,
 					               mdrun->element2,
-					               mdrun->m2,
 					               DEFAULT_ATOM_ATTR,1,
 					               fp->lx,fp->ly,fp->lz,
-					               &o,fp->p_type,
-					               &(fp->p_vals));
+					               &o,
+					               &(fp->p_params),
+					               &(fp->d_params));
+					break;
+
+					default:
+
+					create_lattice(moldyn,
+					               fp->lattice,fp->lc,
+					               mdrun->element1,
+					               DEFAULT_ATOM_ATTR,0,
+					               fp->lx,fp->ly,fp->lz,
+					               NULL,
+					               &(fp->p_params),
+					               &(fp->d_params));
+					break;
 				}
 				moldyn_bc_check(moldyn);
 				change_stage=TRUE;
diff --git a/mdrun.h b/mdrun.h
index d6476e5..fad1f62 100644
--- a/mdrun.h
+++ b/mdrun.h
@@ -71,9 +71,7 @@ typedef struct s_mdrun {
 	u8 pbcz;
 
 	int element1;				// element 1
-	double m1;
 	int element2;				// element 2
-	double m2;
 
 	double lc;				// lattice constant
 	u8 lattice;				// type of lattice
@@ -199,8 +197,8 @@ typedef struct s_fill_params {
 	u8 lattice;
 	int fill_element;
 	u8 fill_brand;
-	u8 p_type;
-	t_part_vals p_vals;
+	t_part_params p_params;
+	t_defect_params d_params;
 } t_fill_params;
 
 /*
diff --git a/moldyn.c b/moldyn.c
index a934194..bc4fdff 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -41,9 +41,11 @@
 #endif
 
 /* pse */
+#define PSE_MASS
 #define PSE_NAME
 #define PSE_COL
 #include "pse.h"
+#undef PSE_MASS
 #undef PSE_NAME
 #undef PSE_COL
 
@@ -512,9 +514,9 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
  * creating lattice functions
  */
 
-int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
+int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin,
-                   u8 p_type,t_part_vals *p_vals) {
+                   t_part_params *p_params,t_defect_params *d_params) {
 
 	int new,count;
 	int ret;
@@ -538,6 +540,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	if(type==FCC) new*=4;
 	if(type==DIAMOND) new*=8;
 
+	/* defects */
+	if(d_params->type) {
+		switch(d_params->stype) {
+			case DEFECT_STYPE_DB_X:
+			case DEFECT_STYPE_DB_Y:
+			case DEFECT_STYPE_DB_Z:
+			case DEFECT_STYPE_DB_R:
+				new*=2;
+				break;
+			default:
+				printf("[moldyn] WARNING: cl unknown defect\n");
+				break;
+		}
+	}
+
 	/* allocate space for atoms */
 	ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom));
 	if(!ptr) {
@@ -572,21 +589,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	switch(type) {
 		case CUBIC:
 			set_nn_dist(moldyn,lc);
-			ret=cubic_init(a,b,c,lc,atom,&orig,p_type,p_vals);
+			ret=cubic_init(a,b,c,lc,atom,&orig,p_params,d_params);
 			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,p_type,p_vals);
+			ret=fcc_init(a,b,c,lc,atom,&orig,p_params,d_params);
 			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,p_type,p_vals);
+			ret=diamond_init(a,b,c,lc,atom,&orig,p_params,d_params);
 			strcpy(name,"diamond");
 			break;
 		default:
@@ -598,7 +615,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	if(ret!=new) {
 		printf("[moldyn] creating %s lattice (lc=%f) incomplete\n",
 		       name,lc);
-		printf("  (ignore in case of partial lattice creation)\n");
+		printf("  (ignore for partial lattice creation)\n");
 		printf("  amount of atoms\n");
 		printf("  - expected: %d\n",new);
 		printf("  - created: %d\n",ret);
@@ -610,7 +627,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 
 	for(new=0;new<ret;new++) {
 		atom[new].element=element;
-		atom[new].mass=mass;
+		atom[new].mass=pse_mass[element];
 		atom[new].attr=attr;
 		atom[new].brand=brand;
 		atom[new].tag=count+new;
@@ -619,6 +636,19 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 #ifdef PTHREADS
 		pthread_mutex_init(&(mutex[new]),NULL);
 #endif
+		if(d_params->type) {
+			new+=1;
+			atom[new].element=d_params->element;
+			atom[new].mass=pse_mass[d_params->element];
+			atom[new].attr=d_params->attr;
+			atom[new].brand=d_params->brand;
+			atom[new].tag=count+new;
+			check_per_bound(moldyn,&(atom[new].r));
+			atom[new].r_0=atom[new].r;
+#ifdef PTHREADS
+			pthread_mutex_init(&(mutex[new]),NULL);
+#endif
+		}
 	}
 
 	/* fix allocation */
@@ -629,6 +659,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	}
 	moldyn->atom=ptr;
 
+// WHAT ABOUT AMUTEX !!!!
+
 #ifdef LOWMEM_LISTS
 	ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int));
 	if(!ptr) {
@@ -644,7 +676,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	return ret;
 }
 
-int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
+int add_atom(t_moldyn *moldyn,int element,u8 brand,u8 attr,
              t_3dvec *r,t_3dvec *v) {
 
 	t_atom *atom;
@@ -687,7 +719,7 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 	atom[count].r=*r;
 	atom[count].v=*v;
 	atom[count].element=element;
-	atom[count].mass=mass;
+	atom[count].mass=pse_mass[element];
 	atom[count].brand=brand;
 	atom[count].tag=count;
 	atom[count].attr=attr;
@@ -729,9 +761,44 @@ int del_atom(t_moldyn *moldyn,int tag) {
 	return 0;
 }
 
+#define	set_atom_positions(pos) \
+	if(d_params->type) {\
+		d_o.x=0; d_o.y=0; d_o.z=0;\
+		d_d.x=0; d_d.y=0; d_d.z=0;\
+		switch(d_params->stype) {\
+			case DEFECT_STYPE_DB_X:\
+				d_o.x=d_params->od;\
+				d_d.x=d_params->dd;\
+				break;\
+			case DEFECT_STYPE_DB_Y:\
+				d_o.y=d_params->od;\
+				d_d.y=d_params->dd;\
+				break;\
+			case DEFECT_STYPE_DB_Z:\
+				d_o.z=d_params->od;\
+				d_d.z=d_params->dd;\
+				break;\
+			case DEFECT_STYPE_DB_R:\
+				break;\
+			default:\
+				printf("[moldyn] WARNING: unknown defect\n");\
+				break;\
+		}\
+		v3_add(&dr,&pos,&d_o);\
+		v3_copy(&(atom[count].r),&dr);\
+		count+=1;\
+		v3_add(&dr,&pos,&d_d);\
+		v3_copy(&(atom[count].r),&dr);\
+		count+=1;\
+	}\
+	else {\
+		v3_copy(&(atom[count].r),&pos);\
+		count+=1;\
+	}
+
 /* cubic init */
 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
-               u8 p_type,t_part_vals *p_vals) {
+               t_part_params *p_params,t_defect_params *d_params) {
 
 	int count;
 	t_3dvec r;
@@ -739,6 +806,9 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 	t_3dvec o;
 	t_3dvec dist;
 	t_3dvec p;
+	t_3dvec d_o;
+	t_3dvec d_d;
+	t_3dvec dr;
 
 	p.x=0; p.y=0; p.z=0;
 
@@ -749,10 +819,10 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 		v3_zero(&o);
 
 	/* shift partition values */
-	if(p_type) {
-		p.x=p_vals->p.x+(a*lc)/2.0;
-		p.y=p_vals->p.y+(b*lc)/2.0;
-		p.z=p_vals->p.z+(c*lc)/2.0;
+	if(p_params->type) {
+		p.x=p_params->p.x+(a*lc)/2.0;
+		p.y=p_params->p.y+(b*lc)/2.0;
+		p.z=p_params->p.z+(c*lc)/2.0;
 	}
 
 	r.x=o.x;
@@ -761,42 +831,39 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 		for(j=0;j<b;j++) {
 			r.z=o.z;
 			for(k=0;k<c;k++) {
-				switch(p_type) {
+				switch(p_params->type) {
 					case PART_INSIDE_R:
 						v3_sub(&dist,&r,&p);
-			if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
-				v3_copy(&(atom[count].r),&r);
-				count+=1;
+			if(v3_absolute_square(&dist)<
+                           (p_params->r*p_params->r)) {
+				set_atom_positions(r);
 			}
 						break;
 					case PART_OUTSIDE_R:
 						v3_sub(&dist,&r,&p);
-			if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
-				v3_copy(&(atom[count].r),&r);
-				count+=1;
+			if(v3_absolute_square(&dist)>=
+                           (p_params->r*p_params->r)) {
+				set_atom_positions(r);
 			}
 						break;
 					case PART_INSIDE_D:
 						v3_sub(&dist,&r,&p);
-			if((fabs(dist.x)<p_vals->d.x)&&
-			   (fabs(dist.y)<p_vals->d.y)&&
-			   (fabs(dist.z)<p_vals->d.z)) {
-				v3_copy(&(atom[count].r),&r);
-				count+=1;
+			if((fabs(dist.x)<p_params->d.x)&&
+			   (fabs(dist.y)<p_params->d.y)&&
+			   (fabs(dist.z)<p_params->d.z)) {
+				set_atom_positions(r);
 			}
 						break;
 					case PART_OUTSIDE_D:
 						v3_sub(&dist,&r,&p);
-			if((fabs(dist.x)>=p_vals->d.x)||
-			   (fabs(dist.y)>=p_vals->d.y)||
-			   (fabs(dist.z)>=p_vals->d.z)) {
-				v3_copy(&(atom[count].r),&r);
-				count+=1;
+			if((fabs(dist.x)>=p_params->d.x)||
+			   (fabs(dist.y)>=p_params->d.y)||
+			   (fabs(dist.z)>=p_params->d.z)) {
+				set_atom_positions(r);
 			}
 						break;
 					default:	
-						v3_copy(&(atom[count].r),&r);
-						count+=1;
+						set_atom_positions(r);
 						break;
 				}
 				r.z+=lc;
@@ -817,7 +884,7 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 
 /* fcc lattice init */
 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
-             u8 p_type,t_part_vals *p_vals) {
+             t_part_params *p_params,t_defect_params *d_params) {
 
 	int count;
 	int i,j,k,l;
@@ -825,6 +892,7 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 	t_3dvec basis[3];
 	t_3dvec dist;
 	t_3dvec p;
+	t_3dvec d_d,d_o,dr;
 
 	p.x=0; p.y=0; p.z=0;
 
@@ -844,10 +912,10 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 	basis[2].z=0.5*lc;
 
 	/* shift partition values */
-	if(p_type) {
-		p.x=p_vals->p.x+(a*lc)/2.0;
-		p.y=p_vals->p.y+(b*lc)/2.0;
-		p.z=p_vals->p.z+(c*lc)/2.0;
+	if(p_params->type) {
+		p.x=p_params->p.x+(a*lc)/2.0;
+		p.y=p_params->p.y+(b*lc)/2.0;
+		p.z=p_params->p.z+(c*lc)/2.0;
 	}
 
 	/* fill up the room */
@@ -858,83 +926,77 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 			r.z=o.z;
 			for(k=0;k<c;k++) {
 				/* first atom */
-				switch(p_type) {
+				switch(p_params->type) {
 					case PART_INSIDE_R:
 						v3_sub(&dist,&r,&p);
-			if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
-				v3_copy(&(atom[count].r),&r);
-				count+=1;
+			if(v3_absolute_square(&dist)<
+			   (p_params->r*p_params->r)) {
+				set_atom_positions(r);
 			}
 						break;
 					case PART_OUTSIDE_R:
 						v3_sub(&dist,&r,&p);
-			if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
-				v3_copy(&(atom[count].r),&r);
-				count+=1;
+			if(v3_absolute_square(&dist)>=
+			   (p_params->r*p_params->r)) {
+				set_atom_positions(r);
 			}
 						break;
 					case PART_INSIDE_D:
 						v3_sub(&dist,&r,&p);
-			if((fabs(dist.x)<p_vals->d.x)&&
-			   (fabs(dist.y)<p_vals->d.y)&&
-			   (fabs(dist.z)<p_vals->d.z)) {
-				v3_copy(&(atom[count].r),&r);
-				count+=1;
+			if((fabs(dist.x)<p_params->d.x)&&
+			   (fabs(dist.y)<p_params->d.y)&&
+			   (fabs(dist.z)<p_params->d.z)) {
+				set_atom_positions(r);
 			}
 						break;
 					case PART_OUTSIDE_D:
 						v3_sub(&dist,&r,&p);
-			if((fabs(dist.x)>=p_vals->d.x)||
-			   (fabs(dist.y)>=p_vals->d.y)||
-			   (fabs(dist.z)>=p_vals->d.z)) {
-				v3_copy(&(atom[count].r),&r);
-				count+=1;
+			if((fabs(dist.x)>=p_params->d.x)||
+			   (fabs(dist.y)>=p_params->d.y)||
+			   (fabs(dist.z)>=p_params->d.z)) {
+				set_atom_positions(r);
 			}
 						break;
 					default:
-						v3_copy(&(atom[count].r),&r);
-						count+=1;
+						set_atom_positions(r);
 						break;
 				}
 				/* the three face centered atoms */
 				for(l=0;l<3;l++) {
 					v3_add(&n,&r,&basis[l]);
-					switch(p_type) {
+					switch(p_params->type) {
 						case PART_INSIDE_R:
 			v3_sub(&dist,&n,&p);
-			if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
-				v3_copy(&(atom[count].r),&n);
-				count+=1;
+			if(v3_absolute_square(&dist)<
+			   (p_params->r*p_params->r)) {
+				set_atom_positions(n);
 			}
 							break;
 						case PART_OUTSIDE_R:
 			v3_sub(&dist,&n,&p);
-			if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
-				v3_copy(&(atom[count].r),&n);
-				count+=1;
+			if(v3_absolute_square(&dist)>=
+			   (p_params->r*p_params->r)) {
+				set_atom_positions(n);
 			}
 							break;
 					case PART_INSIDE_D:
 						v3_sub(&dist,&n,&p);
-			if((fabs(dist.x)<p_vals->d.x)&&
-			   (fabs(dist.y)<p_vals->d.y)&&
-			   (fabs(dist.z)<p_vals->d.z)) {
-				v3_copy(&(atom[count].r),&n);
-				count+=1;
+			if((fabs(dist.x)<p_params->d.x)&&
+			   (fabs(dist.y)<p_params->d.y)&&
+			   (fabs(dist.z)<p_params->d.z)) {
+				set_atom_positions(n);
 			}
 						break;
 					case PART_OUTSIDE_D:
 						v3_sub(&dist,&n,&p);
-			if((fabs(dist.x)>=p_vals->d.x)||
-			   (fabs(dist.y)>=p_vals->d.y)||
-			   (fabs(dist.z)>=p_vals->d.z)) {
-				v3_copy(&(atom[count].r),&n);
-				count+=1;
+			if((fabs(dist.x)>=p_params->d.x)||
+			   (fabs(dist.y)>=p_params->d.y)||
+			   (fabs(dist.z)>=p_params->d.z)) {
+				set_atom_positions(n);
 			}
 						break;
 						default:
-					v3_copy(&(atom[count].r),&n);
-					count+=1;
+							set_atom_positions(n);
 							break;
 					}
 				}
@@ -956,12 +1018,12 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 }
 
 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
-                 u8 p_type,t_part_vals *p_vals) {
+                 t_part_params *p_params,t_defect_params *d_params) {
 
 	int count;
 	t_3dvec o;
 
-	count=fcc_init(a,b,c,lc,atom,origin,p_type,p_vals);
+	count=fcc_init(a,b,c,lc,atom,origin,p_params,d_params);
 
 	o.x=0.25*lc;
 	o.y=0.25*lc;
@@ -969,7 +1031,7 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 
 	if(origin) v3_add(&o,&o,origin);
 
-	count+=fcc_init(a,b,c,lc,&atom[count],&o,p_type,p_vals);
+	count+=fcc_init(a,b,c,lc,&atom[count],&o,p_params,d_params);
 
 	return count;
 }
@@ -1171,12 +1233,12 @@ double pressure_calc(t_moldyn *moldyn) {
 	moldyn->p=2.0*moldyn->ekin+moldyn->virial;
 	moldyn->p/=(3.0*moldyn->volume);
 
-	moldyn->px=2.0*moldyn->ekinx+moldyn->vir.xx;
-	moldyn->px/=moldyn->volume;
-	moldyn->py=2.0*moldyn->ekiny+moldyn->vir.yy;
-	moldyn->py/=moldyn->volume;
-	moldyn->pz=2.0*moldyn->ekinz+moldyn->vir.zz;
-	moldyn->pz/=moldyn->volume;
+	//moldyn->px=2.0*moldyn->ekinx+moldyn->vir.xx;
+	//moldyn->px/=moldyn->volume;
+	//moldyn->py=2.0*moldyn->ekiny+moldyn->vir.yy;
+	//moldyn->py/=moldyn->volume;
+	//moldyn->pz=2.0*moldyn->ekinz+moldyn->vir.zz;
+	//moldyn->pz/=moldyn->volume;
 
 	/* pressure (absolute coordinates) */
 	//moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
@@ -1439,16 +1501,15 @@ int scale_dim_ind(t_moldyn *moldyn,double x,double y,double z) {
 int scale_volume(t_moldyn *moldyn) {
 
 	t_3dvec *dim,*vdim;
-	//double scale;
+	double scale;
 	t_linkcell *lc;
-	double sx,sy,sz;
+	//double sx,sy,sz;
 
 	vdim=&(moldyn->vis.dim);
 	dim=&(moldyn->dim);
 	lc=&(moldyn->lc);
 
 	/* scaling factor */
-	/*
 	if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
 		scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau;
 		scale=pow(scale,ONE_THIRD);
@@ -1456,20 +1517,22 @@ int scale_volume(t_moldyn *moldyn) {
 	else {
 		scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
 	}
-	*/
 
+
+	/*	
 	sx=1.0-(moldyn->p_ref-moldyn->px)*moldyn->p_tc*moldyn->tau;
 	sy=1.0-(moldyn->p_ref-moldyn->py)*moldyn->p_tc*moldyn->tau;
 	sz=1.0-(moldyn->p_ref-moldyn->pz)*moldyn->p_tc*moldyn->tau;
 	sx=pow(sx,ONE_THIRD);
 	sy=pow(sy,ONE_THIRD);
 	sz=pow(sz,ONE_THIRD);
+	*/
 
 	/* scale the atoms and dimensions */
-	//scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
-	//scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
-	scale_atoms_ind(moldyn,sx,sy,sz);
-	scale_dim_ind(moldyn,sx,sy,sz);
+	scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
+	scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
+	//scale_atoms_ind(moldyn,sx,sy,sz);
+	//scale_dim_ind(moldyn,sx,sy,sz);
 
 	/* visualize dimensions */
 	if(vdim->x!=0) {
@@ -1488,12 +1551,12 @@ int scale_volume(t_moldyn *moldyn) {
 		link_cell_shutdown(moldyn);
 		link_cell_init(moldyn,QUIET);
 	} else {
-		//lc->x*=scale;
-		//lc->y*=scale;
-		//lc->z*=scale;
-		lc->x*=sx;
-		lc->y*=sx;
-		lc->z*=sy;
+		lc->x*=scale;
+		lc->y*=scale;
+		lc->z*=scale;
+		//lc->x*=sx;
+		//lc->y*=sx;
+		//lc->z*=sy;
 	}
 
 	return 0;
@@ -1507,16 +1570,16 @@ double e_kin_calc(t_moldyn *moldyn) {
 
 	atom=moldyn->atom;
 	moldyn->ekin=0.0;
-	moldyn->ekinx=0.0;
-	moldyn->ekiny=0.0;
-	moldyn->ekinz=0.0;
+	//moldyn->ekinx=0.0;
+	//moldyn->ekiny=0.0;
+	//moldyn->ekinz=0.0;
 
 	for(i=0;i<moldyn->count;i++) {
 		atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 		moldyn->ekin+=atom[i].ekin;
-		moldyn->ekinx+=0.5*atom[i].mass*atom[i].v.x*atom[i].v.x;
-		moldyn->ekiny+=0.5*atom[i].mass*atom[i].v.y*atom[i].v.y;
-		moldyn->ekinz+=0.5*atom[i].mass*atom[i].v.z*atom[i].v.z;
+		//moldyn->ekinx+=0.5*atom[i].mass*atom[i].v.x*atom[i].v.x;
+		//moldyn->ekiny+=0.5*atom[i].mass*atom[i].v.y*atom[i].v.y;
+		//moldyn->ekinz+=0.5*atom[i].mass*atom[i].v.z*atom[i].v.z;
 	}
 
 	return moldyn->ekin;
@@ -3264,6 +3327,7 @@ int visual_atoms(t_moldyn *moldyn) {
 	t_visual *v;
 	t_atom *atom;
 	t_vb vb;
+	t_3dvec strain;
 #ifdef VISUAL_THREAD
 	t_moldyn *moldyn;
 
@@ -3294,14 +3358,18 @@ int visual_atoms(t_moldyn *moldyn) {
 	        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++) {
+		v3_sub(&strain,&(atom[i].r),&(atom[i].r_0));
+		check_per_bound(moldyn,&strain);
 		// atom type, positions, color and kinetic energy
 		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);
+		                                    //atom[i].ekin);
+		                                    sqrt(v3_absolute_square(&strain)));
+	}
 	
 	// bonds between atoms
 #ifndef VISUAL_THREAD
diff --git a/moldyn.h b/moldyn.h
index d065d43..538e0eb 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -258,17 +258,38 @@ typedef struct s_vb {
 	int fd;
 } t_vb;
 
-typedef struct s_part_vals {
+typedef struct s_part_params {
+	u8 type;
 	double r;
 	t_3dvec p;
 	t_3dvec d;
-} t_part_vals;
+} t_part_params;
 
 #define PART_INSIDE_R   1
 #define PART_OUTSIDE_R  2
 #define PART_INSIDE_D   3
 #define PART_OUTSIDE_D  4
 
+typedef struct s_defect_params {
+	u8 type;
+	u8 stype;
+	double od;
+	double dd;
+	int element;
+	u8 brand;
+	u8 attr;
+} t_defect_params;
+
+#define DEFECT_TYPE_0D	1
+#define DEFECT_TYPE_1D	2
+#define DEFECT_TYPE_2D	3
+#define DEFECT_TYPE_3D	4
+
+#define DEFECT_STYPE_DB_X	1
+#define DEFECT_STYPE_DB_Y	2
+#define DEFECT_STYPE_DB_Z	3
+#define DEFECT_STYPE_DB_R	4
+
 /*
  *
  *  defines
@@ -415,18 +436,18 @@ int moldyn_set_report(t_moldyn *moldyn,char *author,char *title);
 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer);
 int moldyn_log_shutdown(t_moldyn *moldyn);
 
-int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
+int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin,
-                   u8 p_type,t_part_vals *p_vals);
-int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
+                   t_part_params *p_params,t_defect_params *d_params);
+int add_atom(t_moldyn *moldyn,int element,u8 brand,u8 attr,
              t_3dvec *r,t_3dvec *v);
 int del_atom(t_moldyn *moldyn,int tag);
 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
-               u8 p_type,t_part_vals *p_vals);
+               t_part_params *p_params,t_defect_params *d_params);
 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
-             u8 p_type,t_part_vals *p_vals);
+             t_part_params *p_params,t_defect_params *d_params);
 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
-                 u8 p_type,t_part_vals *p_vals);
+                 t_part_params *p_params,t_defect_params *d_params);
 int destroy_atoms(t_moldyn *moldyn);
 
 int thermal_init(t_moldyn *moldyn,u8 equi_init);
diff --git a/remember_me.txt b/remember_me.txt
index 7247fbb..3ae541e 100644
--- a/remember_me.txt
+++ b/remember_me.txt
@@ -11,3 +11,13 @@ c-c
 
 3.1: saves/c_in_si_prec_450_01/s-0566000.save 239458 240660 110 - 100
                                               239703 241605 100 - 100
+
+
+concatenated, differently oriented 100 dumbbells
+################################################
+
+posic_new/saves/c_in_si_prec_450_tot_02/s-0272000.save b 3.08 0.01
+
+ atoms 238338/1 238631/1 - 3.078045
+ atoms 238461/1 244278/1 - 3.085715
+
diff --git a/sic.c b/sic.c
index f769f6f..375c6d5 100644
--- a/sic.c
+++ b/sic.c
@@ -121,7 +121,7 @@ int insert_atoms(t_moldyn *moldyn) {
 					dmin=d;
 			}
 		}
-		add_atom(moldyn,INS_TYPE,INS_MASS,INS_BRAND,
+		add_atom(moldyn,INS_TYPE,INS_BRAND,
 		         ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|\
 			 INS_ATTR,
 		         &r,&v);
@@ -325,23 +325,23 @@ int main(int argc,char **argv) {
 	// diamond
 #ifdef ALBE
  #ifdef INIT_SI
-	create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,M_SI,
+	create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,
 	               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 		       0,LCNTX,LCNTY,LCNTZ,NULL,0,NULL);
  #endif
  #ifdef INIT_C
-	create_lattice(&md,DIAMOND,ALBE_LC_C,C,M_C,
+	create_lattice(&md,DIAMOND,ALBE_LC_C,C,
 	               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 		       1,LCNTX,LCNTY,LCNTZ,NULL,0,NULL);
  #endif
 #else
  #ifdef INIT_SI
-	create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+	create_lattice(&md,DIAMOND,LC_SI,SI,
 	               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 		       0,LCNTX,LCNTY,LCNTZ,NULL,0,NULL);
  #endif
  #ifdef INIT_C
-	create_lattice(&md,DIAMOND,LC_C,SI,M_SI,
+	create_lattice(&md,DIAMOND,LC_C,SI,
 	               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 		       1,LCNTX,LCNTY,LCNTZ,NULL,0,NULL);
  #endif
@@ -351,20 +351,20 @@ int main(int argc,char **argv) {
 #ifdef INIT_3CSIC
  #ifdef ALBE
 	r.x=0.5*0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
-	create_lattice(&md,FCC,ALBE_LC_SIC,SI,M_SI,
+	create_lattice(&md,FCC,ALBE_LC_SIC,SI,
 	               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 	               0,LCNTX,LCNTY,LCNTZ,&r,0,NULL);
 	r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
-	create_lattice(&md,FCC,ALBE_LC_SIC,C,M_C,
+	create_lattice(&md,FCC,ALBE_LC_SIC,C,
 	               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB|ATOM_ATTR_VB,
 	               1,LCNTX,LCNTY,LCNTZ,&r,0,NULL);
  #else
 	r.x=0.5*0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
-	create_lattice(&md,FCC,TM_LC_SIC,SI,M_SI,
+	create_lattice(&md,FCC,TM_LC_SIC,SI,
 	               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 	               0,LCNTX,LCNTY,LCNTZ,&r,0,NULL);
 	r.x+=0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
-	create_lattice(&md,FCC,TM_LC_SIC,C,M_C,
+	create_lattice(&md,FCC,TM_LC_SIC,C,
 	               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 	               1,LCNTX,LCNTY,LCNTZ,&r,0,NULL);
  #endif
diff --git a/visualize b/visualize
index 665cb67..0b8fc61 100755
--- a/visualize
+++ b/visualize
@@ -39,6 +39,7 @@ bx0=""; by0=""; bz0="";
 bx1=""; by1=""; bz1="";
 bcr="";
 clx="0"; cly="0"; clz="0";
+extra=0
 
 # parse argv
 
@@ -58,6 +59,7 @@ while [ "$1" ]; do
 				bx1=$5; by1=$6; bz1=$7;	shift 7;;
 		-B)		bcr=$2;			shift 2;;
 		-C)		lc=$2;			shift 2;;
+		-e)		extra=1; 		shift 1;;
 		*)
 				echo "options:"
 				echo "########"
@@ -159,18 +161,21 @@ EOF
 
 	# atoms
 	if [ -n "$x0" ]; then
-		export x0 y0 z0 x1 y1 z1 radius
+		export x0 y0 z0 x1 y1 z1 radius extra
 		cat $file | grep -v '#' | awk '\
 		BEGIN {
 			x0=ENVIRON["x0"]; y0=ENVIRON["y0"]; z0=ENVIRON["z0"];
 			x1=ENVIRON["x1"]; y1=ENVIRON["y1"]; z1=ENVIRON["z1"];
-			radius=ENVIRON["radius"];
+			radius=ENVIRON["radius"]; extra=ENVIRON["extra"];
 		}
 		{
 			if(($2>=x0)&&($3>=y0)&&($4>=z0)&&\
 			   ($2<=x1)&&($3<=y1)&&($4<=z1)) {
 				print "sphere { <"$2","$4","$3">, "radius" ";
-				print "texture { pigment { color "$5" } ";
+				if(extra)
+		print "texture { pigment { color rgb<"$6/4.4",0,"1-$6/4.4"> } ";
+				else
+		print "texture { pigment { color "$5" } ";
 				print "finish { phong 1 metallic } } }";
 			}
 		}' >> temp.pov
@@ -194,29 +199,30 @@ EOF
 	# boundaries
 	if [ -z "$bx0" ]; then
 
-	if [ -z "$x0" ]; then
+	#if [ -z "$x0" ]; then
 
 	cat $file | grep '# \[D\]' | while read foo bar x1 y1 z1 x2 y2 z2 ; do
-		draw_cyl $x1 $z1 $y1 $x2 $z2 $y2 0.05
+		draw_cyl $x1 $y1 $z1 $x2 $y2 $z2 0.05
 	done
 
-	else
+	#else
+
 		# manually drawing the 3x4 boundaries ...
-		draw_cyl $x0 $y0 $z0 $x1 $y0 $z0
-		draw_cyl $x0 $y0 $z0 $x0 $y1 $z0
-		draw_cyl $x1 $y1 $z0 $x1 $y0 $z0
-		draw_cyl $x0 $y1 $z0 $x1 $y1 $z0
-
-		draw_cyl $x0 $y0 $z1 $x1 $y0 $z1
-		draw_cyl $x0 $y0 $z1 $x0 $y1 $z1
-		draw_cyl $x1 $y1 $z1 $x1 $y0 $z1
-		draw_cyl $x0 $y1 $z1 $x1 $y1 $z1
-
-		draw_cyl $x0 $y0 $z1 $x0 $y0 $z0
-		draw_cyl $x0 $y1 $z1 $x0 $y1 $z0
-		draw_cyl $x1 $y0 $z1 $x1 $y0 $z0
-		draw_cyl $x1 $y1 $z1 $x1 $y1 $z0
-	fi
+#		draw_cyl $x0 $y0 $z0 $x1 $y0 $z0
+#		draw_cyl $x0 $y0 $z0 $x0 $y1 $z0
+#		draw_cyl $x1 $y1 $z0 $x1 $y0 $z0
+#		draw_cyl $x0 $y1 $z0 $x1 $y1 $z0
+
+#		draw_cyl $x0 $y0 $z1 $x1 $y0 $z1
+#		draw_cyl $x0 $y0 $z1 $x0 $y1 $z1
+#		draw_cyl $x1 $y1 $z1 $x1 $y0 $z1
+#		draw_cyl $x0 $y1 $z1 $x1 $y1 $z1
+
+#		draw_cyl $x0 $y0 $z1 $x0 $y0 $z0
+#		draw_cyl $x0 $y1 $z1 $x0 $y1 $z0
+#		draw_cyl $x1 $y0 $z1 $x1 $y0 $z0
+#		draw_cyl $x1 $y1 $z1 $x1 $y1 $z0
+#	fi
 
 	else
 
@@ -287,7 +293,7 @@ EOF
 
 	# mv png
 	$POVRAY temp.pov > /dev/null 2>&1
-	time=`echo $file | awk -F. '{ print $1 }' | awk -F_ '{ print $3 }'`
+	time=`basename $file | awk -F. '{ print $1 }' | awk -F_ '{ print $3 }'`
 	convert $COPTS -draw "text 5,20 't = $time fs'" temp.png temp.png
 	mv temp.png `echo $file | sed 's/\.xyz/\.png/'`