X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Falbe_fast.c;h=46102303d364ca038496753053a16a056c0403eb;hb=452a348351ae8a2255809918a05683d2510655ce;hp=45ebff9203abf296c9d3b6f5e4af454ff7049599;hpb=45f9190fc7489f3f9ff8f0d917945b1d5bc9d584;p=physik%2Fposic.git

diff --git a/potentials/albe_fast.c b/potentials/albe_fast.c
index 45ebff9..4610230 100644
--- a/potentials/albe_fast.c
+++ b/potentials/albe_fast.c
@@ -21,6 +21,10 @@
 #include <omp.h>
 #endif
 
+#ifdef PTHREAD
+#include <pthread.h>
+#endif
+
 #include "../moldyn.h"
 #include "../math/math.h"
 #include "albe.h"
@@ -36,6 +40,116 @@
 				a->virial.xz+=f->x*d->z; \
 				a->virial.yz+=f->y*d->z
 
+#if 0
+#ifdef PTHREADS
+void *k1_calc(void *ptr) {
+
+	/* albe 3 body potential function (first k loop) */
+	
+	t_albe_mult_params *params;
+	t_albe_exchange *exchange;
+	unsigned char brand_i;
+	double Rk,Sk,Sk2,gamma_i,c_i,d_i,h_i,ci2,di2,ci2di2;
+	t_atom *ai,*jtom,*ktom;
+	
+
+	if(kcount>ALBE_MAXN) {
+		printf("FATAL: neighbours = %d\n",kcount);
+		printf("  -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
+	}
+
+	/* ik constants */
+	if(brand_i==ktom->brand) {
+		Rk=params->R[brand_i];
+		Sk=params->S[brand_i];
+		Sk2=params->S2[brand_i];
+		/* albe needs i,k depending c,d,h and gamma values */
+		gamma_i=params->gamma[brand_i];
+		c_i=params->c[brand_i];
+		d_i=params->d[brand_i];
+		h_i=params->h[brand_i];
+		ci2=params->c2[brand_i];
+		di2=params->d2[brand_i];
+		ci2di2=params->c2d2[brand_i];
+	}
+	else {
+		Rk=params->Rmixed;
+		Sk=params->Smixed;
+		Sk2=params->S2mixed;
+		/* albe needs i,k depending c,d,h and gamma values */
+		gamma_i=params->gamma_m;
+		c_i=params->c_mixed;
+		d_i=params->d_mixed;
+		h_i=params->h_mixed;
+		ci2=params->c2_mixed;
+		di2=params->d2_mixed;
+		ci2di2=params->c2d2_m;
+	}
+
+	/* dist_ik, d_ik2 */
+	v3_sub(&dist_ik,&(ktom->r),&(ai->r));
+	if(bc_ik) 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>Sk2) {
+		kcount++;
+		continue;
+	}
+
+	/* d_ik */
+	d_ik=sqrt(d_ik2);
+
+	/* 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..
+	*/
+
+	h_cos=h_i+cos_theta; // + in albe formalism
+	d2_h_cos2=di2+(h_cos*h_cos);
+	frac=ci2/d2_h_cos2;
+	g=gamma_i*(1.0+ci2di2-frac);
+	dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
+
+	/* zeta sum += f_c_ik * g_ijk */
+	if(d_ik<=Rk) {
+		zeta_ij+=g;
+		f_c_ik=1.0;
+		df_c_ik=0.0;
+	}
+	else {
+		s_r=Sk-Rk;
+		arg=M_PI*(d_ik-Rk)/s_r;
+		f_c_ik=0.5+0.5*cos(arg);
+		df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
+		zeta_ij+=f_c_ik*g;
+	}
+
+	/* store even more data for second k loop */
+	exchange->g[kcount]=g;
+	exchange->dg[kcount]=dg;
+	exchange->d_ik[kcount]=d_ik;
+	exchange->cos_theta[kcount]=cos_theta;
+	exchange->f_c_ik[kcount]=f_c_ik;
+	exchange->df_c_ik[kcount]=df_c_ik;
+
+	/* increase k counter */
+	kcount++;
+
+}
+#endif
+#endif
+
 int albe_potential_force_calc(t_moldyn *moldyn) {
 
 	int i,j,k,count;
@@ -275,6 +389,14 @@ int albe_potential_force_calc(t_moldyn *moldyn) {
 						if(ktom==&(itom[i]))
 							continue;
 
+#if 0
+//#ifdef PTHREADS
+	ret=pthread_create(&(k_thread[k]),NULL,k1_calc,k_data);
+	if(ret) {
+		perror("[albe fast] create thread\n");
+		return ret;
+	}
+#else
 
 /* k1 func here ... */
 /* albe 3 body potential function (first k loop) */
@@ -372,6 +494,8 @@ int albe_potential_force_calc(t_moldyn *moldyn) {
 	/* increase k counter */
 	kcount++;
 
+#endif // PTHREADS
+
 #ifdef STATIC_LISTS
 					}
 #elif LOWMEM_LISTS