X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=potentials%2Falbe.c;h=9f9828280acc5e0c908e44d0aaa0ce4dc5ba6222;hp=1ec3938844aa80bcb0043be78133c0ccfd08751a;hb=6e6d7126ea9a845f11637d8e1b8eb2b570ac4dc9;hpb=98df942a2412e6ace72ae221f2e967f24bba7108

diff --git a/potentials/albe.c b/potentials/albe.c
index 1ec3938..9f98282 100644
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -53,7 +53,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 			p->mu[0]=ALBE_MU_SI;
 			p->gamma[0]=ALBE_GAMMA_SI;
 			p->c[0]=ALBE_C_SI;
+			p->c2[0]=p->c[0]*p->c[0];
 			p->d[0]=ALBE_D_SI;
+			p->d2[0]=p->d[0]*p->d[0];
+			p->c2d2[0]=p->c2[0]/p->d2[0];
 			p->h[0]=ALBE_H_SI;
 			switch(element2) {
 				case C:
@@ -67,7 +70,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 					p->mu[1]=ALBE_MU_C;
 					p->gamma[1]=ALBE_GAMMA_C;
 					p->c[1]=ALBE_C_C;
+					p->c2[1]=p->c[1]*p->c[1];
 					p->d[1]=ALBE_D_C;
+					p->d2[1]=p->d[1]*p->d[1];
+					p->c2d2[1]=p->c2[1]/p->d2[1];
 					p->h[1]=ALBE_H_C;
 					/* mixed type: silicon carbide */
 					p->Smixed=ALBE_S_SIC;
@@ -79,7 +85,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 					p->mu_m=ALBE_MU_SIC;
 					p->gamma_m=ALBE_GAMMA_SIC;
 					p->c_mixed=ALBE_C_SIC;
+					p->c2_mixed=p->c_mixed*p->c_mixed;
 					p->d_mixed=ALBE_D_SIC;
+					p->d2_mixed=p->d_mixed*p->d_mixed;
+					p->c2d2_m=p->c2_mixed/p->d2_mixed;
 					p->h_mixed=ALBE_H_SIC;
 					break;
 				default:
@@ -114,34 +123,56 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
 	printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
 	                                  p->lambda_m);
 	printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
-	printf("  gamma  | %f | %f\n",p->gamma[0],p->gamma[1]);
-	printf("  c      | %f | %f\n",p->c[0],p->c[1]);
-	printf("  d      | %f | %f\n",p->d[0],p->d[1]);
-	printf("  h      | %f | %f\n",p->h[0],p->h[1]);
+	printf("  gamma  | %f | %f | %f\n",p->gamma[0],p->gamma[1],p->gamma_m);
+	printf("  c      | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed);
+	printf("  d      | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed);
+	printf("  c2     | %f | %f | %f\n",p->c2[0],p->c2[1],p->c2_mixed);
+	printf("  d2     | %f | %f | %f\n",p->d2[0],p->d2[1],p->d2_mixed);
+	printf("  c2d2   | %f | %f | %f\n",p->c2d2[0],p->c2d2[1],p->c2d2_m);
+	printf("  h      | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed);
 
 	return 0;
 }
 
-/* albe 3 body potential function (first ij loop) */
-int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+/* first i loop */
+int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) {
 
 	t_albe_mult_params *params;
 	t_albe_exchange *exchange;
-	unsigned char brand;
-	double S2;
-	t_3dvec dist_ij;
-	double d_ij2,d_ij;
+	
+	int i;
 
 	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 
-	/* reset zeta sum */
-	exchange->zeta_ij=0.0;
+	/* zero exchange values */
+	memset(exchange->zeta,0,ALBE_MAXN*sizeof(double));
+	for(i=0;i<ALBE_MAXN;i++)
+		memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(t_3dvec));
+	exchange->jcnt=0;
+	exchange->j2cnt=0;
 
-	/*
-	 * set ij depending values
-	 */
+	return 0;
+}
 
+/* first j loop within first i loop */
+int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+	t_albe_mult_params *params;
+	t_albe_exchange *exchange;
+
+	double S2,S,R,d2,d,s_r,arg;
+	t_3dvec dist;
+	int j;
+	u8 brand;
+	
+	params=moldyn->pot_params;
+	exchange=&(params->exchange);
+
+	/* get j counter */
+	j=exchange->jcnt;
+
+	/* set ij depending values */
 	brand=ai->brand;
 	if(brand==aj->brand) {
 		S2=params->S2[brand];
@@ -151,123 +182,132 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	}
 
 	/* dist_ij, d_ij2 */
-	v3_sub(&dist_ij,&(aj->r),&(ai->r));
-	if(bc) check_per_bound(moldyn,&dist_ij);
-	d_ij2=v3_absolute_square(&dist_ij);
+	v3_sub(&dist,&(aj->r),&(ai->r));
+	if(bc) check_per_bound(moldyn,&dist);
+	exchange->dist[j]=dist;
+	d2=v3_absolute_square(&dist);
+	exchange->d2[j]=d2;
 
 	/* if d_ij2 > S2 => no force & potential energy contribution */
-	if(d_ij2>S2) {
+	if(d2>S2) {
 		moldyn->run3bp=0;
+		exchange->skip[j]=1;
+		exchange->jcnt+=1;
 		return 0;
 	}
+	exchange->skip[j]=0;
+
+	/* more ij depending values */
+	if(brand==aj->brand) {
+		R=params->R[brand];
+		S=params->S[brand];
+		/* albe needs i,(j/k) depending c,d,h and gamma values */
+		exchange->gamma_[j]=&(params->gamma[brand]);
+		exchange->c_[j]=&(params->c[brand]);
+		exchange->d_[j]=&(params->d[brand]);
+		exchange->h_[j]=&(params->h[brand]);
+		exchange->c2_[j]=&(params->c2[brand]);
+		exchange->d2_[j]=&(params->d2[brand]);
+		exchange->c2d2_[j]=&(params->c2d2[brand]);
+	}
+	else {
+		R=params->Rmixed;
+		S=params->Smixed;
+		/* albe needs i,(j/k) depending c,d,h and gamma values */
+		exchange->gamma_[j]=&(params->gamma_m);
+		exchange->c_[j]=&(params->c_mixed);
+		exchange->d_[j]=&(params->d_mixed);
+		exchange->h_[j]=&(params->h_mixed);
+		exchange->c2_[j]=&(params->c2_mixed);
+		exchange->d2_[j]=&(params->d2_mixed);
+		exchange->c2d2_[j]=&(params->c2d2_m);
+	}
 
 	/* d_ij */
-	d_ij=sqrt(d_ij2);
+	d=sqrt(d2);
+	exchange->d[j]=d;
+	
+	/* f_c, df_c */
+	if(d<R) {
+		exchange->f_c[j]=1.0;
+		exchange->df_c[j]=0.0;
+	}
+	else {
+		s_r=S-R;
+		arg=M_PI*(d-R)/s_r;
+		exchange->f_c[j]=0.5+0.5*cos(arg);
+		exchange->df_c[j]=0.5*sin(arg)*(M_PI/(s_r*d));
+	}
 
-	/* store values */
-	exchange->dist_ij=dist_ij;
-	exchange->d_ij2=d_ij2;
-	exchange->d_ij=d_ij;
+	/* reset k counter */
+	exchange->kcnt=0;
 
-	/* reset k counter for first k loop */
-	exchange->kcount=0;
-		
 	return 0;
 }
 
-/* albe 3 body potential function (first k loop) */
-int albe_mult_3bp_k1(t_moldyn *moldyn,
-                     t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+/* first k loop within first j loop within first i loop */
+int albe_mult_i0_j0_k0(t_moldyn *moldyn,
+                       t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
 	t_albe_mult_params *params;
 	t_albe_exchange *exchange;
-	unsigned char brand;
-	double R,S,S2;
-	t_3dvec dist_ij,dist_ik;
-	double d_ik2,d_ik,d_ij;
-	double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
-	double f_c_ik,df_c_ik;
-	int kcount;
+
+	int j,k;
+	t_3dvec distj,distk;
+	double dj,dk,djdk_inv,cos_theta;
+	double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j;
+	double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k;
+	t_3dvec dcosdrj,dcosdrk,tmp;
+	t_3dvec *dzjj,*dzkk,*dzjk,*dzkj;
 
 	params=moldyn->pot_params;
 	exchange=&(params->exchange);
-	kcount=exchange->kcount;
 
-	if(kcount>ALBE_MAXN) {
-		printf("FATAL: neighbours = %d\n",kcount);
-		printf("  -> %d %d %d\n",ai->tag,aj->tag,ak->tag);
+	if(aj==ak) {
+		exchange->kcnt+=1;
+		return 0;
 	}
 
-	/* ik constants */
-	brand=ai->brand;
-	if(brand==ak->brand) {
-		R=params->R[brand];
-		S=params->S[brand];
-		S2=params->S2[brand];
-		/* albe needs i,k depending c,d,h and gamma values */
-		exchange->gamma_i=&(params->gamma[brand]);
-		exchange->c_i=&(params->c[brand]);
-		exchange->d_i=&(params->d[brand]);
-		exchange->h_i=&(params->h[brand]);
-	}
-	else {
-		R=params->Rmixed;
-		S=params->Smixed;
-		S2=params->S2mixed;
-		/* albe needs i,k depending c,d,h and gamma values */
-		exchange->gamma_i=&(params->gamma_m);
-		exchange->c_i=&(params->c_mixed);
-		exchange->d_i=&(params->d_mixed);
-		exchange->h_i=&(params->h_mixed);
+	/* k<j & check whether to run k */
+	j=exchange->jcnt;
+	k=exchange->kcnt;
+	if(k>=ALBE_MAXN) {
+		printf("FATAL: too many neighbours! (%d)\n",k);
+		printf("  atom i:%d | j:%d | k:%d\n",ai->tag,aj->tag,ak->tag);
 	}
-	exchange->ci2=*(exchange->c_i)**(exchange->c_i);
-	exchange->di2=*(exchange->d_i)**(exchange->d_i);
-	exchange->ci2di2=exchange->ci2/exchange->di2;
-
-	/* dist_ik, d_ik2 */
-	v3_sub(&dist_ik,&(ak->r),&(ai->r));
-	if(bc) check_per_bound(moldyn,&dist_ik);
-	d_ik2=v3_absolute_square(&dist_ik);
-
-	/* store data for second k loop */
-	exchange->dist_ik[kcount]=dist_ik;
-	exchange->d_ik2[kcount]=d_ik2;
-
-	/* return if not within cutoff */
-	if(d_ik2>S2) {
-		exchange->kcount++;
+	if((k>=j)|(exchange->skip[k])) {
+		exchange->kcnt+=1;
 		return 0;
 	}
 
-	/* d_ik */
-	d_ik=sqrt(d_ik2);
-
-	/* dist_ij, d_ij */
-	dist_ij=exchange->dist_ij;
-	d_ij=exchange->d_ij;
+	/* distances */
+	distj=exchange->dist[j];
+	distk=exchange->dist[k];
+	dj=exchange->d[j];
+	dk=exchange->d[k];
+	djdk_inv=1.0/(dj*dk);
 
 	/* cos theta */
-	cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
-
-	/* g_ijk */
-	h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
-	d2_h_cos2=exchange->di2+(h_cos*h_cos);
-	frac=exchange->ci2/d2_h_cos2;
-	g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
-	dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
-
-	/* zeta sum += f_c_ik * g_ijk */
-	if(d_ik<=R) {
-		exchange->zeta_ij+=g;
-		f_c_ik=1.0;
-		df_c_ik=0.0;
+	cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv;
+
+	/* g(cos(theta)) ij and ik values */
+	h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism
+	d2_h_cos2_j=*exchange->d2_[j]+(h_cos_j*h_cos_j);
+	frac_j=*exchange->c2_[j]/d2_h_cos2_j;
+	gj=1.0+*exchange->c2d2_[j]-frac_j;
+	gj*=*(exchange->gamma_[j]);
+	dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe 
+	if(ak->brand==aj->brand) {
+		gk=gj;
+		dgk=dgj;
 	}
 	else {
-		s_r=S-R;
-		arg=M_PI*(d_ik-R)/s_r;
-		f_c_ik=0.5+0.5*cos(arg);
-		df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
-		exchange->zeta_ij+=f_c_ik*g;
+		h_cos_k=*(exchange->h_[k])+cos_theta;
+		d2_h_cos2_k=*exchange->d2_[k]+(h_cos_k*h_cos_k);
+		frac_k=*exchange->c2_[k]/d2_h_cos2_k;
+		gk=1.0+*exchange->c2d2_[k]-frac_k;
+		gk*=*(exchange->gamma_[k]);
+		dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k;
 	}
 
 #ifdef DEBUG
@@ -284,33 +324,62 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
 	exchange->df_c_ik[kcount]=df_c_ik;
 
 	/* increase k counter */
-	exchange->kcount++;
+	exchange->kcnt+=1;
+		
+	return 0;
+}
+
+/* first j loop within first i loop */
+int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+
+	t_albe_mult_params *params;
+	t_albe_exchange *exchange;
+
+	params=moldyn->pot_params;
+	exchange=&(params->exchange);
+
+	/* increase j counter */
+	exchange->jcnt+=1;
 
 	return 0;
 }
 
-int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
+/* second j loop within first i loop */
+int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
 	t_albe_mult_params *params;
 	t_albe_exchange *exchange;
-	t_3dvec force;
-	double f_a,df_a,b,db,f_c,df_c;
-	double f_r,df_r;
-	double scale;
-	double mu,B;
-	double lambda,A;
-	double d_ij,r0;
-	unsigned char brand;
-	double S,R,s_r,arg;
+
+	int j;
+	double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db;
+	double A,B,mu,lambda,r0;
 	double energy;
+	t_3dvec *dist,force;
+	double scale;
+	u8 brand;
 
 	params=moldyn->pot_params;
 	exchange=&(params->exchange);
 
+	/* get j counter */
+	j=exchange->j2cnt;
+
+	/* skip if j not within cutoff */
+	if(exchange->skip[j]) {
+		moldyn->run3bp=0;
+		exchange->j2cnt+=1;
+		return 0;
+	}
+
+	/* distance */
+	d=exchange->d[j];
+	dist=&(exchange->dist[j]);
+	f_c=exchange->f_c[j];
+	df_c=exchange->df_c[j];
+
+	/* determine parameters to calculate fa, dfa, fr, dfr */
 	brand=aj->brand;
 	if(brand==ai->brand) {
-		S=params->S[brand];
-		R=params->R[brand];
 		B=params->B[brand];
 		A=params->A[brand];
 		r0=params->r0[brand];
@@ -318,8 +387,6 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		lambda=params->lambda[brand];
 	}
 	else {
-		S=params->Smixed;
-		R=params->Rmixed;
 		B=params->Bmixed;
 		A=params->Amixed;
 		r0=params->r0_mixed;
@@ -327,44 +394,36 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		lambda=params->lambda_m;
 	}
 
-	d_ij=exchange->d_ij;
-
-	/* f_c, df_c */
-	if(d_ij<R) {
-		f_c=1.0;
-		df_c=0.0;
-	}
-	else {
-		s_r=S-R;
-		arg=M_PI*(d_ij-R)/s_r;
-		f_c=0.5+0.5*cos(arg);
-		df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
-	}
-
 	/* f_a, df_a */
-	f_a=-B*exp(-mu*(d_ij-r0));
-	df_a=mu*f_a/d_ij;
+	f_a=-B*exp(-mu*(d-r0));
+	df_a=mu*f_a/d;
 
 	/* f_r, df_r */
-	f_r=A*exp(-lambda*(d_ij-r0));
-	df_r=lambda*f_r/d_ij;
+	f_r=A*exp(-lambda*(d-r0));
+	df_r=lambda*f_r/d;
 
 	/* b, db */
-	if(exchange->zeta_ij==0.0) {
-		b=1.0;
-		db=0.0;
-	}
-	else {
-		b=1.0/sqrt(1.0+exchange->zeta_ij);
-		db=-0.5*b/(1.0+exchange->zeta_ij);
-	}
+	b=1.0/sqrt(1.0+exchange->zeta[j]);
+	db=-0.5*b/(1.0+exchange->zeta[j]);
+
+	/* energy contribution */
+	energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
+	moldyn->energy+=energy;
+	ai->e+=energy;
 
-	/* force contribution for atom i */
+	/* force contribution for atom i due to ij bond */
 	scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
-	v3_scale(&force,&(exchange->dist_ij),scale);
+	v3_scale(&force,dist,scale);
 	v3_add(&(ai->f),&(ai->f),&force);
 
-	/* force contribution for atom j */
+#ifdef NDEBUG
+if(ai->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (ij) %.15f\n",force.x,force.y,force.z,ai->tag,aj->tag,exchange->zeta[j]);
+printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
+}
+#endif
+
+	/* force contribution for atom j due to ij bond */
 	v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 	v3_add(&(aj->f),&(aj->f),&force);
 
@@ -385,57 +444,57 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	}
 #endif
 
+	/* virial */
+	virial_calc(ai,&force,dist);
+
 	/* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
 	exchange->pre_dzeta=0.5*f_a*f_c*db;
 
-	/* energy contribution */
-	energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
-	moldyn->energy+=energy;
-	ai->e+=energy;
+	/* force contribution (drj derivative) */
+	v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta);
+	v3_add(&(aj->f),&(aj->f),&force);
+
+#ifdef NDEBUG
+if(aj->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (j der)\n",force.x,force.y,force.z,aj->tag,ai->tag);
+printf("    t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z);
+}
+#endif
+
+	/* virial */
+	virial_calc(ai,&force,dist);
+
+	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
+
+#ifdef NDEBUG
+if(ai->tag==0) {
+printf("force: %.15f %.15f %.15f | %d %d (i contr j der)\n",force.x,force.y,force.z,ai->tag,aj->tag);
+printf("    t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z);
+}
+#endif
 
 	/* reset k counter for second k loop */
-	exchange->kcount=0;
+	exchange->kcnt=0;
 		
 	return 0;
 }
 
-/* albe 3 body potential function (second k loop) */
-int albe_mult_3bp_k2(t_moldyn *moldyn,
-                     t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+/* second k loop within second j loop within first i loop */
+int albe_mult_i0_j2_k0(t_moldyn *moldyn,
+                       t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
 	t_albe_mult_params *params;
 	t_albe_exchange *exchange;
-	int kcount;
-	t_3dvec dist_ik,dist_ij;
-	double d_ik2,d_ik,d_ij2,d_ij;
-	unsigned char brand;
-	double S2;
-	double g,dg,cos_theta;
-	double pre_dzeta;
-	double f_c_ik,df_c_ik;
-	double dijdik_inv,fcdg,dfcg;
-	t_3dvec dcosdrj,dcosdrk;
-	t_3dvec force,tmp;
+
+	int j,k;
+	t_3dvec force;
 
 	params=moldyn->pot_params;
 	exchange=&(params->exchange);
-	kcount=exchange->kcount;
-
-	if(kcount>ALBE_MAXN)
-		printf("FATAL: neighbours!\n");
-
-	/* d_ik2 */
-	d_ik2=exchange->d_ik2[kcount];
 
-	brand=ak->brand;
-	if(brand==ai->brand)
-		S2=params->S2[brand];
-	else
-		S2=params->S2mixed;
-
-	/* return if d_ik > S */
-	if(d_ik2>S2) {
-		exchange->kcount++;
+	if(aj==ak) {
+		exchange->kcnt+=1;
 		return 0;
 	}
 
@@ -502,8 +561,8 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
 	v3_add(&force,&force,&tmp);
 	v3_scale(&force,&force,pre_dzeta);
 
-	/* force contribution */
-	v3_add(&(ak->f),&(ak->f),&force);
+	v3_scale(&force,&force,-1.0);
+	v3_add(&(ai->f),&(ai->f),&force);
 
 #ifdef DEBUG
 	if(ak==&(moldyn->atom[DATOM])) {
@@ -523,10 +582,23 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
 	v3_add(&(ai->f),&(ai->f),&force);
 
 	/* increase k counter */
-	exchange->kcount++;	
+	exchange->kcnt+=1;
 
 	return 0;
+}
+
+int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
+	t_albe_mult_params *params;
+	t_albe_exchange *exchange;
+
+	params=moldyn->pot_params;
+	exchange=&(params->exchange);
+
+	/* increase j counter */
+	exchange->j2cnt+=1;
+
+	return 0;
 }
 
 int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) {