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

diff --git a/moldyn.c b/moldyn.c
index ad4f900..1aff5a8 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -20,53 +20,15 @@
 #include "moldyn.h"
 #include "report/report.h"
 
-/*
- * global variables, pse and atom colors (only needed here)
- */ 
-
-static char *pse_name[]={
-	"*",
-	"H",
-	"He",
-	"Li",
-	"Be",
-	"B",
-	"C",
-	"N",
-	"O",
-	"F",
-	"Ne",
-	"Na",
-	"Mg",
-	"Al",
-	"Si",
-	"P",
-	"S",
-	"Cl",
-	"Ar",
-};
-
-static char *pse_col[]={
-	"*",
-	"White",
-	"He",
-	"Li",
-	"Be",
-	"B",
-	"Gray",
-	"N",
-	"Blue",
-	"F",
-	"Ne",
-	"Na",
-	"Mg",
-	"Al",
-	"Yellow",
-	"P",
-	"S",
-	"Cl",
-	"Ar",
-};
+/* potential includes */
+#include "potentials/harmonic_oscillator.h"
+#include "potentials/lennard_jones.h"
+#include "potentials/albe.h"
+#ifdef TERSOFF_ORIG
+#include "potentials/tersoff_orig.h"
+#else
+#include "potentials/tersoff.h"
+#endif
 
 /*
  * the moldyn functions
@@ -120,24 +82,13 @@ int set_int_alg(t_moldyn *moldyn,u8 algo) {
 int set_cutoff(t_moldyn *moldyn,double cutoff) {
 
 	moldyn->cutoff=cutoff;
+	moldyn->cutoff_square=cutoff*cutoff;
 
 	printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff);
 
 	return 0;
 }
 
-int set_bondlen(t_moldyn *moldyn,double b0,double b1,double bm) {
-
-	moldyn->bondlen[0]=b0*b0;
-	moldyn->bondlen[1]=b1*b1;
-	if(bm<0)
-		moldyn->bondlen[2]=b0*b1;
-	else
-		moldyn->bondlen[2]=bm*bm;
-
-	return 0;
-}
-
 int set_temperature(t_moldyn *moldyn,double t_ref) {
 
 	moldyn->t_ref=t_ref;
@@ -156,6 +107,38 @@ int set_pressure(t_moldyn *moldyn,double p_ref) {
 	return 0;
 }
 
+int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) {
+
+	moldyn->pt_scale&=(~(P_SCALE_MASK));
+	moldyn->pt_scale|=ptype;
+	moldyn->p_tc=ptc;
+
+	printf("[moldyn] p scaling:\n");
+
+	printf("  p: %s",ptype?"yes":"no ");
+	if(ptype)
+		printf(" | type: %02x | factor: %f",ptype,ptc);
+	printf("\n");
+
+	return 0;
+}
+
+int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) {
+
+	moldyn->pt_scale&=(~(T_SCALE_MASK));
+	moldyn->pt_scale|=ttype;
+	moldyn->t_tc=ttc;
+
+	printf("[moldyn] t scaling:\n");
+
+	printf("  t: %s",ttype?"yes":"no ");
+	if(ttype)
+		printf(" | type: %02x | factor: %f",ttype,ttc);
+	printf("\n");
+
+	return 0;
+}
+
 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
 
 	moldyn->pt_scale=(ptype|ttype);
@@ -191,15 +174,12 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
 		moldyn->vis.dim.z=z;
 	}
 
-	moldyn->dv=0.000001*moldyn->volume;
-
 	printf("[moldyn] dimensions in A and A^3 respectively:\n");
 	printf("  x: %f\n",moldyn->dim.x);
 	printf("  y: %f\n",moldyn->dim.y);
 	printf("  z: %f\n",moldyn->dim.z);
 	printf("  volume: %f\n",moldyn->volume);
 	printf("  visualize simulation box: %s\n",visualize?"yes":"no");
-	printf("  delta volume (pressure calc): %f\n",moldyn->dv);
 
 	return 0;
 }
@@ -231,58 +211,49 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
 	return 0;
 }
 
-int set_potential1b(t_moldyn *moldyn,pf_func1b func) {
-
-	moldyn->func1b=func;
-
-	return 0;
-}
-
-int set_potential2b(t_moldyn *moldyn,pf_func2b func) {
-
-	moldyn->func2b=func;
-
-	return 0;
-}
-
-int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) {
-
-	moldyn->func3b_j1=func;
-
-	return 0;
-}
-
-int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) {
-
-	moldyn->func3b_j2=func;
-
-	return 0;
-}
-
-int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) {
-
-	moldyn->func3b_j3=func;
-
-	return 0;
-}
-
-int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) {
+int set_potential(t_moldyn *moldyn,u8 type) {
 
-	moldyn->func3b_k1=func;
-
-	return 0;
-}
-
-int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) {
-
-	moldyn->func3b_k2=func;
-
-	return 0;
-}
-
-int set_potential_params(t_moldyn *moldyn,void *params) {
-
-	moldyn->pot_params=params;
+	switch(type) {
+		case MOLDYN_POTENTIAL_TM:
+			moldyn->func_i0=tersoff_mult_1bp;
+			moldyn->func_j1=tersoff_mult_3bp_j1;
+			moldyn->func_j1_k0=tersoff_mult_3bp_k1;
+			moldyn->func_j1c=tersoff_mult_3bp_j2;
+			moldyn->func_j1_k1=tersoff_mult_3bp_k2;
+			moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
+			break;
+		/*
+		case MOLDYN_POTENTIAL_AM:
+			moldyn->func3b_j1=albe_mult_3bp_j1;
+			moldyn->func3b_k1=albe_mult_3bp_k1;
+			moldyn->func3b_j2=albe_mult_3bp_j2;
+			moldyn->func3b_k2=albe_mult_3bp_k2;
+			moldyn->check_2b_bond=albe_mult_check_2b_bond;
+			break;
+		*/
+		case MOLDYN_POTENTIAL_AM:
+			moldyn->func_i0=albe_mult_i0;
+			moldyn->func_j0=albe_mult_i0_j0;
+			moldyn->func_j0_k0=albe_mult_i0_j0_k0;
+			moldyn->func_j0e=albe_mult_i0_j1;
+			moldyn->func_j1=albe_mult_i0_j2;
+			moldyn->func_j1_k0=albe_mult_i0_j2_k0;
+			moldyn->func_j1c=albe_mult_i0_j3;
+			moldyn->check_2b_bond=albe_mult_check_2b_bond;
+			break;
+		case MOLDYN_POTENTIAL_HO:
+			moldyn->func_j0=harmonic_oscillator;
+			moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
+			break;
+		case MOLDYN_POTENTIAL_LJ:
+			moldyn->func_j0=lennard_jones;
+			moldyn->check_2b_bond=lennard_jones_check_2b_bond;
+			break;
+		default:
+			printf("[moldyn] set potential: unknown type %02x\n",
+			       type);
+			return -1;
+	}
 
 	return 0;
 }
@@ -370,12 +341,25 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 			dprintf(moldyn->tfd,"# temperature log file\n");
 			printf("temperature (%d)\n",timer);
 			break;
+		case LOG_VOLUME:
+			moldyn->vwrite=timer;
+			snprintf(filename,127,"%s/volume",moldyn->vlsdir);
+			moldyn->vfd=open(filename,
+			                 O_WRONLY|O_CREAT|O_EXCL,
+			                 S_IRUSR|S_IWUSR);
+			if(moldyn->vfd<0) {
+				perror("[moldyn] volume log file\n");
+				return moldyn->vfd;
+			}
+			dprintf(moldyn->vfd,"# volume log file\n");
+			printf("volume (%d)\n",timer);
+			break;
 		case SAVE_STEP:
 			moldyn->swrite=timer;
 			printf("save file (%d)\n",timer);
 			break;
 		case VISUAL_STEP:
-			moldyn->vwrite=timer;
+			moldyn->awrite=timer;
 			ret=visual_init(moldyn,moldyn->vlsdir);
 			if(ret<0) {
 				printf("[moldyn] visual init failure\n");
@@ -507,6 +491,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	t_3dvec orig;
 	void *ptr;
 	t_atom *atom;
+	char name[16];
 
 	new=a*b*c;
 	count=moldyn->count;
@@ -541,18 +526,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 		case CUBIC:
 			set_nn_dist(moldyn,lc);
 			ret=cubic_init(a,b,c,lc,atom,&orig);
+			strcpy(name,"cubic");
 			break;
 		case FCC:
 			if(!origin)
 				v3_scale(&orig,&orig,0.5);
 			set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
 			ret=fcc_init(a,b,c,lc,atom,&orig);
+			strcpy(name,"fcc");
 			break;
 		case DIAMOND:
 			if(!origin)
 				v3_scale(&orig,&orig,0.25);
 			set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
 			ret=diamond_init(a,b,c,lc,atom,&orig);
+			strcpy(name,"diamond");
 			break;
 		default:
 			printf("unknown lattice type (%02x)\n",type);
@@ -569,7 +557,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	}
 
 	moldyn->count+=new;
-	printf("[moldyn] created lattice with %d atoms\n",new);
+	printf("[moldyn] created %s lattice with %d atoms\n",name,new);
 
 	for(ret=0;ret<new;ret++) {
 		atom[ret].element=element;
@@ -595,7 +583,7 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 	int count;
 	
 	atom=moldyn->atom;
-	count=(moldyn->count)++;
+	count=(moldyn->count)++;	// asshole style!
 
 	ptr=realloc(atom,(count+1)*sizeof(t_atom));
 	if(!ptr) {
@@ -605,6 +593,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 	moldyn->atom=ptr;
 
 	atom=moldyn->atom;
+
+	/* initialize new atom */
+	memset(&(atom[count]),0,sizeof(t_atom));
 	atom[count].r=*r;
 	atom[count].v=*v;
 	atom[count].element=element;
@@ -913,16 +904,26 @@ double ideal_gas_law_pressure(t_moldyn *moldyn) {
 double virial_sum(t_moldyn *moldyn) {
 
 	int i;
-	double v;
 	t_virial *virial;
 
 	/* virial (sum over atom virials) */
-	v=0.0;
+	moldyn->virial=0.0;
+	moldyn->vir.xx=0.0;
+	moldyn->vir.yy=0.0;
+	moldyn->vir.zz=0.0;
+	moldyn->vir.xy=0.0;
+	moldyn->vir.xz=0.0;
+	moldyn->vir.yz=0.0;
 	for(i=0;i<moldyn->count;i++) {
 		virial=&(moldyn->atom[i].virial);
-		v+=(virial->xx+virial->yy+virial->zz);
+		moldyn->virial+=(virial->xx+virial->yy+virial->zz);
+		moldyn->vir.xx+=virial->xx;
+		moldyn->vir.yy+=virial->yy;
+		moldyn->vir.zz+=virial->zz;
+		moldyn->vir.xy+=virial->xy;
+		moldyn->vir.xz+=virial->xz;
+		moldyn->vir.yz+=virial->yz;
 	}
-	moldyn->virial=v;
 
 	/* global virial (absolute coordinates) */
 	virial=&(moldyn->gvir);
@@ -954,12 +955,44 @@ double pressure_calc(t_moldyn *moldyn) {
 	return moldyn->p;
 }
 
+int average_reset(t_moldyn *moldyn) {
+
+	printf("[moldyn] average reset\n");
+
+	/* update skip value */
+	moldyn->avg_skip=moldyn->total_steps;
+
+	/* kinetic energy */
+	moldyn->k_sum=0.0;
+	moldyn->k2_sum=0.0;
+	
+	/* potential energy */
+	moldyn->v_sum=0.0;
+	moldyn->v2_sum=0.0;
+
+	/* temperature */
+	moldyn->t_sum=0.0;
+
+	/* virial */
+	moldyn->virial_sum=0.0;
+	moldyn->gv_sum=0.0;
+
+	/* pressure */
+	moldyn->p_sum=0.0;
+	moldyn->gp_sum=0.0;
+	moldyn->tp_sum=0.0;
+
+	return 0;
+}
+
 int average_and_fluctuation_calc(t_moldyn *moldyn) {
 
+	int denom;
+
 	if(moldyn->total_steps<moldyn->avg_skip)
 		return 0;
 
-	int denom=moldyn->total_steps+1-moldyn->avg_skip;
+	denom=moldyn->total_steps+1-moldyn->avg_skip;
 
 	/* assume up to date energies, temperature, pressure etc */
 
@@ -992,6 +1025,8 @@ int average_and_fluctuation_calc(t_moldyn *moldyn) {
 	moldyn->p_avg=moldyn->p_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;
 
 	return 0;
 }
@@ -1031,9 +1066,11 @@ printf("  --> <dV2> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->coun
 
 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 
-	t_3dvec dim,*tp;
-	double u_up,u_down,dv;
-	double scale,p;
+	t_3dvec dim;
+	//t_3dvec *tp;
+	double h,dv;
+	double y0,y1;
+	double su,sd;
 	t_atom *store;
 
 	/*
@@ -1043,54 +1080,56 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	 *
 	 */
 
-	scale=0.00001;
-	dv=8*scale*scale*scale*moldyn->volume;
-
+	/* store atomic configuration + dimension */
 	store=malloc(moldyn->count*sizeof(t_atom));
 	if(store==NULL) {
 		printf("[moldyn] allocating store mem failed\n");
 		return -1;
 	}
-
-	/* save unscaled potential energy + atom/dim configuration */
 	memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
 	dim=moldyn->dim;
 
+	/* x1, y1 */
+	sd=0.00001;
+	h=(1.0-sd)*(1.0-sd)*(1.0-sd);
+	su=pow(2.0-h,ONE_THIRD)-1.0;
+	dv=(1.0-h)*moldyn->volume;
+
 	/* scale up dimension and atom positions */
-	scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
-	scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
+	scale_dim(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
+	scale_atoms(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
 	potential_force_calc(moldyn);
-	u_up=moldyn->energy;
+	y1=moldyn->energy;
 
 	/* restore atomic configuration + dim */
 	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 	moldyn->dim=dim;
 
 	/* scale down dimension and atom positions */
-	scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
-	scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
+	scale_dim(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
+	scale_atoms(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
 	potential_force_calc(moldyn);
-	u_down=moldyn->energy;
+	y0=moldyn->energy;
 	
 	/* calculate pressure */
-	p=-(u_up-u_down)/dv;
-printf("-------> %.10f %.10f %f\n",u_up/EV/moldyn->count,u_down/EV/moldyn->count,p/BAR);
+	moldyn->tp=-(y1-y0)/(2.0*dv);
 
-	/* restore atomic configuration + dim */
+	/* restore atomic configuration */
 	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 	moldyn->dim=dim;
-
-	/* restore energy */
-	potential_force_calc(moldyn);
-
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
+	//potential_force_calc(moldyn);
 
-	return p;
+	/* free store buffer */
+	if(store)
+		free(store);
+
+	return moldyn->tp;
 }
 
 double get_pressure(t_moldyn *moldyn) {
@@ -1151,13 +1190,12 @@ int scale_volume(t_moldyn *moldyn) {
 
 	/* scaling factor */
 	if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
-		scale=1.0-(moldyn->p_ref-moldyn->p)/moldyn->p_tc;
+		scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc;
 		scale=pow(scale,ONE_THIRD);
 	}
 	else {
 		scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
 	}
-moldyn->debug=scale;
 
 	/* scale the atoms and dimensions */
 	scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
@@ -1258,25 +1296,57 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 	lc->y=moldyn->dim.y/lc->ny;
 	lc->nz=moldyn->dim.z/moldyn->cutoff;
 	lc->z=moldyn->dim.z/lc->nz;
-
 	lc->cells=lc->nx*lc->ny*lc->nz;
+
+#ifdef STATIC_LISTS
+	lc->subcell=malloc(lc->cells*sizeof(int*));
+#else
 	lc->subcell=malloc(lc->cells*sizeof(t_list));
+#endif
+
+	if(lc->subcell==NULL) {
+		perror("[moldyn] cell init (malloc)");
+		return -1;
+	}
 
 	if(lc->cells<27)
 		printf("[moldyn] FATAL: less then 27 subcells!\n");
 
 	if(vol) {
-		printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
+#ifdef STATIC_LISTS
+		printf("[moldyn] initializing 'static' linked cells (%d)\n",
+		       lc->cells);
+#else
+		printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
+		       lc->cells);
+#endif
 		printf("  x: %d x %f A\n",lc->nx,lc->x);
 		printf("  y: %d x %f A\n",lc->ny,lc->y);
 		printf("  z: %d x %f A\n",lc->nz,lc->z);
 	}
 
+#ifdef STATIC_LISTS
+	/* list init */
+	for(i=0;i<lc->cells;i++) {
+		lc->subcell[i]=malloc((MAX_ATOMS_PER_LIST+1)*sizeof(int));
+		if(lc->subcell[i]==NULL) {
+			perror("[moldyn] list init (malloc)");
+			return -1;
+		}
+		/*
+		if(i==0)
+			printf(" ---> %d malloc %p (%p)\n",
+			       i,lc->subcell[0],lc->subcell);
+		*/
+	}
+#else
 	for(i=0;i<lc->cells;i++)
 		list_init_f(&(lc->subcell[i]));
+#endif
+
+	/* update the list */
+	link_cell_update(moldyn);
 
- 	link_cell_update(moldyn);
-	
 	return 0;
 }
 
@@ -1286,7 +1356,9 @@ int link_cell_update(t_moldyn *moldyn) {
 	int nx,ny;
 	t_atom *atom;
 	t_linkcell *lc;
-	double x,y,z;
+#ifdef STATIC_LISTS
+	int p;
+#endif
 
 	atom=moldyn->atom;
 	lc=&(moldyn->lc);
@@ -1294,25 +1366,50 @@ int link_cell_update(t_moldyn *moldyn) {
 	nx=lc->nx;
 	ny=lc->ny;
 
-	x=moldyn->dim.x/2;
-	y=moldyn->dim.y/2;
-	z=moldyn->dim.z/2;
-
 	for(i=0;i<lc->cells;i++)
+#ifdef STATIC_LISTS
+		memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
+#else
 		list_destroy_f(&(lc->subcell[i]));
-	
+#endif
+
 	for(count=0;count<moldyn->count;count++) {
 		i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
 		j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
 		k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
+	
+#ifdef STATIC_LISTS
+		p=0;
+		while(lc->subcell[i+j*nx+k*nx*ny][p]!=0)
+			p++;
+
+		if(p>=MAX_ATOMS_PER_LIST) {
+			printf("[moldyn] FATAL: amount of atoms too high!\n");
+			return -1;
+		}
+
+		lc->subcell[i+j*nx+k*nx*ny][p]=count;
+#else
 		list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
 		                     &(atom[count]));
+		/*
+		if(j==0&&k==0)
+			printf(" ---> %d %d malloc %p (%p)\n",
+			       i,count,lc->subcell[i].current,lc->subcell);
+		*/
+#endif
 	}
 
 	return 0;
 }
 
-int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
+int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
+#ifdef STATIC_LISTS
+                              int **cell
+#else
+                              t_list *cell
+#endif
+			     ) {
 
 	t_linkcell *lc;
 	int a;
@@ -1330,6 +1427,10 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
 	count2=27;
 	a=nx*ny;
 
+	if(i>=nx||j>=ny||k>=nz)
+		printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n",
+		       i,nx,j,ny,k,nz);
+
 	cell[0]=lc->subcell[i+j*nx+k*a];
 	for(ci=-1;ci<=1;ci++) {
 		bx=0;
@@ -1375,8 +1476,14 @@ int link_cell_shutdown(t_moldyn *moldyn) {
 
 	lc=&(moldyn->lc);
 
-	for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
-		list_destroy_f(&(moldyn->lc.subcell[i]));
+	for(i=0;i<lc->cells;i++) {
+#ifdef STATIC_LISTS
+		free(lc->subcell[i]);
+#else
+		//printf(" ---> %d free %p\n",i,lc->subcell[i].start);
+		list_destroy_f(&(lc->subcell[i]));
+#endif
+	}
 
 	free(lc->subcell);
 
@@ -1434,7 +1541,7 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
 int moldyn_integrate(t_moldyn *moldyn) {
 
 	int i;
-	unsigned int e,m,s,v,p,t;
+	unsigned int e,m,s,v,p,t,a;
 	t_3dvec momentum;
 	t_moldyn_schedule *sched;
 	t_atom *atom;
@@ -1456,12 +1563,12 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	m=moldyn->mwrite;
 	s=moldyn->swrite;
 	v=moldyn->vwrite;
+	a=moldyn->awrite;
 	p=moldyn->pwrite;
 	t=moldyn->twrite;
 
 	/* sqaure of some variables */
 	moldyn->tau_square=moldyn->tau*moldyn->tau;
-	moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
 
 	/* get current time */
 	gettimeofday(&t1,NULL);
@@ -1474,18 +1581,19 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
 	/* some stupid checks before we actually start calculating bullshit */
 	if(moldyn->cutoff>0.5*moldyn->dim.x)
-		printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
+		printf("[moldyn] WARNING: cutoff > 0.5 x dim.x\n");
 	if(moldyn->cutoff>0.5*moldyn->dim.y)
-		printf("[moldyn] warning: cutoff > 0.5 x dim.y\n");
+		printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n");
 	if(moldyn->cutoff>0.5*moldyn->dim.z)
-		printf("[moldyn] warning: cutoff > 0.5 x dim.z\n");
+		printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n");
 	ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
 	if(ds>0.05*moldyn->nnd)
-		printf("[moldyn] warning: forces too high / tau too small!\n");
+		printf("[moldyn] WARNING: forces too high / tau too small!\n");
 
 	/* zero absolute time */
-	moldyn->time=0.0;
-	moldyn->total_steps=0;
+	// should have right values!
+	//moldyn->time=0.0;
+	//moldyn->total_steps=0;
 
 	/* debugging, ignore */
 	moldyn->debug=0;
@@ -1517,6 +1625,13 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		temperature_calc(moldyn);
 		virial_sum(moldyn);
 		pressure_calc(moldyn);
+		/*
+		thermodynamic_pressure_calc(moldyn);
+		printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
+		       moldyn->tp/BAR);
+		*/
+
+		/* calculate fluctuations + averages */
 		average_and_fluctuation_calc(moldyn);
 
 		/* p/t scaling */
@@ -1546,9 +1661,10 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		if(p) {
 			if(!(moldyn->total_steps%p)) {
 				dprintf(moldyn->pfd,
-				        "%f %f %f %f %f\n",moldyn->time,
+				        "%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->gp/BAR,moldyn->gp_avg/BAR,
+					 moldyn->tp/BAR,moldyn->tp_avg/BAR);
 			}
 		}
 		if(t) {
@@ -1558,6 +1674,12 @@ int moldyn_integrate(t_moldyn *moldyn) {
 				        moldyn->time,moldyn->t,moldyn->t_avg);
 			}
 		}
+		if(v) {
+			if(!(moldyn->total_steps%v)) {
+				dprintf(moldyn->vfd,
+				        "%f %f\n",moldyn->time,moldyn->volume);
+			}
+		}
 		if(s) {
 			if(!(moldyn->total_steps%s)) {
 				snprintf(dir,128,"%s/s-%07.f.save",
@@ -1573,8 +1695,8 @@ int moldyn_integrate(t_moldyn *moldyn) {
 				close(fd);
 			}	
 		}
-		if(v) {
-			if(!(moldyn->total_steps%v)) {
+		if(a) {
+			if(!(moldyn->total_steps%a)) {
 				visual_atoms(moldyn);
 			}
 		}
@@ -1584,13 +1706,14 @@ int moldyn_integrate(t_moldyn *moldyn) {
 			/* get current time */
 			gettimeofday(&t2,NULL);
 
-	printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
-	       sched->count,i,
-	       moldyn->t,moldyn->t_avg,
-	       moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
-	       moldyn->volume,
-	       (int)(t2.tv_sec-t1.tv_sec));
-	fflush(stdout);
+printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
+       sched->count,i,moldyn->total_steps,
+       moldyn->t,moldyn->t_avg,
+       moldyn->p/BAR,moldyn->p_avg/BAR,
+       moldyn->volume,
+       (int)(t2.tv_sec-t1.tv_sec));
+
+			fflush(stdout);
 
 			/* copy over time */
 			t1=t2;
@@ -1604,8 +1727,8 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
 		/* check for hooks */
 		if(sched->hook) {
-			printf("\n ## schedule hook %d/%d start ##\n",
-			       sched->count+1,sched->total_sched-1);
+			printf("\n ## schedule hook %d start ##\n",
+			       sched->count);
 			sched->hook(moldyn,sched->hook_params);
 			printf(" ## schedule hook end ##\n");
 		}
@@ -1633,6 +1756,9 @@ 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;
 		v3_scale(&delta,&(atom[i].v),tau);
@@ -1646,6 +1772,9 @@ int velocity_verlet(t_moldyn *moldyn) {
 		v3_add(&(atom[i].v),&(atom[i].v),&delta);
 	}
 
+	/* criticial check */
+	moldyn_bc_check(moldyn);
+
 	/* neighbour list update */
 	link_cell_update(moldyn);
 
@@ -1653,6 +1782,9 @@ int velocity_verlet(t_moldyn *moldyn) {
 	potential_force_calc(moldyn);
 
 	for(i=0;i<count;i++) {
+		/* check whether fixed atom */
+		if(atom[i].attr&ATOM_ATTR_FP)
+			continue;
 		/* 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);
@@ -1676,15 +1808,24 @@ int potential_force_calc(t_moldyn *moldyn) {
 	t_atom *itom,*jtom,*ktom;
 	t_virial *virial;
 	t_linkcell *lc;
+#ifdef STATIC_LISTS
+	int *neighbour_i[27];
+	int p,q;
+	t_atom *atom;
+#else
 	t_list neighbour_i[27];
 	t_list neighbour_i2[27];
 	t_list *this,*that;
+#endif
 	u8 bc_ij,bc_ik;
 	int dnlc;
 
 	count=moldyn->count;
 	itom=moldyn->atom;
 	lc=&(moldyn->lc);
+#ifdef STATIC_LISTS
+	atom=moldyn->atom;
+#endif
 
 	/* reset energy */
 	moldyn->energy=0.0;
@@ -1719,8 +1860,8 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 		/* single particle potential/force */
 		if(itom[i].attr&ATOM_ATTR_1BP)
-			if(moldyn->func1b)
-				moldyn->func1b(moldyn,&(itom[i]));
+			if(moldyn->func_i0)
+				moldyn->func_i0(moldyn,&(itom[i]));
 
 		if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
 			continue;
@@ -1735,57 +1876,152 @@ 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;
+#ifdef STATIC_LISTS
+				p=0;
+
+				while(neighbour_i[j][p]!=0) {
+
+					jtom=&(atom[neighbour_i[j][p]]);
+					p++;
+#else
 				this=&(neighbour_i[j]);
 				list_reset_f(this);
 
 				if(this->start==NULL)
 					continue;
 
-				bc_ij=(j<dnlc)?0:1;
-
 				do {
 					jtom=this->current->data;
+#endif
 
 					if(jtom==&(itom[i]))
 						continue;
 
+					/* reset 3bp run */
+					moldyn->run3bp=1;
+
 					if((jtom->attr&ATOM_ATTR_2BP)&
 					   (itom[i].attr&ATOM_ATTR_2BP)) {
-						moldyn->func2b(moldyn,
-						               &(itom[i]),
-						               jtom,
-					        	       bc_ij);
+						moldyn->func_j0(moldyn,
+						                &(itom[i]),
+						                jtom,
+					        	        bc_ij);
 					}
+
+					/* 3 body potential/force */
+
+					/* in j loop, 3bp run can be skipped */
+					if(!(moldyn->run3bp))
+						continue;
+
+					if(!(itom[i].attr&ATOM_ATTR_3BP))
+						continue;
+
+					if(!(jtom->attr&ATOM_ATTR_3BP))
+						continue;
+
+					if(moldyn->func_j0_k0==NULL)
+						continue;
+
+				/* first loop over atoms k in first j loop */
+				for(k=0;k<27;k++) {
+
+					bc_ik=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+					q=0;
+
+					while(neighbour_i[j][q]!=0) {
+
+						ktom=&(atom[neighbour_i[k][q]]);
+						q++;
+#else
+					that=&(neighbour_i2[k]);
+					list_reset_f(that);
+					
+					if(that->start==NULL)
+						continue;
+
+					do {
+						ktom=that->current->data;
+#endif
+
+						if(!(ktom->attr&ATOM_ATTR_3BP))
+							continue;
+
+						//if(ktom==jtom)
+						//	continue;
+
+						if(ktom==&(itom[i]))
+							continue;
+
+						moldyn->func_j0_k0(moldyn,
+						                   &(itom[i]),
+						                   jtom,
+						                   ktom,
+						                   bc_ik|bc_ij);
+#ifdef STATIC_LISTS
+					}
+#else
+					} while(list_next_f(that)!=\
+					        L_NO_NEXT_ELEMENT);
+#endif
+
+				}
+
+				/* finish of first j loop after first k loop */
+				if(moldyn->func_j0e)
+					moldyn->func_j0e(moldyn,
+					                 &(itom[i]),
+					                 jtom,
+					                 bc_ij);
+
+#ifdef STATIC_LISTS
+				}
+#else
 				} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
 
 			}
 		}
 
-		/* 3 body potential/force */
+		/* continued 3 body potential/force */
 
 		if(!(itom[i].attr&ATOM_ATTR_3BP))
 			continue;
 
-		/* copy the neighbour lists */
-		memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
-
 		/* second loop over atoms j */
 		for(j=0;j<27;j++) {
 
+			bc_ij=(j<dnlc)?0:1;
+#ifdef STATIC_LISTS
+			p=0;
+
+			while(neighbour_i[j][p]!=0) {
+
+				jtom=&(atom[neighbour_i[j][p]]);
+				p++;
+#else
 			this=&(neighbour_i[j]);
 			list_reset_f(this);
 
 			if(this->start==NULL)
 				continue;
 
-			bc_ij=(j<dnlc)?0:1;
-
 			do {
+
 				jtom=this->current->data;
+#endif
 
 				if(jtom==&(itom[i]))
 					continue;
@@ -1796,111 +2032,136 @@ 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++) {
 
+					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;
 
-					bc_ik=(k<dnlc)?0:1;
-
 					do {
-
 						ktom=that->current->data;
+#endif
 
 						if(!(ktom->attr&ATOM_ATTR_3BP))
 							continue;
 
-						if(ktom==jtom)
-							continue;
+						//if(ktom==jtom)
+						//	continue;
 
 						if(ktom==&(itom[i]))
 							continue;
 
-						moldyn->func3b_k1(moldyn,
-						                  &(itom[i]),
-						                  jtom,
-						                  ktom,
-						                  bc_ik|bc_ij);
-
+						moldyn->func_j1_k0(moldyn,
+						                   &(itom[i]),
+						                   jtom,
+						                   ktom,
+						                   bc_ik|bc_ij);
+#ifdef STATIC_LISTS
+					}
+#else
 					} while(list_next_f(that)!=\
 					        L_NO_NEXT_ELEMENT);
+#endif
 
 				}
 
 				}
 
-				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++) {
 
+					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;
 
-					bc_ik=(k<dnlc)?0:1;
-
 					do {
-
 						ktom=that->current->data;
+#endif
 
 						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
+					}
+#else
 					} while(list_next_f(that)!=\
 					        L_NO_NEXT_ELEMENT);
+#endif
 
 				}
 				
 				}
 
-				/* 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
+			}
+#else
 			} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
 		
 		}
-		
+
 #ifdef DEBUG
 	//printf("\n\n");
 #endif
@@ -1914,20 +2175,27 @@ int potential_force_calc(t_moldyn *moldyn) {
 	//printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
 	if(moldyn->time>DSTART&&moldyn->time<DEND) {
 		printf("force:\n");
-		printf("  x: %0.40f\n",moldyn->atom[5832].f.x);
-		printf("  y: %0.40f\n",moldyn->atom[5832].f.y);
-		printf("  z: %0.40f\n",moldyn->atom[5832].f.z);
+		printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
+		printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
+		printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
 	}
 #endif
 
-	/* calculate global virial */
+	/* some postprocessing */
 	for(i=0;i<count;i++) {
-		moldyn->gvir.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x;
-		moldyn->gvir.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y;
-		moldyn->gvir.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z;
-		moldyn->gvir.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x;
-		moldyn->gvir.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x;
-		moldyn->gvir.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y;
+		/* calculate global virial */
+		moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
+		moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
+		moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
+		moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
+		moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
+		moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
+
+		/* check forces regarding the given timestep */
+		if(v3_norm(&(itom[i].f))>\
+		   0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
+			printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
+			       i);
 	}
 
 	return 0;
@@ -2042,6 +2310,8 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 
 	int fd;
 	int cnt,size;
+	int fsize;
+	int corr;
 
 	fd=open(file,O_RDONLY);
 	if(fd<0) {
@@ -2049,21 +2319,58 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 		return fd;
 	}
 
+	fsize=lseek(fd,0,SEEK_END);
+	lseek(fd,0,SEEK_SET);
+
 	size=sizeof(t_moldyn);
-	cnt=read(fd,moldyn,size);
-	if(cnt!=size) {
-		perror("[moldyn] load save file read (moldyn)");
-		return cnt;
+
+	while(size) {
+		cnt=read(fd,moldyn,size);
+		if(cnt<0) {
+			perror("[moldyn] load save file read (moldyn)");
+			return cnt;
+		}
+		size-=cnt;
 	}
 
 	size=moldyn->count*sizeof(t_atom);
-	cnt=read(fd,moldyn->atom,size);
-	if(cnt!=size) {
-		perror("[moldyn] load save file read (atoms)");
-		return cnt;
+
+	/* correcting possible atom data offset */
+	corr=0;
+	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("  - current pos: %d\n",sizeof(t_moldyn));
+		printf("  - atom size: %d\n",size);
+		printf("  - file size: %d\n",fsize);
+		printf("  => correction: %d\n",corr);
+		lseek(fd,corr,SEEK_CUR);
+	}
+
+	moldyn->atom=(t_atom *)malloc(size);
+	if(moldyn->atom==NULL) {
+		perror("[moldyn] load save file malloc (atoms)");
+		return -1;
+	}
+
+	while(size) {
+		cnt=read(fd,moldyn->atom,size);
+		if(cnt<0) {
+			perror("[moldyn] load save file read (atoms)");
+			return cnt;
+		}
+		size-=cnt;
 	}
 
-	// hooks
+	// hooks etc ...
+
+	return 0;
+}
+
+int moldyn_free_save_file(t_moldyn *moldyn) {
+
+	free(moldyn->atom);
 
 	return 0;
 }
@@ -2075,6 +2382,75 @@ int moldyn_load(t_moldyn *moldyn) {
 	return 0;
 }
 
+/*
+ * function to find/callback all combinations of 2 body bonds
+ */
+
+int process_2b_bonds(t_moldyn *moldyn,void *data,
+                     int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+		                    void *data,u8 bc)) {
+
+	t_linkcell *lc;
+#ifdef STATIC_LISTS
+	int *neighbour[27];
+	int p;
+#else
+	t_list neighbour[27];
+#endif
+	u8 bc;
+	t_atom *itom,*jtom;
+	int i,j;
+	t_list *this;
+
+	lc=&(moldyn->lc);
+	itom=moldyn->atom;
+	
+	for(i=0;i<moldyn->count;i++) {
+		/* neighbour indexing */
+		link_cell_neighbour_index(moldyn,
+		                          (itom[i].r.x+moldyn->dim.x/2)/lc->x,
+		                          (itom[i].r.y+moldyn->dim.y/2)/lc->x,
+		                          (itom[i].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]!=0) {
+
+				jtom=&(moldyn->atom[neighbour[j][p]]);
+				p++;
+#else
+			this=&(neighbour[j]);
+			list_reset_f(this);
+
+			if(this->start==NULL)
+				continue;
+
+			do {
+
+				jtom=this->current->data;
+#endif
+
+				/* process bond */
+				process(moldyn,&(itom[i]),jtom,data,bc);
+
+#ifdef STATIC_LISTS
+			}
+#else
+			} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
+		}
+	}
+
+	return 0;
+
+}
+
 /*
  * post processing functions
  */
@@ -2090,7 +2466,8 @@ int get_line(int fd,char *line,int max) {
 		ret=read(fd,line+count,1);
 		if(ret<=0) return ret;
 		if(line[count]=='\n') {
-			line[count]='\0';
+			memset(line+count,0,max-count-1);
+			//line[count]='\0';
 			return count+1;
 		}
 		count+=1;
@@ -2103,58 +2480,267 @@ int pair_correlation_init(t_moldyn *moldyn,double dr) {
 	return 0;
 }
 
-int calculate_pair_correlation(t_moldyn *moldyn,double dr) {
+int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 
-	int slots;
-	int *stat;
 	int i;
-	t_linkcell *lc;
-	t_list neighbour[27];
-	t_atom *itom,*jtom;
-	t_list *this;
+	t_atom *atom;
+	t_3dvec dist;
+	double d2;
+	int a_cnt;
+	int b_cnt;
 
-	lc=&(moldyn->lc);
+	atom=moldyn->atom;
+	dc[0]=0;
+	dc[1]=0;
+	dc[2]=0;
+	a_cnt=0;
+	b_cnt=0;
 
-	slots=(int)(moldyn->cutoff/dr);
+	for(i=0;i<moldyn->count;i++) {
 
-	stat=(int *)malloc(3*slots*sizeof(int));
-	if(stat==NULL) {
-		perror("[moldyn] pair correlation malloc");
-		return -1;
+		v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
+		check_per_bound(moldyn,&dist);
+		d2=v3_absolute_square(&dist);
+
+		if(atom[i].brand) {
+			b_cnt+=1;
+			dc[1]+=d2;
+		}
+		else {
+			a_cnt+=1;
+			dc[0]+=d2;
+		}
+
+		dc[2]+=d2;
 	}
 
-	link_cell_init(moldyn,VERBOSE);
-	
-	for(i=0;i<lc->cells;i++) {
-		/* check for atoms */
-		itom=lc->subcell[i].current->data;
-		if(itom==NULL)
-			continue;
+	dc[0]*=(1.0/(6.0*moldyn->time*a_cnt));
+	dc[1]*=(1.0/(6.0*moldyn->time*b_cnt));
+	dc[2]*=(1.0/(6.0*moldyn->time*moldyn->count));
+		
+	return 0;
+}
 
-		/* pick first atom and do neighbour indexing */
-		link_cell_neighbour_index(moldyn,
-		                          (itom->r.x+moldyn->dim.x/2)/lc->x,
-		                          (itom->r.y+moldyn->dim.y/2)/lc->x,
-		                          (itom->r.z+moldyn->dim.z/2)/lc->x,
-		                          neighbour);
+int bonding_analyze(t_moldyn *moldyn,double *cnt) {
+
+	return 0;
+}
 
-		/* calculation of g(r) */
-		do {
-			itom=lc->subcell[i].current->data;
-// GO ON HERE ...
-		} while(list_next_f(this)!=L_NO_NEXT_ELEMENT););
+int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
+                                       t_atom *jtom,void *data,u8 bc) {
+
+	t_3dvec dist;
+	double d;
+	int s;
+	t_pcc *pcc;
+
+	/* only count pairs once,
+	 * skip same atoms */
+	if(itom->tag>=jtom->tag)
+		return 0;
+
+	/*
+	 * pair correlation calc
+	 */
 
+	/* get pcc data */
+	pcc=data;
+
+	/* distance */
+	v3_sub(&dist,&(jtom->r),&(itom->r));
+	if(bc) check_per_bound(moldyn,&dist);
+	d=v3_absolute_square(&dist);
+
+	/* ignore if greater cutoff */
+	if(d>moldyn->cutoff_square)
+		return 0;
+
+	/* fill the slots */
+	d=sqrt(d);
+	s=(int)(d/pcc->dr);
+
+	/* should never happen but it does 8) -
+	 * related to -ffloat-store problem! */
+	if(s>=pcc->o1) {
+		printf("[moldyn] WARNING: pcc (%d/%d)",
+		       s,pcc->o1);
+		printf("\n");
+		s=pcc->o1-1;
 	}
 
-	link_cell_shutdown(moldyn);
+	if(itom->brand!=jtom->brand) {
+		/* mixed */
+		pcc->stat[s]+=1;
+	}
+	else {
+		/* type a - type a bonds */
+		if(itom->brand==0)
+			pcc->stat[s+pcc->o1]+=1;
+		else
+		/* type b - type b bonds */
+			pcc->stat[s+pcc->o2]+=1;
+	}
 
 	return 0;
 }
 
-int analyze_bonds(t_moldyn *moldyn) {
+int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
 
+	t_pcc pcc;
+	double norm;
+	int i;
+
+	pcc.dr=dr;
+	pcc.o1=moldyn->cutoff/dr;
+	pcc.o2=2*pcc.o1;
+
+	if(pcc.o1*dr<=moldyn->cutoff)
+		printf("[moldyn] WARNING: pcc (low #slots)\n");
+
+	printf("[moldyn] pair correlation calc info:\n");
+	printf("  time: %f\n",moldyn->time);
+	printf("  count: %d\n",moldyn->count);
+	printf("  cutoff: %f\n",moldyn->cutoff);
+	printf("  temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg);
+
+	if(ptr!=NULL) {
+		pcc.stat=(double *)ptr;
+	}
+	else {
+		pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double));
+		if(pcc.stat==NULL) {
+			perror("[moldyn] pair correlation malloc");
+			return -1;
+		}
+	}
+
+	memset(pcc.stat,0,3*pcc.o1*sizeof(double));
+
+	/* process */
+	process_2b_bonds(moldyn,&pcc,calculate_pair_correlation_process);
+
+	/* normalization */
+	for(i=1;i<pcc.o1;i++) {
+		 // normalization: 4 pi r^2 dr
+		 // here: not double counting pairs -> 2 pi r r dr
+		 // ... and actually it's a constant times r^2
+		norm=i*i*dr*dr;
+		pcc.stat[i]/=norm;
+		pcc.stat[pcc.o1+i]/=norm;
+		pcc.stat[pcc.o2+i]/=norm;
+	}
+	/* */
+
+	if(ptr==NULL) {
+		/* todo: store/print pair correlation function */
+		free(pcc.stat);
+	}
+
+	return 0;
+}
+
+int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+                         void *data,u8 bc) {
+
+	t_ba *ba;
+	t_3dvec dist;
+	double d;
+
+	if(itom->tag>=jtom->tag)
+		return 0;
+
+	/* distance */
+	v3_sub(&dist,&(jtom->r),&(itom->r));
+	if(bc) check_per_bound(moldyn,&dist);
+	d=v3_absolute_square(&dist);
+
+	/* ignore if greater or equal cutoff */
+	if(d>moldyn->cutoff_square)
+		return 0;
+
+	/* check for potential bond */
+	if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
+		return 0;
+
+	/* now count this bonding ... */
+	ba=data;
+
+	/* increase total bond counter
+	 * ... double counting!
+	 */
+	ba->tcnt+=2;
+
+	if(itom->brand==0)
+		ba->acnt[jtom->tag]+=1;
+	else
+		ba->bcnt[jtom->tag]+=1;
 	
-	
+	if(jtom->brand==0)
+		ba->acnt[itom->tag]+=1;
+	else
+		ba->bcnt[itom->tag]+=1;
+
+	return 0;
+}
+
+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;
+	t_ba ba;
+	int i;
+	t_atom *atom;
+
+	ba.acnt=malloc(moldyn->count*sizeof(int));
+	if(ba.acnt==NULL) {
+		perror("[moldyn] bond analyze malloc (a)");
+		return -1;
+	}
+	memset(ba.acnt,0,moldyn->count*sizeof(int));
+
+	ba.bcnt=malloc(moldyn->count*sizeof(int));
+	if(ba.bcnt==NULL) {
+		perror("[moldyn] bond analyze malloc (b)");
+		return -1;
+	}
+	memset(ba.bcnt,0,moldyn->count*sizeof(int));
+
+	ba.tcnt=0;
+	qcnt=0;
+	ccnt=0;
+	cset=0;
+
+	atom=moldyn->atom;
+
+	process_2b_bonds(moldyn,&ba,bond_analyze_process);
+
+	for(i=0;i<moldyn->count;i++) {
+		if(atom[i].brand==0) {
+			if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
+				qcnt+=4;
+		}
+		else {
+			if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
+				qcnt+=4;
+				ccnt+=1;
+			}
+			cset+=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;
+	}
+	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);
+	}
 
 	return 0;
 }
@@ -2170,34 +2756,48 @@ int visual_init(t_moldyn *moldyn,char *filebase) {
 	return 0;
 }
 
+int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+                         void *data,u8 bc) {
+
+	t_vb *vb;
+
+	vb=data;
+
+	if(itom->tag>=jtom->tag)
+		return 0;
+	
+	if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
+		return 0;
+
+	if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB))
+		dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n",
+		        itom->r.x,itom->r.y,itom->r.z,
+		        jtom->r.x,jtom->r.y,jtom->r.z);
+
+	return 0;
+}
+
 int visual_atoms(t_moldyn *moldyn) {
 
-	int i,j,fd;
+	int i;
  	char file[128+64];
 	t_3dvec dim;
 	double help;
 	t_visual *v;
 	t_atom *atom;
-	t_atom *btom;
-	t_linkcell *lc;
-	t_list neighbour[27];
-	u8 bc;
-	t_3dvec dist;
-	double d2;
-	u8 brand;
+	t_vb vb;
 
 	v=&(moldyn->vis);
 	dim.x=v->dim.x;
 	dim.y=v->dim.y;
 	dim.z=v->dim.z;
 	atom=moldyn->atom;
-	lc=&(moldyn->lc);
 
 	help=(dim.x+dim.y);
 
 	sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
-	fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
-	if(fd<0) {
+	vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+	if(vb.fd<0) {
 		perror("open visual save file fd");
 		return -1;
 	}
@@ -2205,107 +2805,65 @@ int visual_atoms(t_moldyn *moldyn) {
 	/* write the actual data file */
 
 	// povray header
-	dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n",
+	dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n",
 	        moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
 
 	// atomic configuration
-	for(i=0;i<moldyn->count;i++) {
+	for(i=0;i<moldyn->count;i++)
 		// atom type, positions, color and kinetic energy
-		dprintf(fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
-		                                 atom[i].r.x,
-		                                 atom[i].r.y,
-		                                 atom[i].r.z,
-		                                 pse_col[atom[i].element],
-		                                 atom[i].ekin);
-
-		/*
-		 * bond detection should usually be done by potential
-		 * functions. brrrrr! EVIL!
-		 * 
-		 * todo: potentials need to export a 'find_bonds' function!
-		 */
-
-		// bonds between atoms
-		if(!(atom[i].attr&ATOM_ATTR_VB))
-			continue;
-		link_cell_neighbour_index(moldyn,
-		                          (atom[i].r.x+moldyn->dim.x/2)/lc->x,
-		                          (atom[i].r.y+moldyn->dim.y/2)/lc->y,
-		                          (atom[i].r.z+moldyn->dim.z/2)/lc->z,
-		                          neighbour);
-		for(j=0;j<27;j++) {
-			list_reset_f(&neighbour[j]);
-			if(neighbour[j].start==NULL)
-				continue;
-			bc=j<lc->dnlc?0:1;
-			do {
-				btom=neighbour[j].current->data;
-				if(btom==&atom[i])	// skip identical atoms
-					continue;
-				//if(btom<&atom[i])	// skip half of them
-				//	continue;
-				v3_sub(&dist,&(atom[i].r),&(btom->r));
-				if(bc) check_per_bound(moldyn,&dist);
-				d2=v3_absolute_square(&dist);
-				brand=atom[i].brand;
-				if(brand==btom->brand) {
-					if(d2>moldyn->bondlen[brand])
-						continue;
-				}
-				else {
-					if(d2>moldyn->bondlen[2])
-						continue;
-				}
-				dprintf(fd,"# [B] %f %f %f %f %f %f\n",
-				        atom[i].r.x,atom[i].r.y,atom[i].r.z,
-		       		        btom->r.x,btom->r.y,btom->r.z);
-			} while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT);
-		}
-	}
-
+		dprintf(vb.fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
+		                                    atom[i].r.x,
+		                                    atom[i].r.y,
+		                                    atom[i].r.z,
+		                                    pse_col[atom[i].element],
+		                                    atom[i].ekin);
+	
+	// bonds between atoms
+	process_2b_bonds(moldyn,&vb,visual_bonds_process);
+	
 	// boundaries
 	if(dim.x) {
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,-dim.y/2,-dim.z/2,
 		        dim.x/2,-dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,-dim.y/2,-dim.z/2,
 		        -dim.x/2,dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        dim.x/2,dim.y/2,-dim.z/2,
 		        dim.x/2,-dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,dim.y/2,-dim.z/2,
 		        dim.x/2,dim.y/2,-dim.z/2);
 
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,-dim.y/2,dim.z/2,
 		        dim.x/2,-dim.y/2,dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,-dim.y/2,dim.z/2,
 		        -dim.x/2,dim.y/2,dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        dim.x/2,dim.y/2,dim.z/2,
 		        dim.x/2,-dim.y/2,dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,dim.y/2,dim.z/2,
 		        dim.x/2,dim.y/2,dim.z/2);
 
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,-dim.y/2,dim.z/2,
 		        -dim.x/2,-dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        -dim.x/2,dim.y/2,dim.z/2,
 		        -dim.x/2,dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        dim.x/2,-dim.y/2,dim.z/2,
 		        dim.x/2,-dim.y/2,-dim.z/2);
-		dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+		dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
 		        dim.x/2,dim.y/2,dim.z/2,
 		        dim.x/2,dim.y/2,-dim.z/2);
 	}
 
-	close(fd);
+	close(vb.fd);
 
 	return 0;
 }