X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=bc4fdffdee29be7341c67692f1a3282bc049927b;hb=1965279af348fbb9b69459d01516bbcd52b6f240;hp=ed30815a6acb2b0d9cf151f990a18d3579067483;hpb=815b8841eb76564b2f90421d4eae180253453106;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index ed30815..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
 
@@ -90,6 +92,7 @@ int moldyn_shutdown(t_moldyn *moldyn) {
 #ifdef PTHREADS
 	for(i=0;i<moldyn->count;i++)
 		pthread_mutex_destroy(&(amutex[i]));
+	free(amutex);
 	pthread_mutex_destroy(&emutex);
 #endif
 
@@ -511,8 +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,
-                   u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) {
+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,
+                   t_part_params *p_params,t_defect_params *d_params) {
 
 	int new,count;
 	int ret;
@@ -536,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) {
@@ -570,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);
+			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);
+			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);
+			ret=diamond_init(a,b,c,lc,atom,&orig,p_params,d_params);
 			strcpy(name,"diamond");
 			break;
 		default:
@@ -594,36 +613,70 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 
 	/* debug */
 	if(ret!=new) {
-		printf("[moldyn] creating lattice failed\n");
+		printf("[moldyn] creating %s lattice (lc=%f) incomplete\n",
+		       name,lc);
+		printf("  (ignore for partial lattice creation)\n");
 		printf("  amount of atoms\n");
 		printf("  - expected: %d\n",new);
 		printf("  - created: %d\n",ret);
-		return -1;
 	}
 
-	moldyn->count+=new;
-	printf("[moldyn] created %s lattice with %d atoms\n",name,new);
+	moldyn->count+=ret;
+	if(ret==new)
+		printf("[moldyn] created %s lattice with %d atoms\n",name,ret);
 
-	for(ret=0;ret<new;ret++) {
-		atom[ret].element=element;
-		atom[ret].mass=mass;
-		atom[ret].attr=attr;
-		atom[ret].brand=brand;
-		atom[ret].tag=count+ret;
-		check_per_bound(moldyn,&(atom[ret].r));
-		atom[ret].r_0=atom[ret].r;
+	for(new=0;new<ret;new++) {
+		atom[new].element=element;
+		atom[new].mass=pse_mass[element];
+		atom[new].attr=attr;
+		atom[new].brand=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
+		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[ret]),NULL);
+			pthread_mutex_init(&(mutex[new]),NULL);
 #endif
+		}
 	}
 
+	/* fix allocation */
+	ptr=realloc(moldyn->atom,moldyn->count*sizeof(t_atom));
+	if(!ptr) {
+		perror("[moldyn] realloc (create lattice - alloc fix)");
+		return -1;
+	}
+	moldyn->atom=ptr;
+
+// WHAT ABOUT AMUTEX !!!!
+
+#ifdef LOWMEM_LISTS
+	ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int));
+	if(!ptr) {
+		perror("[moldyn] list realloc (create lattice)");
+		return -1;
+	}
+	moldyn->lc.subcell->list=ptr;
+#endif
+
 	/* update total system mass */
 	total_mass_calc(moldyn);
 
 	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;
@@ -666,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;
@@ -708,13 +761,56 @@ 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) {
+int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
+               t_part_params *p_params,t_defect_params *d_params) {
 
 	int count;
 	t_3dvec r;
 	int i,j,k;
 	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;
 
 	count=0;
 	if(origin)
@@ -722,14 +818,54 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 	else
 		v3_zero(&o);
 
+	/* shift partition values */
+	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;
 	for(i=0;i<a;i++) {
 		r.y=o.y;
 		for(j=0;j<b;j++) {
 			r.z=o.z;
 			for(k=0;k<c;k++) {
-				v3_copy(&(atom[count].r),&r);
-				count+=1;
+				switch(p_params->type) {
+					case PART_INSIDE_R:
+						v3_sub(&dist,&r,&p);
+			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_params->r*p_params->r)) {
+				set_atom_positions(r);
+			}
+						break;
+					case PART_INSIDE_D:
+						v3_sub(&dist,&r,&p);
+			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_params->d.x)||
+			   (fabs(dist.y)>=p_params->d.y)||
+			   (fabs(dist.z)>=p_params->d.z)) {
+				set_atom_positions(r);
+			}
+						break;
+					default:	
+						set_atom_positions(r);
+						break;
+				}
 				r.z+=lc;
 			}
 			r.y+=lc;
@@ -747,12 +883,18 @@ 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) {
+int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
+             t_part_params *p_params,t_defect_params *d_params) {
 
 	int count;
 	int i,j,k,l;
 	t_3dvec o,r,n;
 	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;
 
 	count=0;
 	if(origin)
@@ -769,6 +911,13 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 	basis[2].y=0.5*lc;
 	basis[2].z=0.5*lc;
 
+	/* shift partition values */
+	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 */
 	r.x=o.x;
 	for(i=0;i<a;i++) {
@@ -777,15 +926,81 @@ 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 */
-				v3_copy(&(atom[count].r),&r);
-				count+=1;
-				r.z+=lc;
+				switch(p_params->type) {
+					case PART_INSIDE_R:
+						v3_sub(&dist,&r,&p);
+			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_params->r*p_params->r)) {
+				set_atom_positions(r);
+			}
+						break;
+					case PART_INSIDE_D:
+						v3_sub(&dist,&r,&p);
+			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_params->d.x)||
+			   (fabs(dist.y)>=p_params->d.y)||
+			   (fabs(dist.z)>=p_params->d.z)) {
+				set_atom_positions(r);
+			}
+						break;
+					default:
+						set_atom_positions(r);
+						break;
+				}
 				/* the three face centered atoms */
 				for(l=0;l<3;l++) {
 					v3_add(&n,&r,&basis[l]);
-					v3_copy(&(atom[count].r),&n);
-					count+=1;
+					switch(p_params->type) {
+						case PART_INSIDE_R:
+			v3_sub(&dist,&n,&p);
+			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_params->r*p_params->r)) {
+				set_atom_positions(n);
+			}
+							break;
+					case PART_INSIDE_D:
+						v3_sub(&dist,&n,&p);
+			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_params->d.x)||
+			   (fabs(dist.y)>=p_params->d.y)||
+			   (fabs(dist.z)>=p_params->d.z)) {
+				set_atom_positions(n);
+			}
+						break;
+						default:
+							set_atom_positions(n);
+							break;
+					}
 				}
+				r.z+=lc;
 			}
 			r.y+=lc;
 		}
@@ -802,12 +1017,13 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 	return count;
 }
 
-int diamond_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,
+                 t_part_params *p_params,t_defect_params *d_params) {
 
 	int count;
 	t_3dvec o;
 
-	count=fcc_init(a,b,c,lc,atom,origin);
+	count=fcc_init(a,b,c,lc,atom,origin,p_params,d_params);
 
 	o.x=0.25*lc;
 	o.y=0.25*lc;
@@ -815,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);
+	count+=fcc_init(a,b,c,lc,&atom[count],&o,p_params,d_params);
 
 	return count;
 }
@@ -893,7 +1109,9 @@ double temperature_calc(t_moldyn *moldyn) {
 
 	/* assume up to date kinetic energy, which is 3/2 N k_B T */
 
-	moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
+	if(moldyn->count)
+		moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
+	else moldyn->t=0.0;
 
 	return moldyn->t;
 }
@@ -993,8 +1211,8 @@ double virial_sum(t_moldyn *moldyn) {
 	}
 
 	/* global virial (absolute coordinates) */
-	virial=&(moldyn->gvir);
-	moldyn->gv=virial->xx+virial->yy+virial->zz;
+	//virial=&(moldyn->gvir);
+	//moldyn->gv=virial->xx+virial->yy+virial->zz;
 
 	return moldyn->virial;
 }
@@ -1015,9 +1233,16 @@ 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;
+
 	/* pressure (absolute coordinates) */
-	moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
-	moldyn->gp/=(3.0*moldyn->volume);
+	//moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
+	//moldyn->gp/=(3.0*moldyn->volume);
 
 	return moldyn->p;
 }
@@ -1042,11 +1267,11 @@ int average_reset(t_moldyn *moldyn) {
 
 	/* virial */
 	moldyn->virial_sum=0.0;
-	moldyn->gv_sum=0.0;
+	//moldyn->gv_sum=0.0;
 
 	/* pressure */
 	moldyn->p_sum=0.0;
-	moldyn->gp_sum=0.0;
+	//moldyn->gp_sum=0.0;
 	moldyn->tp_sum=0.0;
 
 	return 0;
@@ -1084,14 +1309,14 @@ int average_and_fluctuation_calc(t_moldyn *moldyn) {
 	/* virial */
 	moldyn->virial_sum+=moldyn->virial;
 	moldyn->virial_avg=moldyn->virial_sum/denom;
-	moldyn->gv_sum+=moldyn->gv;
-	moldyn->gv_avg=moldyn->gv_sum/denom;
+	//moldyn->gv_sum+=moldyn->gv;
+	//moldyn->gv_avg=moldyn->gv_sum/denom;
 
 	/* pressure */
 	moldyn->p_sum+=moldyn->p;
 	moldyn->p_avg=moldyn->p_sum/denom;
-	moldyn->gp_sum+=moldyn->gp;
-	moldyn->gp_avg=moldyn->gp_sum/denom;
+	//moldyn->gp_sum+=moldyn->gp;
+	//moldyn->gp_avg=moldyn->gp_sum/denom;
 	moldyn->tp_sum+=moldyn->tp;
 	moldyn->tp_avg=moldyn->tp_sum/denom;
 
@@ -1245,11 +1470,40 @@ int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
 	return 0;
 }
 
+int scale_atoms_ind(t_moldyn *moldyn,double x,double y,double z) {
+
+	int i;
+	t_3dvec *r;
+
+	for(i=0;i<moldyn->count;i++) {
+		r=&(moldyn->atom[i].r);
+		r->x*=x;
+		r->y*=y;
+		r->z*=z;
+	}
+
+	return 0;
+}
+
+int scale_dim_ind(t_moldyn *moldyn,double x,double y,double z) {
+
+	t_3dvec *dim;
+
+	dim=&(moldyn->dim);
+
+	dim->x*=x;
+	dim->y*=y;
+	dim->z*=z;
+
+	return 0;
+}
+
 int scale_volume(t_moldyn *moldyn) {
 
 	t_3dvec *dim,*vdim;
 	double scale;
 	t_linkcell *lc;
+	//double sx,sy,sz;
 
 	vdim=&(moldyn->vis.dim);
 	dim=&(moldyn->dim);
@@ -1264,9 +1518,21 @@ int scale_volume(t_moldyn *moldyn) {
 		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);
 
 	/* visualize dimensions */
 	if(vdim->x!=0) {
@@ -1288,6 +1554,9 @@ int scale_volume(t_moldyn *moldyn) {
 		lc->x*=scale;
 		lc->y*=scale;
 		lc->z*=scale;
+		//lc->x*=sx;
+		//lc->y*=sx;
+		//lc->z*=sy;
 	}
 
 	return 0;
@@ -1301,10 +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;
 
 	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;
 	}
 
 	return moldyn->ekin;
@@ -1739,6 +2014,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	}
 #endif
 
+#ifdef PTHREADS
+	printf("##################\n");
+	printf("# USING PTHREADS #\n");
+	printf("##################\n");
+#endif
 	/* tell the world */
 	printf("[moldyn] integration start, go get a coffee ...\n");
 
@@ -1804,7 +2084,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 				dprintf(moldyn->pfd,
 				        "%f %f %f %f %f %f %f\n",moldyn->time,
 				         moldyn->p/BAR,moldyn->p_avg/BAR,
-				         moldyn->gp/BAR,moldyn->gp_avg/BAR,
+				         moldyn->p/BAR,moldyn->p_avg/BAR,
 					 moldyn->tp/BAR,moldyn->tp_avg/BAR);
 			}
 		}
@@ -1818,7 +2098,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		if(v) {
 			if(!(moldyn->total_steps%v)) {
 				dprintf(moldyn->vfd,
-				        "%f %f\n",moldyn->time,moldyn->volume);
+				        "%f %f %f %f %f\n",moldyn->time,
+				                           moldyn->volume,
+				                           moldyn->dim.x,
+				                           moldyn->dim.y,
+				                           moldyn->dim.z);
 			}
 		}
 		if(s) {
@@ -2516,6 +2800,16 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 		size-=cnt;
 	}
 
+#ifdef PTHREADS
+	amutex=malloc(moldyn->count*sizeof(pthread_mutex_t));
+	if(amutex==NULL) {
+		perror("[moldyn] load save file (mutexes)");
+		return -1;
+	}
+	for(cnt=0;cnt<moldyn->count;cnt++)
+		pthread_mutex_init(&(amutex[cnt]),NULL);
+#endif
+
 	// hooks etc ...
 
 	return 0;
@@ -3033,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;
 
@@ -3063,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