From: hackbard <hackbard@sage.physik.uni-augsburg.de>
Date: Sat, 13 Jun 2009 13:54:53 +0000 (+0200)
Subject: removed tersoff 1bp function, added PDEBUG define, check per bound with boundary... 
X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=commitdiff_plain;h=0d2f9a11030dff3583104dac5d4dcb9f040a1327

removed tersoff 1bp function, added PDEBUG define, check per bound with boundary conditions (for diffusion), adapted diffusion calc, modified DEBUG output in albe and tersoff pots, more mods to tersoff pot calc,
---

diff --git a/moldyn.c b/moldyn.c
index 7698c24..12cb49b 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -258,7 +258,7 @@ 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;
@@ -802,6 +802,8 @@ int del_atom(t_moldyn *moldyn,int tag) {
 			case DEFECT_STYPE_DB_Z:\
 				d_o.z=d_params->od;\
 				d_d.z=d_params->dd;\
+d_d.x=0.9;\
+d_d.y=0.9;\
 				break;\
 			case DEFECT_STYPE_DB_R:\
 				break;\
@@ -2071,11 +2073,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);
@@ -2182,15 +2184,20 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		}
 
 		/* display progress */
+#ifndef PDEBUG
 		if(!(i%10)) {
+#endif
 			/* get current time */
 			gettimeofday(&t2,NULL);
 
 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));
 
@@ -2198,7 +2205,9 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
 
 			/* copy over time */
 			t1=t2;
+#ifndef PDEBUG
 		}
+#endif
 
 		/* increase absolute time */
 		moldyn->time+=moldyn->tau;
@@ -2246,6 +2255,7 @@ int velocity_verlet(t_moldyn *moldyn) {
 		v3_add(&(atom[i].r),&(atom[i].r),&delta);
 		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)] */
@@ -2261,7 +2271,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
@@ -2713,6 +2727,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
  */
@@ -3048,6 +3107,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;
@@ -3061,8 +3121,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) {
diff --git a/moldyn.h b/moldyn.h
index 538e0eb..6310416 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -45,6 +45,7 @@ typedef struct s_atom {
 	u8 brand;		/* brand id */
 	int tag;		/* atom unique id (number of atom) */
 	u8 attr;		/* attributes */
+	int pbc[3];		/* pb crossing in x, y and z direction */
 } t_atom;
 
 #define ATOM_ATTR_FP	0x01	/* fixed position (bulk material) */
@@ -290,6 +291,11 @@ typedef struct s_defect_params {
 #define DEFECT_STYPE_DB_Z	3
 #define DEFECT_STYPE_DB_R	4
 
+typedef struct s_offset_params {
+	t_3dvec o;
+	u8 use;
+} t_offset_params;
+
 /*
  *
  *  defines
@@ -496,6 +502,7 @@ 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)
 //	__attribute__((always_inline));
 int check_per_bound(t_moldyn *moldyn,t_3dvec *a);
+int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a);
 //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a)
 //	__attribute__((always_inline));
 
diff --git a/potentials/albe.c b/potentials/albe.c
index 08ce4fa..1ec3938 100644
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -270,6 +270,11 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
 		exchange->zeta_ij+=f_c_ik*g;
 	}
 
+#ifdef DEBUG
+	if(ai==&(moldyn->atom[DATOM])) 
+		printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik);
+#endif
+
 	/* store even more data for second k loop */
 	exchange->g[kcount]=g;
 	exchange->dg[kcount]=dg;
@@ -367,19 +372,17 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	virial_calc(ai,&force,&(exchange->dist_ij));
 
 #ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
 	if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
 		printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
 		printf("  adding %f %f %f\n",force.x,force.y,force.z);
-		if(ai==&(moldyn->atom[0]))
+		if(ai==&(moldyn->atom[DATOM]))
 			printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
-		if(aj==&(moldyn->atom[0]))
+		if(aj==&(moldyn->atom[DATOM]))
 			printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 		printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 		                                    f_c,b,f_a,f_r);
 		printf("          %f %f %f\n",exchange->zeta_ij,.0,.0);
 	}
-}
 #endif
 
 	/* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
@@ -477,7 +480,6 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
 	v3_add(&(aj->f),&(aj->f),&force);
 
 #ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
 	if(aj==&(moldyn->atom[DATOM])) {
 		printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 		printf("  adding %f %f %f\n",force.x,force.y,force.z);
@@ -485,7 +487,6 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 		printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 		printf("    d ij ik = %f %f\n",d_ij,d_ik);
 	}
-}
 #endif
 
 	/* virial */
@@ -505,7 +506,6 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 	v3_add(&(ak->f),&(ak->f),&force);
 
 #ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
 	if(ak==&(moldyn->atom[DATOM])) {
 		printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
 		printf("  adding %f %f %f\n",force.x,force.y,force.z);
@@ -513,7 +513,6 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 		printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 		printf("    d ij ik = %f %f\n",d_ij,d_ik);
 	}
-}
 #endif
 
 	/* virial */
diff --git a/potentials/tersoff.c b/potentials/tersoff.c
index 5039ff7..4ca7141 100644
--- a/potentials/tersoff.c
+++ b/potentials/tersoff.c
@@ -85,20 +85,28 @@ int tersoff_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 			p->h[1]=TM_H_C;
 			break;
 		default:
-			printf("[tersoff] WARNING: element1\n");
+			printf("[tersoff] WARNING: element2\n");
 			return -1;
 	}
 
 	printf("[tersoff] parameter completion\n");
 	p->S2[0]=p->S[0]*p->S[0];
 	p->S2[1]=p->S[1]*p->S[1];
-	p->Smixed=sqrt(p->S[0]*p->S[1]);
-	p->S2mixed=p->Smixed*p->Smixed;
+	p->S2mixed=p->S[0]*p->S[1];
+	p->Smixed=sqrt(p->S2mixed);
 	p->Rmixed=sqrt(p->R[0]*p->R[1]);
 	p->Amixed=sqrt(p->A[0]*p->A[1]);
 	p->Bmixed=sqrt(p->B[0]*p->B[1]);
 	p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
 	p->mu_m=0.5*(p->mu[0]+p->mu[1]);
+	p->betaini[0]=pow(p->beta[0],p->n[0]);
+	p->betaini[1]=pow(p->beta[1],p->n[1]);
+	p->ci2[0]=p->c[0]*p->c[0];
+	p->ci2[1]=p->c[1]*p->c[1];
+	p->di2[0]=p->d[0]*p->d[0];
+	p->di2[1]=p->d[1]*p->d[1];
+	p->ci2di2[0]=p->ci2[0]/p->di2[0];
+	p->ci2di2[1]=p->ci2[1]/p->di2[1];
 
 	printf("[tersoff] mult parameter info:\n");
 	printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
@@ -118,36 +126,6 @@ int tersoff_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 	return 0;
 }
 
-/* tersoff 1 body part */
-int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
-
-	int brand;
-	t_tersoff_mult_params *params;
-	t_tersoff_exchange *exchange;
-	
-	brand=ai->brand;
-	params=moldyn->pot_params;
-	exchange=&(params->exchange);
-
-	/*
-	 * simple: point constant parameters only depending on atom i to
-	 *         their right values
-	 */
-
-	exchange->beta_i=&(params->beta[brand]);
-	exchange->n_i=&(params->n[brand]);
-	exchange->c_i=&(params->c[brand]);
-	exchange->d_i=&(params->d[brand]);
-	exchange->h_i=&(params->h[brand]);
-
-	exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
-	exchange->ci2=params->c[brand]*params->c[brand];
-	exchange->di2=params->d[brand]*params->d[brand];
-	exchange->ci2di2=exchange->ci2/exchange->di2;
-
-	return 0;
-}
-	
 /* tersoff 2 body part */
 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
@@ -239,7 +217,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 #endif
 
 	/* virial */
-	virial_calc(aj,&force,&dist_ij);
+	virial_calc(ai,&force,&dist_ij);
 
 	/* energy 2bp contribution */
 	energy=f_r*f_c;
@@ -364,10 +342,10 @@ int tersoff_mult_3bp_k1(t_moldyn *moldyn,
 	cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 
 	/* g_ijk */
-	h_cos=*(exchange->h_i)-cos_theta;
-	d2_h_cos2=exchange->di2+(h_cos*h_cos);
-	frac=exchange->ci2/d2_h_cos2;
-	g=1.0+exchange->ci2di2-frac;
+	h_cos=params->h[brand]-cos_theta;
+	d2_h_cos2=params->di2[brand]+(h_cos*h_cos);
+	frac=params->ci2[brand]/d2_h_cos2;
+	g=1.0+params->ci2di2[brand]-frac;
 	dg=-2.0*frac*h_cos/d2_h_cos2;
 
 	/* zeta sum += f_c_ik * g_ijk */
@@ -385,11 +363,8 @@ int tersoff_mult_3bp_k1(t_moldyn *moldyn,
 	}
 
 #ifdef DEBUG
-	if((ai==&(moldyn->atom[0]))|
-	   (aj==&(moldyn->atom[864]))|
-	   (ak==&(moldyn->atom[1003]))) {
-		printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos);
-	}
+	if(ai==&(moldyn->atom[DATOM]))
+		printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik);
 #endif
 
 	/* store even more data for second k loop */
@@ -425,8 +400,8 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 
-	brand=aj->brand;
-	if(brand==ai->brand) {
+	brand=ai->brand;
+	if(brand==aj->brand) {
 		S=params->S[brand];
 		R=params->R[brand];
 		B=params->B[brand];
@@ -473,8 +448,8 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		db=0.0;
 	}
 	else {
-		ni=*(exchange->n_i);
-		tmp=exchange->betaini*pow(exchange->zeta_ij,ni-1.0);
+		ni=params->n[brand];
+		tmp=params->betaini[brand]*pow(exchange->zeta_ij,ni-1.0);
 		b=(1.0+exchange->zeta_ij*tmp);
 		db=chi*pow(b,-1.0/(2.0*ni)-1.0);
 		b=db*b;
@@ -489,12 +464,12 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	v3_add(&(aj->f),&(aj->f),&force);
 
 #ifdef DEBUG
-	if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
+	if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
 		printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
 		printf("adding %f %f %f\n",force.x,force.y,force.z);
-		if(ai==&(moldyn->atom[0]))
+		if(ai==&(moldyn->atom[DATOM]))
 			printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
-		if(aj==&(moldyn->atom[0]))
+		if(aj==&(moldyn->atom[DATOM]))
 			printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
 		printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 		                                    f_c,b,f_a,f_r);
@@ -503,7 +478,7 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 #endif
 
 	/* virial */
-	virial_calc(aj,&force,&(exchange->dist_ij));
+	virial_calc(ai,&force,&(exchange->dist_ij));
 
 	/* dzeta prefactor = - 0.5 f_c f_a db */
 	exchange->pre_dzeta=-0.5*f_a*f_c*db;
@@ -534,7 +509,7 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
 	double pre_dzeta;
 	double f_c_ik,df_c_ik;
 	double dijdik_inv,fcdg,dfcg;
-	t_3dvec dcosdri,dcosdrj,dcosdrk;
+	t_3dvec dcosdrj,dcosdrk;
 	t_3dvec force,tmp;
 
 	params=moldyn->pot_params;
@@ -588,30 +563,11 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
 	v3_scale(&dcosdrk,&dist_ij,dijdik_inv);
 	v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 	v3_add(&dcosdrk,&dcosdrk,&tmp);
-	v3_add(&dcosdri,&dcosdrj,&dcosdrk);
-	v3_scale(&dcosdri,&dcosdri,-1.0);
 
 	/* f_c_ik * dg, df_c_ik * g */
 	fcdg=f_c_ik*dg;
 	dfcg=df_c_ik*g;
 
-	/* derivative wrt i */
-	v3_scale(&force,&dist_ik,dfcg);
-	v3_scale(&tmp,&dcosdri,fcdg);
-	v3_add(&force,&force,&tmp);
-	v3_scale(&force,&force,pre_dzeta);
-
-	/* force contribution */
-	v3_add(&(ai->f),&(ai->f),&force);
-	
-#ifdef DEBUG
-	if(ai==&(moldyn->atom[0])) {
-		printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
-		printf("adding %f %f %f\n",force.x,force.y,force.z);
-		printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
-	}
-#endif
-
 	/* derivative wrt j */
 	v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 
@@ -619,17 +575,22 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
 	v3_add(&(aj->f),&(aj->f),&force);
 
 #ifdef DEBUG
-	if(aj==&(moldyn->atom[0])) {
+	if(aj==&(moldyn->atom[DATOM])) {
 		printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
-		printf("adding %f %f %f\n",force.x,force.y,force.z);
-		printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+		printf("  adding %f %f %f\n",force.x,force.y,force.z);
+		printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+		printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
+		printf("    d ij ik = %f %f\n",d_ij,d_ik);
 	}
 #endif
 
 	/* virial */
-	v3_scale(&force,&force,-1.0);
 	virial_calc(ai,&force,&dist_ij);
 
+	/* force contribution to atom i */
+	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
+
 	/* derivative wrt k */
 	v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
 	v3_scale(&tmp,&dcosdrk,fcdg);
@@ -640,19 +601,24 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
 	v3_add(&(ak->f),&(ak->f),&force);
 
 #ifdef DEBUG
-	if(ak==&(moldyn->atom[0])) {
+	if(ak==&(moldyn->atom[DATOM])) {
 		printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
-		printf("adding %f %f %f\n",force.x,force.y,force.z);
-		printf("total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
+		printf("  adding %f %f %f\n",force.x,force.y,force.z);
+		printf("  total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
+		printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
+		printf("    d ij ik = %f %f\n",d_ij,d_ik);
 	}
 #endif
 
 	/* virial */
-	v3_scale(&force,&force,-1.0);
 	virial_calc(ai,&force,&dist_ik);
+
+	/* force contribution to atom i */
+	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
 	
 	/* increase k counter */
-	exchange->kcount++;	
+	exchange->kcount++;
 
 	return 0;
 
diff --git a/potentials/tersoff.h b/potentials/tersoff.h
index a23503a..2e4f242 100644
--- a/potentials/tersoff.h
+++ b/potentials/tersoff.h
@@ -28,17 +28,6 @@ typedef struct s_tersoff_echange {
 	double dg[TERSOFF_MAXN];
 	double cos_theta[TERSOFF_MAXN];
 
-	double *beta_i;
-	double *n_i;
-	double *c_i;
-	double *d_i;
-	double *h_i;
-
-	double ci2;
-	double di2;
-	double ci2di2;
-	double betaini;
-
 	double zeta_ij;
 	double pre_dzeta;
 
@@ -70,12 +59,17 @@ typedef struct s_tersoff_mult_params {
 	double d[2];
 	double h[2];
 
+	double ci2[2];
+	double di2[2];
+	double ci2di2[2];
+	double betaini[2];
+
 	t_tersoff_exchange exchange;	/* exchange between 2bp and 3bp calc */
 } t_tersoff_mult_params;
 
 /* function prototypes */
 int tersoff_mult_set_params(t_moldyn *moldyn,int element1,int element2);
-int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai);
+//int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai);
 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
 int tersoff_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
 int tersoff_mult_3bp_k1(t_moldyn *moldyn,