X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=moldyn.c;h=1d6e0b7dec9e3032df43f33baa692ba078bb5f7a;hp=7a9709e107a69c0f735eff1f01006611bcbb30ec;hb=HEAD;hpb=e3e1d8e2f7fc65748c8da0accb9b3e6b0b1eb049

diff --git a/moldyn.c b/moldyn.c
index 7a9709e..1d6e0b7 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -23,6 +23,10 @@
 #include <omp.h>
 #endif
 
+#if defined PTHREADS || defined VISUAL_THREAD
+#include <pthread.h>
+#endif
+
 #include "moldyn.h"
 #include "report/report.h"
 
@@ -30,6 +34,7 @@
 #include "potentials/harmonic_oscillator.h"
 #include "potentials/lennard_jones.h"
 #include "potentials/albe.h"
+#include "potentials/albe_orig.h"
 #ifdef TERSOFF_ORIG
 #include "potentials/tersoff_orig.h"
 #else
@@ -37,12 +42,25 @@
 #endif
 
 /* pse */
+#define PSE_MASS
 #define PSE_NAME
 #define PSE_COL
 #include "pse.h"
+#undef PSE_MASS
 #undef PSE_NAME
 #undef PSE_COL
 
+#ifdef PTHREADS
+/* global mutexes */
+pthread_mutex_t *amutex;
+pthread_mutex_t emutex;
+#endif
+
+/* fully constrained relaxation technique - global pointers */
+u8 crtt;
+u8 *constraints;
+double *trafo_angle;
+
 /*
  * the moldyn functions
  */
@@ -62,13 +80,28 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 	rand_init(&(moldyn->random),NULL,1);
 	moldyn->random.status|=RAND_STAT_VERBOSE;
 
+#ifdef PTHREADS
+	pthread_mutex_init(&emutex,NULL);
+#endif
+
 	return 0;
 }
 
 int moldyn_shutdown(t_moldyn *moldyn) {
 
+#ifdef PTHREADS
+	int i;
+#endif
+
 	printf("[moldyn] shutdown\n");
 
+#ifdef PTHREADS
+	for(i=0;i<moldyn->count;i++)
+		pthread_mutex_destroy(&(amutex[i]));
+	free(amutex);
+	pthread_mutex_destroy(&emutex);
+#endif
+
 	moldyn_log_shutdown(moldyn);
 	link_cell_shutdown(moldyn);
 	rand_close(&(moldyn->random));
@@ -231,26 +264,36 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 
 	switch(type) {
 		case MOLDYN_POTENTIAL_TM:
-			moldyn->func1b=tersoff_mult_1bp;
+			//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;
 			moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
 			break;
+		case MOLDYN_POTENTIAL_AO:
+			moldyn->func_j1=albe_orig_mult_3bp_j1;
+			moldyn->func_j1_k0=albe_orig_mult_3bp_k1;
+			moldyn->func_j1c=albe_orig_mult_3bp_j2;
+			moldyn->func_j1_k1=albe_orig_mult_3bp_k2;
+			moldyn->check_2b_bond=albe_orig_mult_check_2b_bond;
+			break;
 		case MOLDYN_POTENTIAL_AM:
-			moldyn->func3b_j1=albe_mult_3bp_j1;
-			moldyn->func3b_k1=albe_mult_3bp_k1;
-			moldyn->func3b_j2=albe_mult_3bp_j2;
-			moldyn->func3b_k2=albe_mult_3bp_k2;
+			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:
@@ -487,8 +530,10 @@ 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,
+                   t_offset_params *o_params) {
 
 	int new,count;
 	int ret;
@@ -496,6 +541,9 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	void *ptr;
 	t_atom *atom;
 	char name[16];
+#ifdef PTHREADS
+	pthread_mutex_t *mutex;
+#endif
 
 	new=a*b*c;
 	count=moldyn->count;
@@ -509,6 +557,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) {
@@ -518,6 +581,16 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	moldyn->atom=ptr;
 	atom=&(moldyn->atom[count]);
 
+#ifdef PTHREADS
+	ptr=realloc(amutex,(count+new)*sizeof(pthread_mutex_t));
+	if(!ptr) {
+		perror("[moldyn] mutex realloc (add atom)");
+		return -1;
+	}
+	amutex=ptr;
+	mutex=&(amutex[count]);
+#endif
+
 	/* no atoms on the boundaries (only reason: it looks better!) */
 	if(!origin) {
 		orig.x=0.5*lc;
@@ -532,22 +605,28 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 
 	switch(type) {
 		case CUBIC:
+			if(o_params->use)
+				v3_add(&orig,&orig,&(o_params->o));
 			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);
+			if(o_params->use)
+				v3_add(&orig,&orig,&(o_params->o));
 			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);
+			if(o_params->use)
+				v3_add(&orig,&orig,&(o_params->o));
 			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:
@@ -557,25 +636,62 @@ 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);
+	}
+
+	moldyn->count+=ret;
+	if(ret==new)
+		printf("[moldyn] created %s lattice with %d atoms\n",name,ret);
+
+	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[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;
 
-	moldyn->count+=new;
-	printf("[moldyn] created %s lattice with %d atoms\n",name,new);
+// WHAT ABOUT AMUTEX !!!!
 
-	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;
+#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);
@@ -583,7 +699,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;
@@ -600,6 +716,25 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 	}
 	moldyn->atom=ptr;
 
+#ifdef LOWMEM_LISTS
+	ptr=realloc(moldyn->lc.subcell->list,(count+1)*sizeof(int));
+	if(!ptr) {
+		perror("[moldyn] list realloc (add atom)");
+		return -1;
+	}
+	moldyn->lc.subcell->list=ptr;
+#endif
+
+#ifdef PTHREADS
+	ptr=realloc(amutex,(count+1)*sizeof(pthread_mutex_t));
+	if(!ptr) {
+		perror("[moldyn] mutex realloc (add atom)");
+		return -1;
+	}
+	amutex=ptr;
+	pthread_mutex_init(&(amutex[count]),NULL);
+#endif
+
 	atom=moldyn->atom;
 
 	/* initialize new atom */
@@ -607,7 +742,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;
@@ -624,6 +759,9 @@ int del_atom(t_moldyn *moldyn,int tag) {
 
 	t_atom *new,*old;
 	int cnt;
+#if defined LOWMEM_LISTS || defined PTHREADS
+	void *ptr;
+#endif
 
 	old=moldyn->atom;
 
@@ -646,16 +784,81 @@ int del_atom(t_moldyn *moldyn,int tag) {
 
 	free(old);
 
+#ifdef LOWMEM_LISTS
+	ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int));
+	if(!ptr) {
+		perror("[moldyn] list realloc (del atom)");
+		return -1;
+	}
+	moldyn->lc.subcell->list=ptr;
+	// update lists
+	link_cell_update(moldyn);
+#endif
+
+#ifdef PTHREADS
+	ptr=realloc(amutex,moldyn->count*sizeof(pthread_mutex_t));
+	if(!ptr) {
+		perror("[moldyn] mutex realloc (add atom)");
+		return -1;
+	}
+	amutex=ptr;
+	pthread_mutex_destroy(&(amutex[moldyn->count+1]));
+#endif
+
+
 	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)
@@ -663,14 +866,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;
@@ -688,12 +931,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)
@@ -710,6 +959,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++) {
@@ -718,15 +974,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;
 		}
@@ -743,12 +1065,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;
@@ -756,7 +1079,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;
 }
@@ -834,7 +1157,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;
 }
@@ -934,8 +1259,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;
 }
@@ -956,9 +1281,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;
 }
@@ -983,11 +1315,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;
@@ -1025,14 +1357,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;
 
@@ -1186,11 +1518,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);
@@ -1205,9 +1566,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) {
@@ -1229,6 +1602,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;
@@ -1242,10 +1618,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;
@@ -1293,7 +1675,9 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
 int link_cell_init(t_moldyn *moldyn,u8 vol) {
 
 	t_linkcell *lc;
+#ifndef LOWMEM_LISTS
 	int i;
+#endif
 
 	lc=&(moldyn->lc);
 
@@ -1308,8 +1692,8 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 
 #ifdef STATIC_LISTS
 	lc->subcell=malloc(lc->cells*sizeof(int*));
-#elif LOWMEM_LIST
-	lc->subcell=malloc(t_lowmem_list);
+#elif LOWMEM_LISTS
+	lc->subcell=malloc(sizeof(t_lowmem_list));
 #else
 	lc->subcell=malloc(lc->cells*sizeof(t_list));
 #endif
@@ -1327,7 +1711,7 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 #ifdef STATIC_LISTS
 		printf("[moldyn] initializing 'static' linked cells (%d)\n",
 		       lc->cells);
-#elif LOWMEM_LIST
+#elif LOWMEM_LISTS
 		printf("[moldyn] initializing 'lowmem' linked cells (%d)\n",
 		       lc->cells);
 #else
@@ -1353,7 +1737,7 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 			       i,lc->subcell[0],lc->subcell);
 		*/
 	}
-#elif LOWMEM_LIST
+#elif LOWMEM_LISTS
 	lc->subcell->head=malloc(lc->cells*sizeof(int));
 	if(lc->subcell->head==NULL) {
 		perror("[moldyn] head init (malloc)");
@@ -1378,23 +1762,25 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 int link_cell_update(t_moldyn *moldyn) {
 
 	int count,i,j,k;
-	int nx,ny;
+	int nx,nxy;
 	t_atom *atom;
 	t_linkcell *lc;
 #ifdef STATIC_LISTS
 	int p;
+#elif LOWMEM_LISTS
+	int p;
 #endif
 
 	atom=moldyn->atom;
 	lc=&(moldyn->lc);
 
 	nx=lc->nx;
-	ny=lc->ny;
+	nxy=nx*lc->ny;
 
 	for(i=0;i<lc->cells;i++)
 #ifdef STATIC_LISTS
 		memset(lc->subcell[i],-1,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
-#elif LOWMEM_LIST
+#elif LOWMEM_LISTS
 		lc->subcell->head[i]=-1;
 #else
 		list_destroy_f(&(lc->subcell[i]));
@@ -1407,7 +1793,7 @@ int link_cell_update(t_moldyn *moldyn) {
 	
 #ifdef STATIC_LISTS
 		p=0;
-		while(lc->subcell[i+j*nx+k*nx*ny][p]!=-1)
+		while(lc->subcell[i+j*nx+k*nxy][p]!=-1)
 			p++;
 
 		if(p>=MAX_ATOMS_PER_LIST) {
@@ -1415,12 +1801,13 @@ int link_cell_update(t_moldyn *moldyn) {
 			return -1;
 		}
 
-		lc->subcell[i+j*nx+k*nx*ny][p]=count;
-#elif LOWMEM_LIST
-		lc->subcell->list[count]=list->subcell->head[i+j*nx+k*nx*ny];
-		list->subcell->head=count;
+		lc->subcell[i+j*nx+k*nxy][p]=count;
+#elif LOWMEM_LISTS
+		p=i+j*nx+k*nxy;
+		lc->subcell->list[count]=lc->subcell->head[p];
+		lc->subcell->head[p]=count;
 #else
-		list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
+		list_add_immediate_f(&(lc->subcell[i+j*nx+k*nxy]),
 		                     &(atom[count]));
 		/*
 		if(j==0&&k==0)
@@ -1436,6 +1823,8 @@ int link_cell_update(t_moldyn *moldyn) {
 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
 #ifdef STATIC_LISTS
                               int **cell
+#elif LOWMEM_LISTS
+                              int *cell
 #else
                               t_list *cell
 #endif
@@ -1461,7 +1850,11 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
 		printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n",
 		       i,nx,j,ny,k,nz);
 
+#ifndef LOWMEM_LISTS
 	cell[0]=lc->subcell[i+j*nx+k*a];
+#else
+	cell[0]=lc->subcell->head[i+j*nx+k*a];
+#endif
 	for(ci=-1;ci<=1;ci++) {
 		bx=0;
 		x=i+ci;
@@ -1485,10 +1878,19 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
 				}
 				if(!(ci|cj|ck)) continue;
 				if(bx|by|bz) {
+#ifndef LOWMEM_LISTS
 					cell[--count2]=lc->subcell[x+y*nx+z*a];
+#else
+				cell[--count2]=lc->subcell->head[x+y*nx+z*a];
+#endif
+					
 				}
 				else {
+#ifndef LOWMEM_LISTS
 					cell[count1++]=lc->subcell[x+y*nx+z*a];
+#else
+				cell[count1++]=lc->subcell->head[x+y*nx+z*a];
+#endif
 				}
 			}
 		}
@@ -1501,11 +1903,19 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
 
 int link_cell_shutdown(t_moldyn *moldyn) {
 
+#ifndef LOWMEM_LISTS
 	int i;
+#endif
 	t_linkcell *lc;
 
 	lc=&(moldyn->lc);
 
+#if LOWMEM_LISTS
+	free(lc->subcell->head);
+	free(lc->subcell->list);
+
+#else
+
 	for(i=0;i<lc->cells;i++) {
 #ifdef STATIC_LISTS
 		free(lc->subcell[i]);
@@ -1514,6 +1924,7 @@ int link_cell_shutdown(t_moldyn *moldyn) {
 		list_destroy_f(&(lc->subcell[i]));
 #endif
 	}
+#endif
 
 	free(lc->subcell);
 
@@ -1582,6 +1993,17 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	struct timeval t1,t2;
 	//double tp;
 
+#ifdef VISUAL_THREAD
+	u8 first,change;
+	pthread_t io_thread;
+	int ret;
+	t_moldyn md_copy;
+	t_atom *atom_copy;
+
+	first=1;
+	change=0;
+#endif
+
 	sched=&(moldyn->schedule);
 	atom=moldyn->atom;
 
@@ -1630,6 +2052,21 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	/* debugging, ignore */
 	moldyn->debug=0;
 
+	/* zero & init moldyn copy */
+#ifdef VISUAL_THREAD
+	memset(&md_copy,0,sizeof(t_moldyn));
+	atom_copy=malloc(moldyn->count*sizeof(t_atom));
+	if(atom_copy==NULL) {
+		perror("[moldyn] malloc atom copy (init)");
+		return -1;
+	}
+#endif
+
+#ifdef PTHREADS
+	printf("##################\n");
+	printf("# USING PTHREADS #\n");
+	printf("##################\n");
+#endif
 	/* tell the world */
 	printf("[moldyn] integration start, go get a coffee ...\n");
 
@@ -1657,11 +2094,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		temperature_calc(moldyn);
 		virial_sum(moldyn);
 		pressure_calc(moldyn);
-		/*
+#ifdef PDEBUG	
 		thermodynamic_pressure_calc(moldyn);
 		printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
 		       moldyn->tp/BAR);
-		*/
+#endif
 
 		/* calculate fluctuations + averages */
 		average_and_fluctuation_calc(moldyn);
@@ -1676,10 +2113,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		if(e) {
 			if(!(moldyn->total_steps%e))
 				dprintf(moldyn->efd,
-				        "%f %f %f %f\n",
+				        "%f %f %f %f %f %f\n",
 				        moldyn->time,moldyn->ekin/energy_scale,
 				        moldyn->energy/energy_scale,
-				        get_total_energy(moldyn)/energy_scale);
+				        get_total_energy(moldyn)/energy_scale,
+				        moldyn->ekin/EV,moldyn->energy/EV);
 		}
 		if(m) {
 			if(!(moldyn->total_steps%m)) {
@@ -1695,7 +2133,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);
 			}
 		}
@@ -1709,12 +2147,16 @@ 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) {
 			if(!(moldyn->total_steps%s)) {
-				snprintf(dir,128,"%s/s-%07.f.save",
+				snprintf(dir,128,"%s/s-%08.f.save",
 				         moldyn->vlsdir,moldyn->time);
 				fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
 				        S_IRUSR|S_IWUSR);
@@ -1729,20 +2171,54 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		}
 		if(a) {
 			if(!(moldyn->total_steps%a)) {
+#ifdef VISUAL_THREAD
+	/* check whether thread has not terminated yet */
+	if(!first) {
+		ret=pthread_join(io_thread,NULL);
+	}
+	first=0;
+	/* prepare and start thread */
+	if(moldyn->count!=md_copy.count) {
+		free(atom_copy);
+		change=1;
+	}
+	memcpy(&md_copy,moldyn,sizeof(t_moldyn));
+	if(change) {
+		atom_copy=malloc(moldyn->count*sizeof(t_atom));
+		if(atom_copy==NULL) {
+			perror("[moldyn] malloc atom copy (change)");
+			return -1;
+		}
+	}
+	md_copy.atom=atom_copy;
+	memcpy(atom_copy,moldyn->atom,moldyn->count*sizeof(t_atom));
+	change=0;
+	ret=pthread_create(&io_thread,NULL,visual_atoms,&md_copy);
+	if(ret) {
+		perror("[moldyn] create visual atoms thread\n");
+		return -1;
+	}
+#else
 				visual_atoms(moldyn);
+#endif
 			}
 		}
 
 		/* display progress */
-		//if(!(moldyn->total_steps%10)) {
+#ifndef PDEBUG
+		if(!(i%10)) {
+#endif
 			/* get current time */
 			gettimeofday(&t2,NULL);
 
-printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
+printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
        sched->count,i,moldyn->total_steps,
        moldyn->t,moldyn->t_avg,
+#ifndef PDEBUG
        moldyn->p/BAR,moldyn->p_avg/BAR,
-       //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#else
+       moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#endif
        moldyn->volume,
        (int)(t2.tv_sec-t1.tv_sec));
 
@@ -1750,7 +2226,9 @@ printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
 
 			/* copy over time */
 			t1=t2;
-		//}
+#ifndef PDEBUG
+		}
+#endif
 
 		/* increase absolute time */
 		moldyn->time+=moldyn->tau;
@@ -1771,6 +2249,59 @@ printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
 
 	}
 
+	/* writing a final save file! */
+	if(s) {
+		snprintf(dir,128,"%s/s-final.save",moldyn->vlsdir);
+		fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,S_IRUSR|S_IWUSR);
+		if(fd<0) perror("[moldyn] save fd open");
+		else {
+			write(fd,moldyn,sizeof(t_moldyn));
+			write(fd,moldyn->atom,
+			      moldyn->count*sizeof(t_atom));
+		}
+		close(fd);
+	}
+	/* writing a final visual file! */
+	if(a)
+		visual_atoms(moldyn);
+
+	return 0;
+}
+
+/* basis trafo */
+
+#define FORWARD		0
+#define BACKWARD	1
+
+int basis_trafo(t_3dvec *r,u8 dir,double z,double x) {
+
+	t_3dvec tmp;
+
+	if(dir==FORWARD) {
+		if(z!=0.0) {
+			v3_copy(&tmp,r);
+			r->x=cos(z)*tmp.x-sin(z)*tmp.y;
+			r->y=sin(z)*tmp.x+cos(z)*tmp.y;
+		}
+		if(x!=0.0) {
+			v3_copy(&tmp,r);
+			r->y=cos(x)*tmp.y-sin(x)*tmp.z;
+			r->z=sin(x)*tmp.y+cos(x)*tmp.z;
+		}
+	}
+	else {
+		if(x!=0.0) {
+			v3_copy(&tmp,r);
+			r->y=cos(-x)*tmp.y-sin(-x)*tmp.z;
+			r->z=sin(-x)*tmp.y+cos(-x)*tmp.z;
+		}
+		if(z!=0.0) {
+			v3_copy(&tmp,r);
+			r->x=cos(-z)*tmp.x-sin(-z)*tmp.y;
+			r->y=sin(-z)*tmp.x+cos(-z)*tmp.y;
+		}
+	}
+
 	return 0;
 }
 
@@ -1789,15 +2320,47 @@ int velocity_verlet(t_moldyn *moldyn) {
 	tau_square=moldyn->tau_square;
 
 	for(i=0;i<count;i++) {
+
 		/* check whether fixed atom */
 		if(atom[i].attr&ATOM_ATTR_FP)
 			continue;
+
 		/* new positions */
 		h=0.5/atom[i].mass;
+
+		/* constraint relaxation */
+		if(crtt) {
+			// forces
+			basis_trafo(&(atom[i].f),FORWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			if(constraints[3*i])
+				atom[i].f.x=0;
+			if(constraints[3*i+1])
+				atom[i].f.y=0;
+			if(constraints[3*i+2])
+				atom[i].f.z=0;
+			basis_trafo(&(atom[i].f),BACKWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			// velocities
+			basis_trafo(&(atom[i].v),FORWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			if(constraints[3*i])
+				atom[i].v.x=0;
+			if(constraints[3*i+1])
+				atom[i].v.y=0;
+			if(constraints[3*i+2])
+				atom[i].v.z=0;
+			basis_trafo(&(atom[i].v),BACKWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+		}
+
+#ifndef QUENCH
 		v3_scale(&delta,&(atom[i].v),tau);
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
+#endif
 		v3_scale(&delta,&(atom[i].f),h*tau_square);
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
+		//check_per_bound_and_care_for_pbc(moldyn,&(atom[i]));
 		check_per_bound(moldyn,&(atom[i].r));
 
 		/* velocities [actually v(t+tau/2)] */
@@ -1813,7 +2376,11 @@ int velocity_verlet(t_moldyn *moldyn) {
 
 	/* forces depending on chosen potential */
 #ifndef ALBE_FAST
-	potential_force_calc(moldyn);
+	// if albe, use albe force calc routine
+	//if(moldyn->func3b_j1==albe_mult_3bp_j1)
+	//	albe_potential_force_calc(moldyn);
+	//else
+		potential_force_calc(moldyn);
 #else
 	albe_potential_force_calc(moldyn);
 #endif
@@ -1822,6 +2389,33 @@ int velocity_verlet(t_moldyn *moldyn) {
 		/* check whether fixed atom */
 		if(atom[i].attr&ATOM_ATTR_FP)
 			continue;
+
+		/* constraint relaxation */
+		if(crtt) {
+			// forces
+			basis_trafo(&(atom[i].f),FORWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			if(constraints[3*i])
+				atom[i].f.x=0;
+			if(constraints[3*i+1])
+				atom[i].f.y=0;
+			if(constraints[3*i+2])
+				atom[i].f.z=0;
+			basis_trafo(&(atom[i].f),BACKWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			// velocities
+			basis_trafo(&(atom[i].v),FORWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+			if(constraints[3*i])
+				atom[i].v.x=0;
+			if(constraints[3*i+1])
+				atom[i].v.y=0;
+			if(constraints[3*i+2])
+				atom[i].v.z=0;
+			basis_trafo(&(atom[i].v),BACKWARD,
+			            trafo_angle[2*i],trafo_angle[2*i+1]);
+		}
+
 		/* again velocities [actually v(t+tau)] */
 		v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 		v3_add(&(atom[i].v),&(atom[i].v),&delta);
@@ -1849,6 +2443,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 	int *neighbour_i[27];
 	int p,q;
 	t_atom *atom;
+#elif LOWMEM_LISTS
+	int neighbour_i[27];
+	int p,q;
 #else
 	t_list neighbour_i[27];
 	t_list neighbour_i2[27];
@@ -1901,8 +2498,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;
@@ -1917,8 +2514,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;
@@ -1929,6 +2532,13 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 					jtom=&(atom[neighbour_i[j][p]]);
 					p++;
+#elif LOWMEM_LISTS
+				p=neighbour_i[j];
+
+				while(p!=-1) {
+
+					jtom=&(itom[p]);
+					p=lc->subcell->list[p];
 #else
 				this=&(neighbour_i[j]);
 				list_reset_f(this);
@@ -1943,15 +2553,74 @@ int potential_force_calc(t_moldyn *moldyn) {
 					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
 					}
 #ifdef STATIC_LISTS
 				}
+#elif LOWMEM_LISTS
+				}
 #else
 				} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
@@ -1959,7 +2628,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 			}
 		}
 
-		/* 3 body potential/force */
+		/* continued 3 body potential/force */
 
 		if(!(itom[i].attr&ATOM_ATTR_3BP))
 			continue;
@@ -1967,6 +2636,8 @@ int potential_force_calc(t_moldyn *moldyn) {
 		/* copy the neighbour lists */
 #ifdef STATIC_LISTS
 		/* no copy needed for static lists */
+#elif LOWMEM_LISTS
+		/* no copy needed for lowmem lists */
 #else
 		memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
 #endif
@@ -1982,6 +2653,13 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 				jtom=&(atom[neighbour_i[j][p]]);
 				p++;
+#elif LOWMEM_LISTS
+				p=neighbour_i[j];
+
+				while(p!=-1) {
+
+					jtom=&(itom[p]);
+					p=lc->subcell->list[p];
 #else
 			this=&(neighbour_i[j]);
 			list_reset_f(this);
@@ -2003,18 +2681,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++) {
 
@@ -2026,6 +2704,13 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 						ktom=&(atom[neighbour_i[k][q]]);
 						q++;
+#elif LOWMEM_LISTS
+					q=neighbour_i[k];
+
+					while(q!=-1) {
+
+						ktom=&(itom[q]);
+						q=lc->subcell->list[q];
 #else
 					that=&(neighbour_i2[k]);
 					list_reset_f(that);
@@ -2040,8 +2725,8 @@ 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;
@@ -2051,8 +2736,11 @@ int potential_force_calc(t_moldyn *moldyn) {
 						                  jtom,
 						                  ktom,
 						                  bc_ik|bc_ij);
+
 #ifdef STATIC_LISTS
 					}
+#elif LOWMEM_LISTS
+					}
 #else
 					} while(list_next_f(that)!=\
 					        L_NO_NEXT_ELEMENT);
@@ -2062,14 +2750,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++) {
 
@@ -2081,6 +2770,13 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 						ktom=&(atom[neighbour_i[k][q]]);
 						q++;
+#elif LOWMEM_LISTS
+					q=neighbour_i[k];
+
+					while(q!=-1) {
+
+						ktom=&(itom[q]);
+						q=lc->subcell->list[q];
 #else
 					that=&(neighbour_i2[k]);
 					list_reset_f(that);
@@ -2095,20 +2791,22 @@ 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
 					}
+#elif LOWMEM_LISTS
+					}
 #else
 					} while(list_next_f(that)!=\
 					        L_NO_NEXT_ELEMENT);
@@ -2118,20 +2816,22 @@ 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
 			}
+#elif LOWMEM_LISTS
+			}
 #else
 			} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
 		
 		}
-		
+
 #ifdef DEBUG
 	//printf("\n\n");
 #endif
@@ -2223,6 +2923,51 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
 	return 0;
 }
         
+int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a) {
+	
+	double x,y,z;
+	t_3dvec *dim;
+
+	dim=&(moldyn->dim);
+
+	x=dim->x/2;
+	y=dim->y/2;
+	z=dim->z/2;
+
+	if(moldyn->status&MOLDYN_STAT_PBX) {
+		if(a->r.x>=x) {
+			a->pbc[0]+=1;
+			a->r.x-=dim->x;
+		}
+		else if(-a->r.x>x) {
+			a->pbc[0]-=1;
+			a->r.x+=dim->x;
+		}
+	}
+	if(moldyn->status&MOLDYN_STAT_PBY) {
+		if(a->r.y>=y) {
+			a->pbc[1]+=1;
+			a->r.y-=dim->y;
+		}
+		else if(-a->r.y>y) {
+			a->pbc[1]-=1;
+			a->r.y+=dim->y;
+		}
+	}
+	if(moldyn->status&MOLDYN_STAT_PBZ) {
+		if(a->r.z>=z) {
+			a->pbc[2]+=1;
+			a->r.z-=dim->z;
+		}
+		else if(-a->r.z>z) {
+			a->pbc[2]-=1;
+			a->r.z+=dim->z;
+		}
+	}
+
+	return 0;
+}
+        
 /*
  * debugging / critical check functions
  */
@@ -2313,7 +3058,7 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 	if(fsize!=sizeof(t_moldyn)+size) {
 		corr=fsize-sizeof(t_moldyn)-size;
 		printf("[moldyn] WARNING: lsf (illegal file size)\n");
-		printf("  moifying offset:\n");
+		printf("  modifying offset:\n");
 		printf("  - current pos: %d\n",sizeof(t_moldyn));
 		printf("  - atom size: %d\n",size);
 		printf("  - file size: %d\n",fsize);
@@ -2336,6 +3081,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;
@@ -2367,6 +3122,9 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
 #ifdef STATIC_LISTS
 	int *neighbour[27];
 	int p;
+#elif LOWMEM_LISTS
+	int neighbour[27];
+	int p;
 #else
 	t_list neighbour[27];
 	t_list *this;
@@ -2397,6 +3155,13 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
 
 				jtom=&(moldyn->atom[neighbour[j][p]]);
 				p++;
+#elif LOWMEM_LISTS
+			p=neighbour[j];
+
+			while(p!=-1) {
+
+				jtom=&(itom[p]);
+				p=lc->subcell->list[p];
 #else
 			this=&(neighbour[j]);
 			list_reset_f(this);
@@ -2414,6 +3179,8 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
 
 #ifdef STATIC_LISTS
 			}
+#elif LOWMEM_LISTS
+			}
 #else
 			} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
@@ -2424,6 +3191,84 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
 
 }
 
+/*
+ * function to find neighboured atoms
+ */
+
+int process_neighbours(t_moldyn *moldyn,void *data,t_atom *atom,
+                       int (*process)(t_moldyn *moldyn,t_atom *atom,t_atom *natom,
+		                      void *data,u8 bc)) {
+
+	t_linkcell *lc;
+#ifdef STATIC_LISTS
+	int *neighbour[27];
+	int p;
+#elif LOWMEM_LISTS
+	int neighbour[27];
+	int p;
+#else
+	t_list neighbour[27];
+	t_list *this;
+#endif
+	u8 bc;
+	t_atom *natom;
+	int j;
+
+	lc=&(moldyn->lc);
+	
+	/* neighbour indexing */
+	link_cell_neighbour_index(moldyn,
+	                          (atom->r.x+moldyn->dim.x/2)/lc->x,
+	                          (atom->r.y+moldyn->dim.y/2)/lc->x,
+	                          (atom->r.z+moldyn->dim.z/2)/lc->x,
+	                          neighbour);
+
+	for(j=0;j<27;j++) {
+
+		bc=(j<lc->dnlc)?0:1;
+
+#ifdef STATIC_LISTS
+		p=0;
+
+		while(neighbour[j][p]!=-1) {
+
+			natom=&(moldyn->atom[neighbour[j][p]]);
+			p++;
+#elif LOWMEM_LISTS
+		p=neighbour[j];
+
+		while(p!=-1) {
+
+			natom=&(moldyn->atom[p]);
+			p=lc->subcell->list[p];
+#else
+		this=&(neighbour[j]);
+		list_reset_f(this);
+
+		if(this->start==NULL)
+			continue;
+
+		do {
+
+			natom=this->current->data;
+#endif
+
+			/* process bond */
+			process(moldyn,atom,natom,data,bc);
+
+#ifdef STATIC_LISTS
+		}
+#elif LOWMEM_LISTS
+		}
+#else
+		} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
+	}
+
+	return 0;
+
+}
+
 /*
  * post processing functions
  */
@@ -2458,6 +3303,7 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 	int i;
 	t_atom *atom;
 	t_3dvec dist;
+	t_3dvec final_r;
 	double d2;
 	int a_cnt;
 	int b_cnt;
@@ -2471,8 +3317,12 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 
 	for(i=0;i<moldyn->count;i++) {
 
-		v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
-		check_per_bound(moldyn,&dist);
+		/* care for pb crossing */
+		final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x;
+		final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y;
+		final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z;
+		/* calculate distance */
+		v3_sub(&dist,&final_r,&(atom[i].r_0));
 		d2=v3_absolute_square(&dist);
 
 		if(atom[i].brand) {
@@ -2494,6 +3344,57 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 	return 0;
 }
 
+int calculate_msd(t_moldyn *moldyn,double *msd) {
+
+	int i;
+	t_atom *atom;
+	t_3dvec dist;
+	t_3dvec final_r;
+	double d2;
+	int a_cnt;
+	int b_cnt;
+
+	atom=moldyn->atom;
+	msd[0]=0;
+	msd[1]=0;
+	msd[2]=0;
+	a_cnt=0;
+	b_cnt=0;
+
+	for(i=0;i<moldyn->count;i++) {
+
+		/* care for pb crossing */
+		if(atom[i].pbc[0]|atom[i].pbc[1]|atom[i].pbc[2]) {
+			printf("[moldyn] msd pb crossings for atom %d\n",i);
+			printf("  x: %d y: %d z: %d\n",
+			       atom[i].pbc[0],atom[i].pbc[1],atom[i].pbc[2]);
+		}
+		final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x;
+		final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y;
+		final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z;
+		/* calculate distance */
+		v3_sub(&dist,&final_r,&(atom[i].r_0));
+		d2=v3_absolute_square(&dist);
+
+		if(atom[i].brand) {
+			b_cnt+=1;
+			msd[1]+=d2;
+		}
+		else {
+			a_cnt+=1;
+			msd[0]+=d2;
+		}
+
+		msd[2]+=d2;
+	}
+
+	msd[0]/=a_cnt;
+	msd[1]/=b_cnt;
+	msd[2]/=moldyn->count;
+		
+	return 0;
+}
+
 int bonding_analyze(t_moldyn *moldyn,double *cnt) {
 
 	return 0;
@@ -2639,9 +3540,8 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 	ba=data;
 
 	/* increase total bond counter
-	 * ... double counting!
 	 */
-	ba->tcnt+=2;
+	ba->tcnt+=1;
 
 	if(itom->brand==0)
 		ba->acnt[jtom->tag]+=1;
@@ -2658,10 +3558,11 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 
 int bond_analyze(t_moldyn *moldyn,double *quality) {
 
-	// by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
-
-	int qcnt;
-	int ccnt,cset;
+	int qcnt4;
+	int qcnt3;
+	int ccnt4;
+	int ccnt3;
+	int bcnt;
 	t_ba ba;
 	int i;
 	t_atom *atom;
@@ -2681,9 +3582,9 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
 	memset(ba.bcnt,0,moldyn->count*sizeof(int));
 
 	ba.tcnt=0;
-	qcnt=0;
-	ccnt=0;
-	cset=0;
+	qcnt4=0; qcnt3=0;
+	ccnt4=0; ccnt3=0;
+	bcnt=0;
 
 	atom=moldyn->atom;
 
@@ -2692,27 +3593,30 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
 	for(i=0;i<moldyn->count;i++) {
 		if(atom[i].brand==0) {
 			if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
-				qcnt+=4;
+				qcnt4+=4;
+			if((ba.acnt[i]==0)&(ba.bcnt[i]==3))
+				qcnt3+=3;
 		}
 		else {
 			if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
-				qcnt+=4;
-				ccnt+=1;
+				qcnt4+=4;
+				ccnt4+=1;
+			}
+			if((ba.acnt[i]==3)&(ba.bcnt[i]==0)) {
+				qcnt3+=4;
+				ccnt3+=1;
 			}
-			cset+=1;
+			bcnt+=1;
 		}
 	}
 
-	printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
-	printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
-
 	if(quality) {
-		quality[0]=1.0*ccnt/cset;
-		quality[1]=1.0*qcnt/ba.tcnt;
+		quality[0]=1.0*ccnt4/bcnt;
+		quality[1]=1.0*ccnt3/bcnt;
 	}
 	else {
-		printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
-		printf("[moldyn] bond analyze:   tot_q=%f\n",1.0*qcnt/ba.tcnt);
+		printf("[moldyn] bond analyze: %f %f\n",
+		       1.0*ccnt4/bcnt,1.0*ccnt3/bcnt);
 	}
 
 	return 0;
@@ -2750,7 +3654,11 @@ int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
 	return 0;
 }
 
+#ifdef VISUAL_THREAD
+void *visual_atoms(void *ptr) {
+#else
 int visual_atoms(t_moldyn *moldyn) {
+#endif
 
 	int i;
  	char file[128+64];
@@ -2759,6 +3667,12 @@ int visual_atoms(t_moldyn *moldyn) {
 	t_visual *v;
 	t_atom *atom;
 	t_vb vb;
+	t_3dvec strain;
+#ifdef VISUAL_THREAD
+	t_moldyn *moldyn;
+
+	moldyn=ptr;
+#endif
 
 	v=&(moldyn->vis);
 	dim.x=v->dim.x;
@@ -2768,31 +3682,39 @@ int visual_atoms(t_moldyn *moldyn) {
 
 	help=(dim.x+dim.y);
 
-	sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
+	sprintf(file,"%s/atomic_conf_%08.f.xyz",v->fb,moldyn->time);
 	vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
 	if(vb.fd<0) {
 		perror("open visual save file fd");
+#ifndef VISUAL_THREAD
 		return -1;
+#endif
 	}
 
 	/* write the actual data file */
 
 	// povray header
-	dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n",
+	dprintf(vb.fd,"# [P] %d %08.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++) {
+		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
 	process_2b_bonds(moldyn,&vb,visual_bonds_process);
+#endif
 	
 	// boundaries
 	if(dim.x) {
@@ -2838,8 +3760,15 @@ int visual_atoms(t_moldyn *moldyn) {
 
 	close(vb.fd);
 
+#ifdef VISUAL_THREAD
+	pthread_exit(NULL);
+
+}
+#else
+
 	return 0;
 }
+#endif
 
 /*
  * fpu cntrol functions