X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=4d736179550d68c7402e89b4412dec3dfc9d40bb;hb=0b96eb313c9bfec6272b1f8de0d99c4ce26d1686;hp=1ab88eb98a4a02c91ac7ef367c99bab2b2f7407a;hpb=caa3bc828974c35df2462fde737c31c0a618ee4e;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 1ab88eb..4d73617 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -82,7 +82,7 @@ int set_pressure(t_moldyn *moldyn,double p_ref) {
 
 	moldyn->p_ref=p_ref;
 
-	printf("[moldyn] pressure [atm]: %f\n",moldyn->p_ref/ATM);
+	printf("[moldyn] pressure [bar]: %f\n",moldyn->p_ref/BAR);
 
 	return 0;
 }
@@ -162,34 +162,58 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
 	return 0;
 }
 
-int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params) {
+int set_potential1b(t_moldyn *moldyn,pf_func1b func) {
 
 	moldyn->func1b=func;
-	moldyn->pot1b_params=params;
 
 	return 0;
 }
 
-int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) {
+int set_potential2b(t_moldyn *moldyn,pf_func2b func) {
 
 	moldyn->func2b=func;
-	moldyn->pot2b_params=params;
 
 	return 0;
 }
 
-int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) {
+int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) {
 
-	moldyn->func2b_post=func;
-	moldyn->pot2b_params=params;
+	moldyn->func3b_j1=func;
 
 	return 0;
 }
 
-int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) {
+int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) {
 
-	moldyn->func3b=func;
-	moldyn->pot3b_params=params;
+	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) {
+
+	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;
 
 	return 0;
 }
@@ -243,6 +267,32 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 			dprintf(moldyn->efd,"# total momentum log file\n");
 			printf("total momentum (%d)\n",timer);
 			break;
+		case LOG_PRESSURE:
+			moldyn->pwrite=timer;
+			snprintf(filename,127,"%s/pressure",moldyn->vlsdir);
+			moldyn->pfd=open(filename,
+			                 O_WRONLY|O_CREAT|O_EXCL,
+			                 S_IRUSR|S_IWUSR);
+			if(moldyn->pfd<0) {
+				perror("[moldyn] pressure log file\n");
+				return moldyn->pfd;
+			}
+			dprintf(moldyn->pfd,"# pressure log file\n");
+			printf("pressure (%d)\n",timer);
+			break;
+		case LOG_TEMPERATURE:
+			moldyn->twrite=timer;
+			snprintf(filename,127,"%s/temperature",moldyn->vlsdir);
+			moldyn->tfd=open(filename,
+			                 O_WRONLY|O_CREAT|O_EXCL,
+			                 S_IRUSR|S_IWUSR);
+			if(moldyn->tfd<0) {
+				perror("[moldyn] temperature log file\n");
+				return moldyn->tfd;
+			}
+			dprintf(moldyn->tfd,"# temperature log file\n");
+			printf("temperature (%d)\n",timer);
+			break;
 		case SAVE_STEP:
 			moldyn->swrite=timer;
 			printf("save file (%d)\n",timer);
@@ -265,18 +315,52 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 				perror("[moldyn] report fd open");	
 				return moldyn->rfd;
 			}
-			snprintf(filename,127,"%s/plot.scr",moldyn->vlsdir);
-			moldyn->pfd=open(filename,
-			                 O_WRONLY|O_CREAT|O_EXCL,
-			                 S_IRUSR|S_IWUSR);
-			if(moldyn->pfd<0) {
-				perror("[moldyn] plot fd open");
-				return moldyn->pfd;
+			printf("report -> ");
+			if(moldyn->efd) {
+				snprintf(filename,127,"%s/e_plot.scr",
+				         moldyn->vlsdir);
+				moldyn->epfd=open(filename,
+			 	                 O_WRONLY|O_CREAT|O_EXCL,
+			  	                 S_IRUSR|S_IWUSR);
+				if(moldyn->epfd<0) {
+					perror("[moldyn] energy plot fd open");
+					return moldyn->epfd;
+				}
+				dprintf(moldyn->epfd,e_plot_script);
+				close(moldyn->epfd);
+				printf("energy ");
+			}
+			if(moldyn->pfd) {
+				snprintf(filename,127,"%s/pressure_plot.scr",
+				         moldyn->vlsdir);
+				moldyn->ppfd=open(filename,
+			 	                  O_WRONLY|O_CREAT|O_EXCL,
+			  	                  S_IRUSR|S_IWUSR);
+				if(moldyn->ppfd<0) {
+					perror("[moldyn] p plot fd open");
+					return moldyn->ppfd;
+				}
+				dprintf(moldyn->ppfd,pressure_plot_script);
+				close(moldyn->ppfd);
+				printf("pressure ");
+			}
+			if(moldyn->tfd) {
+				snprintf(filename,127,"%s/temperature_plot.scr",
+				         moldyn->vlsdir);
+				moldyn->tpfd=open(filename,
+			 	                  O_WRONLY|O_CREAT|O_EXCL,
+			  	                  S_IRUSR|S_IWUSR);
+				if(moldyn->tpfd<0) {
+					perror("[moldyn] t plot fd open");
+					return moldyn->tpfd;
+				}
+				dprintf(moldyn->tpfd,temperature_plot_script);
+				close(moldyn->tpfd);
+				printf("temperature ");
 			}
 			dprintf(moldyn->rfd,report_start,
 			        moldyn->rauthor,moldyn->rtitle);
-			dprintf(moldyn->pfd,plot_script);
-			close(moldyn->pfd);
+			printf("\n");
 			break;
 		default:
 			printf("unknown log type: %02x\n",type);
@@ -291,13 +375,35 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
 	char sc[256];
 
 	printf("[moldyn] log shutdown\n");
-	if(moldyn->efd) close(moldyn->efd);
+	if(moldyn->efd) {
+		close(moldyn->efd);
+		if(moldyn->rfd) {
+			dprintf(moldyn->rfd,report_energy);
+			snprintf(sc,255,"cd %s && gnuplot e_plot.scr",
+			         moldyn->vlsdir);
+			system(sc);
+		}
+	}
 	if(moldyn->mfd) close(moldyn->mfd);
+	if(moldyn->pfd) {
+		close(moldyn->pfd);
+		if(moldyn->rfd)
+			dprintf(moldyn->rfd,report_pressure);
+			snprintf(sc,255,"cd %s && gnuplot pressure_plot.scr",
+			         moldyn->vlsdir);
+			system(sc);
+	}
+	if(moldyn->tfd) {
+		close(moldyn->tfd);
+		if(moldyn->rfd)
+			dprintf(moldyn->rfd,report_temperature);
+			snprintf(sc,255,"cd %s && gnuplot temperature_plot.scr",
+			         moldyn->vlsdir);
+			system(sc);
+	}
 	if(moldyn->rfd) {
 		dprintf(moldyn->rfd,report_end);
 		close(moldyn->rfd);
-		snprintf(sc,255,"cd %s && gnuplot plot.scr",moldyn->vlsdir);
-		system(sc);
 		snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir);
 		system(sc);
 		snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir);
@@ -592,6 +698,8 @@ 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);
+	moldyn->t_sum+=moldyn->t;
+	moldyn->mean_t=moldyn->t_sum/moldyn->total_steps;
 
 	return moldyn->t;
 }
@@ -673,9 +781,11 @@ double pressure_calc(t_moldyn *moldyn) {
 	t_virial *virial;
 
 	/*
-	 * P = 1/(3V) sum_i ( p_i^2 / 2m + f_i r_i )
-	 *
-	 * virial = f_i r_i
+	 * PV = NkT + <W>
+	 * W = 1/3 sum_i f_i r_i
+	 * virial = sum_i f_i r_i
+	 * 
+	 * => P = (2 Ekin + virial) / (3V)
 	 */
 
 	v=0.0;
@@ -687,6 +797,16 @@ double pressure_calc(t_moldyn *moldyn) {
 	/* assume up to date kinetic energy */
 	moldyn->p=2.0*moldyn->ekin+v;
 	moldyn->p/=(3.0*moldyn->volume);
+	moldyn->p_sum+=moldyn->p;
+	moldyn->mean_p=moldyn->p_sum/moldyn->total_steps;
+
+	/* pressure from 'absolute coordinates' virial */
+	virial=&(moldyn->virial);
+	v=virial->xx+virial->yy+virial->zz;
+	moldyn->gp=2.0*moldyn->ekin+v;
+	moldyn->gp/=(3.0*moldyn->volume);
+	moldyn->gp_sum+=moldyn->gp;
+	moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps;
 
 	return moldyn->p;
 }	
@@ -695,9 +815,20 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 
 	t_3dvec dim,*tp;
 	double u,p;
-	double scale;
+	double scale,dv;
 	t_atom *store;
 
+	/*
+	 * dU = - p dV
+	 *
+	 * => p = - dU/dV
+	 *
+	 * dV: dx,y,z = 0.001 x,y,z
+	 */
+
+	scale=1.00000000000001;
+printf("\n\nP-DEBUG:\n");
+
 	tp=&(moldyn->tp);
 	store=malloc(moldyn->count*sizeof(t_atom));
 	if(store==NULL) {
@@ -711,13 +842,13 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	dim=moldyn->dim;
 
 	/* derivative with respect to x direction */
-	scale=1.0+moldyn->dv/(moldyn->dim.y*moldyn->dim.z);
 	scale_dim(moldyn,scale,TRUE,0,0);
 	scale_atoms(moldyn,scale,TRUE,0,0);
+	dv=0.00000000000001*moldyn->dim.x*moldyn->dim.y*moldyn->dim.z;
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
 	potential_force_calc(moldyn);
-	tp->x=(moldyn->energy-u)/moldyn->dv;
+	tp->x=(moldyn->energy-u)/dv;
 	p=tp->x*tp->x;
 
 	/* restore atomic configuration + dim */
@@ -725,13 +856,13 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	moldyn->dim=dim;
 
 	/* derivative with respect to y direction */
-	scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.z);
 	scale_dim(moldyn,scale,0,TRUE,0);
 	scale_atoms(moldyn,scale,0,TRUE,0);
+	dv=0.00000000000001*moldyn->dim.y*moldyn->dim.x*moldyn->dim.z;
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
 	potential_force_calc(moldyn);
-	tp->y=(moldyn->energy-u)/moldyn->dv;
+	tp->y=(moldyn->energy-u)/dv;
 	p+=tp->y*tp->y;
 
 	/* restore atomic configuration + dim */
@@ -739,37 +870,19 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	moldyn->dim=dim;
 
 	/* derivative with respect to z direction */
-	scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.y);
 	scale_dim(moldyn,scale,0,0,TRUE);
 	scale_atoms(moldyn,scale,0,0,TRUE);
+	dv=0.00000000000001*moldyn->dim.z*moldyn->dim.x*moldyn->dim.y;
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
 	potential_force_calc(moldyn);
-	tp->z=(moldyn->energy-u)/moldyn->dv;
+	tp->z=(moldyn->energy-u)/dv;
 	p+=tp->z*tp->z;
 
 	/* restore atomic configuration + dim */
 	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 	moldyn->dim=dim;
 
-	printf("dU/dV komp addiert = %f %f %f\n",tp->x,tp->y,tp->z);
-
-	scale=1.0+pow(moldyn->dv/moldyn->volume,ONE_THIRD);
-
-printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x);
-	scale_dim(moldyn,scale,1,1,1);
-	scale_atoms(moldyn,scale,1,1,1);
-	link_cell_shutdown(moldyn);
-	link_cell_init(moldyn,QUIET);
-	potential_force_calc(moldyn);
-printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x);
-
-	printf("dU/dV einfach = %f\n",((moldyn->energy-u)/moldyn->dv)/ATM);
-
-	/* restore atomic configuration + dim */
-	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
-	moldyn->dim=dim;
-
 	/* restore energy */
 	moldyn->energy=u;
 
@@ -863,7 +976,7 @@ moldyn->debug=scale;
 
 }
 
-double get_e_kin(t_moldyn *moldyn) {
+double e_kin_calc(t_moldyn *moldyn) {
 
 	int i;
 	t_atom *atom;
@@ -877,11 +990,6 @@ double get_e_kin(t_moldyn *moldyn) {
 	return moldyn->ekin;
 }
 
-double update_e_kin(t_moldyn *moldyn) {
-
-	return(get_e_kin(moldyn));
-}
-
 double get_total_energy(t_moldyn *moldyn) {
 
 	return(moldyn->ekin+moldyn->energy);
@@ -942,7 +1050,12 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 	if(lc->cells<27)
 		printf("[moldyn] FATAL: less then 27 subcells!\n");
 
-	if(vol) printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
+	if(vol) {
+		printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
+		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);
+	}
 
 	for(i=0;i<lc->cells;i++)
 		list_init_f(&(lc->subcell[i]));
@@ -977,7 +1090,7 @@ int link_cell_update(t_moldyn *moldyn) {
 		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);
-		list_add_immediate_f(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]),
+		list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
 		                     &(atom[count]));
 	}
 
@@ -1106,14 +1219,15 @@ 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;
-	t_3dvec p;
+	unsigned int e,m,s,v,p,t;
+	t_3dvec momentum;
 	t_moldyn_schedule *sched;
 	t_atom *atom;
 	int fd;
 	char dir[128];
 	double ds;
 	double energy_scale;
+	//double tp;
 
 	sched=&(moldyn->schedule);
 	atom=moldyn->atom;
@@ -1126,6 +1240,8 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	m=moldyn->mwrite;
 	s=moldyn->swrite;
 	v=moldyn->vwrite;
+	p=moldyn->pwrite;
+	t=moldyn->twrite;
 
 	/* sqaure of some variables */
 	moldyn->tau_square=moldyn->tau*moldyn->tau;
@@ -1150,6 +1266,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
 	/* zero absolute time */
 	moldyn->time=0.0;
+	moldyn->total_steps=0;
 
 	/* debugging, ignore */
 	moldyn->debug=0;
@@ -1173,10 +1290,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		moldyn->integrate(moldyn);
 
 		/* calculate kinetic energy, temperature and pressure */
-		update_e_kin(moldyn);
+		e_kin_calc(moldyn);
 		temperature_calc(moldyn);
 		pressure_calc(moldyn);
-		//thermodynamic_pressure_calc(moldyn);
+		//tp=thermodynamic_pressure_calc(moldyn);
+//printf("thermodynamic p: %f %f %f - %f\n",moldyn->tp.x/BAR,moldyn->tp.y/BAR,moldyn->tp.z/BAR,tp/BAR);
 
 		/* p/t scaling */
 		if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
@@ -1195,9 +1313,26 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		}
 		if(m) {
 			if(!(i%m)) {
-				p=get_total_p(moldyn);
+				momentum=get_total_p(moldyn);
 				dprintf(moldyn->mfd,
-				        "%f %f\n",moldyn->time,v3_norm(&p));
+				        "%f %f %f %f %f\n",moldyn->time,
+				        momentum.x,momentum.y,momentum.z,
+				        v3_norm(&momentum));
+			}
+		}
+		if(p) {
+			if(!(i%p)) {
+				dprintf(moldyn->pfd,
+				        "%f %f %f %f %f\n",moldyn->time,
+				         moldyn->p/BAR,moldyn->mean_p/BAR,
+				         moldyn->gp/BAR,moldyn->mean_gp/BAR);
+			}
+		}
+		if(t) {
+			if(!(i%t)) {
+				dprintf(moldyn->tfd,
+				        "%f %f %f\n",
+				        moldyn->time,moldyn->t,moldyn->mean_t);
 			}
 		}
 		if(s) {
@@ -1218,15 +1353,23 @@ int moldyn_integrate(t_moldyn *moldyn) {
 			if(!(i%v)) {
 				visual_atoms(&(moldyn->vis),moldyn->time,
 				             moldyn->atom,moldyn->count);
-				printf("\rsched: %d, steps: %d, T: %f, P: %f V: %f",
-				       sched->count,i,
-				       moldyn->t,moldyn->p/ATM,moldyn->volume);
-				fflush(stdout);
 			}
 		}
 
+		/* display progress */
+		if(!(i%10)) {
+			printf("\rsched: %d, steps: %d, T: %f, P: %f %f V: %f",
+			       sched->count,i,
+			       moldyn->mean_t,
+			       moldyn->mean_p/BAR,
+			       moldyn->mean_gp/BAR,
+			       moldyn->volume);
+			fflush(stdout);
+		}
+
 		/* increase absolute time */
 		moldyn->time+=moldyn->tau;
+		moldyn->total_steps+=1;
 
 	}
 
@@ -1313,6 +1456,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 	/* reset energy */
 	moldyn->energy=0.0;
 
+	/* reset global virial */
+	memset(&(moldyn->virial),0,sizeof(t_virial));
+
 	/* reset force, site energy and virial of every atom */
 	for(i=0;i<count;i++) {
 
@@ -1333,7 +1479,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 	}
 
-	/* get energy,force and virial of every atom */
+	/* get energy, force and virial of every atom */
+
+	/* first (and only) loop over atoms i */
 	for(i=0;i<count;i++) {
 
 		/* single particle potential/force */
@@ -1353,6 +1501,45 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 		dnlc=lc->dnlc;
 
+		/* first loop over atoms j */
+		if(moldyn->func2b) {
+			for(j=0;j<27;j++) {
+
+				this=&(neighbour_i[j]);
+				list_reset_f(this);
+
+				if(this->start==NULL)
+					continue;
+
+				bc_ij=(j<dnlc)?0:1;
+
+				do {
+					jtom=this->current->data;
+
+					if(jtom==&(itom[i]))
+						continue;
+
+					if((jtom->attr&ATOM_ATTR_2BP)&
+					   (itom[i].attr&ATOM_ATTR_2BP)) {
+						moldyn->func2b(moldyn,
+						               &(itom[i]),
+						               jtom,
+					        	       bc_ij);
+					}
+				} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+
+			}
+		}
+
+		/* 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++) {
 
 			this=&(neighbour_i[j]);
@@ -1369,25 +1556,70 @@ int potential_force_calc(t_moldyn *moldyn) {
 				if(jtom==&(itom[i]))
 					continue;
 
-				if((jtom->attr&ATOM_ATTR_2BP)&
-				   (itom[i].attr&ATOM_ATTR_2BP)) {
-					moldyn->func2b(moldyn,
-					               &(itom[i]),
-					               jtom,
-					               bc_ij);
-				}
+				if(!(jtom->attr&ATOM_ATTR_3BP))
+					continue;
 
-				/* 3 body potential/force */
+				/* reset 3bp run */
+				moldyn->run3bp=1;
 
-				if(!(itom[i].attr&ATOM_ATTR_3BP)||
-				   !(jtom->attr&ATOM_ATTR_3BP))
+				if(moldyn->func3b_j1)
+					moldyn->func3b_j1(moldyn,
+					                  &(itom[i]),
+					                  jtom,
+					                  bc_ij);
+
+				/* in first j loop, 3bp run can be skipped */
+				if(!(moldyn->run3bp))
 					continue;
+			
+				/* first loop over atoms k */
+				if(moldyn->func3b_k1) {
+
+				for(k=0;k<27;k++) {
+
+					that=&(neighbour_i2[k]);
+					list_reset_f(that);
+					
+					if(that->start==NULL)
+						continue;
+
+					bc_ik=(k<dnlc)?0:1;
+
+					do {
+
+						ktom=that->current->data;
+
+						if(!(ktom->attr&ATOM_ATTR_3BP))
+							continue;
+
+						if(ktom==jtom)
+							continue;
+
+						if(ktom==&(itom[i]))
+							continue;
+
+						moldyn->func3b_k1(moldyn,
+						                  &(itom[i]),
+						                  jtom,
+						                  ktom,
+						                  bc_ik|bc_ij);
+
+					} while(list_next_f(that)!=\
+					        L_NO_NEXT_ELEMENT);
+
+				}
+
+				}
 
-				/* copy the neighbour lists */
-				memcpy(neighbour_i2,neighbour_i,
-				       27*sizeof(t_list));
+				if(moldyn->func3b_j2)
+					moldyn->func3b_j2(moldyn,
+					                  &(itom[i]),
+					                  jtom,
+					                  bc_ij);
+
+				/* second loop over atoms k */
+				if(moldyn->func3b_k2) {
 
-				/* get neighbours of i */
 				for(k=0;k<27;k++) {
 
 					that=&(neighbour_i2[k]);
@@ -1411,37 +1643,53 @@ int potential_force_calc(t_moldyn *moldyn) {
 						if(ktom==&(itom[i]))
 							continue;
 
-						moldyn->func3b(moldyn,
-						               &(itom[i]),
-						               jtom,
-						               ktom,
-						               bc_ik|bc_ij);
+						moldyn->func3b_k2(moldyn,
+						                  &(itom[i]),
+						                  jtom,
+						                  ktom,
+						                  bc_ik|bc_ij);
 
 					} while(list_next_f(that)!=\
 					        L_NO_NEXT_ELEMENT);
 
 				}
+				
+				}
 
 				/* 2bp post function */
-				if(moldyn->func2b_post) {
-					moldyn->func2b_post(moldyn,
-					                    &(itom[i]),
-					                    jtom,bc_ij);
+				if(moldyn->func3b_j3) {
+					moldyn->func3b_j3(moldyn,
+					                  &(itom[i]),
+					                  jtom,bc_ij);
 				}
 					
 			} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 		
 		}
+		
+#ifdef DEBUG
+	//printf("\n\n");
+#endif
+#ifdef VDEBUG
+	printf("\n\n");
+#endif
 
 	}
 
 #ifdef DEBUG
-printf("\n\n");
-#endif
-#ifdef VDEBUG
-printf("\n\n");
+	printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
 #endif
 
+	/* calculate global virial */
+	for(i=0;i<count;i++) {
+		moldyn->virial.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x;
+		moldyn->virial.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y;
+		moldyn->virial.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z;
+		moldyn->virial.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x;
+		moldyn->virial.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x;
+		moldyn->virial.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y;
+	}
+
 	return 0;
 }
 
@@ -1449,7 +1697,8 @@ printf("\n\n");
  * virial calculation
  */
 
-inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
+//inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
+int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
 
 	a->virial.xx+=f->x*d->x;
 	a->virial.yy+=f->y*d->y;
@@ -1465,7 +1714,8 @@ inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
  * periodic boundayr checking
  */
 
-inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
+//inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
+int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
 	
 	double x,y,z;
 	t_3dvec *dim;
@@ -1492,45 +1742,6 @@ inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
 	return 0;
 }
         
-
-/*
- * example potentials
- */
-
-/* harmonic oscillator potential and force */
-
-int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
-
-	t_ho_params *params;
-	t_3dvec force,distance;
-	double d,f;
-	double sc,equi_dist;
-
-	params=moldyn->pot2b_params;
-	sc=params->spring_constant;
-	equi_dist=params->equilibrium_distance;
-
-	if(ai<aj) return 0;
-
-	v3_sub(&distance,&(aj->r),&(ai->r));
-	
-	if(bc) check_per_bound(moldyn,&distance);
-	d=v3_norm(&distance);
-	if(d<=moldyn->cutoff) {
-		moldyn->energy+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
-		/* f = -grad E; grad r_ij = -1 1/r_ij distance */
-		f=sc*(1.0-equi_dist/d);
-		v3_scale(&force,&distance,f);
-		v3_add(&(ai->f),&(ai->f),&force);
-		virial_calc(ai,&force,&distance);
-		virial_calc(aj,&force,&distance); /* f and d signe switched */
-		v3_scale(&force,&distance,-f);
-		v3_add(&(aj->f),&(aj->f),&force);
-	}
-
-	return 0;
-}
-
 /*
  * debugging / critical check functions
  */