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

diff --git a/moldyn.c b/moldyn.c
index e936b11..8d50a69 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -1951,6 +1951,11 @@ int potential_force_calc(t_moldyn *moldyn) {
 #endif
        u8 bc_ij,bc_ik;
        int dnlc;
+#ifdef PTHREADS
+       int ret;
+       pthread_t kthread[27];
+       t_kdata kdata[27];
+#endif
 
        count=moldyn->count;
        itom=moldyn->atom;
@@ -1959,6 +1964,10 @@ int potential_force_calc(t_moldyn *moldyn) {
        atom=moldyn->atom;
 #endif
 
+#ifdef PTHREADS
+       memset(kdata,0,27*sizeof(t_kdata));
+#endif
+
        /* reset energy */
        moldyn->energy=0.0;
 
@@ -2127,7 +2136,9 @@ int potential_force_calc(t_moldyn *moldyn) {
                                        continue;
                        
                                /* first loop over atoms k */
+#ifndef PTHREADS
                                if(moldyn->func3b_k1) {
+#endif
 
                                for(k=0;k<27;k++) {
 
@@ -2166,11 +2177,25 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                if(ktom==&(itom[i]))
                                                        continue;
 
+#ifdef PTHREADS
+                                               kdata[k].moldyn=moldyn;
+                                               kdata[k].ai=&(itom[i]);
+                                               kdata[k].aj=jtom;
+                                               kdata[k].ak=ktom;
+                                               kdata[k].bc=bc_ik;
+       ret=pthread_create(&(kthread[k]),NULL,moldyn->func3b_k1,&(kdata[k]));
+       if(ret) {
+               perror("[moldyn] create k1 thread");
+               return ret;
+       }
+#else
                                                moldyn->func3b_k1(moldyn,
                                                                  &(itom[i]),
                                                                  jtom,
                                                                  ktom,
                                                                  bc_ik|bc_ij);
+#endif
+
 #ifdef STATIC_LISTS
                                        }
 #elif LOWMEM_LISTS
@@ -2182,7 +2207,9 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                }
 
+#ifndef PTHREADS
                                }
+#endif
 
                                if(moldyn->func3b_j2)
                                        moldyn->func3b_j2(moldyn,

diff --git a/moldyn.h b/moldyn.h
index 31794c8..6fda9ce 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -117,8 +117,12 @@ typedef struct s_moldyn {
        int (*func3b_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
        int (*func3b_j2)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
        int (*func3b_j3)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+#ifdef PTHREADS
+       void *(*func3b_k1)(void *ptr);
+#else
        int (*func3b_k1)(struct s_moldyn *moldyn,
                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bck);
+#endif
        int (*func3b_k2)(struct s_moldyn *moldyn,
                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bck);
        void *pot_params;

diff --git a/potentials/albe.c b/potentials/albe.c
index 08ce4fa..b00e173 100644 (file)
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -176,8 +176,12 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 }
 
 /* albe 3 body potential function (first k loop) */
+#ifdef PTHREADS
+void *albe_mult_3bp_k1(void *ptr) {
+#else
 int albe_mult_3bp_k1(t_moldyn *moldyn,
                      t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
+#endif
 
        t_albe_mult_params *params;
        t_albe_exchange *exchange;
@@ -188,6 +192,19 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
        double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
        double f_c_ik,df_c_ik;
        int kcount;
+#ifdef PTHREADS
+       t_kdata *kdata;
+       t_moldyn *moldyn;
+       t_atom *ai,*aj,*ak;
+       u8 bc;
+
+       kdata=ptr;
+       moldyn=kdata->moldyn;
+       ai=kdata->ai;
+       aj=kdata->aj;
+       ak=kdata->ak;
+       bc=kdata->bc;
+#endif
 
        params=moldyn->pot_params;
        exchange=&(params->exchange);
@@ -281,8 +298,12 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
        /* increase k counter */
        exchange->kcount++;
 
+#ifdef PTHREADS
+}
+#else
        return 0;
 }
+#endif
 
 int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 

diff --git a/potentials/albe.h b/potentials/albe.h
index 87f45a0..3d0d7e5 100644 (file)
--- a/potentials/albe.h
+++ b/potentials/albe.h
@@ -80,11 +80,23 @@ typedef struct s_albe_mult_params {
        t_albe_exchange exchange;       /* exchange between 2bp and 3bp calc */
 } t_albe_mult_params;
 
+#ifdef PTHREADS
+typedef struct s_kdata {
+       t_moldyn *moldyn;
+       t_atom *ai,*aj,*ak;
+       unsigned char bc;
+} t_kdata;
+#endif
+
 /* function prototypes */
 int albe_mult_set_params(t_moldyn *moldyn,int element1,int elemnt2);
 int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+#ifdef PTHREADS
+void *albe_mult_3bp_k1(void *ptr);
+#else
 int albe_mult_3bp_k1(t_moldyn *moldyn,
                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
+#endif
 int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
 int albe_mult_3bp_k2(t_moldyn *moldyn,
                         t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);

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