pthread imp started for orig albe (more easier in the beginning)

[physik/posic.git] / potentials / albe_fast.c

diff --git a/potentials/albe_fast.c b/potentials/albe_fast.c
index 45ebff9..4610230 100644 (file)
--- 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