From 452a348351ae8a2255809918a05683d2510655ce Mon Sep 17 00:00:00 2001 From: hackbard Date: Fri, 26 Sep 2008 14:10:58 +0200 Subject: [PATCH] pthread imp started for orig albe (more easier in the beginning) --- moldyn.c | 27 +++++++++ moldyn.h | 4 ++ potentials/albe.c | 21 +++++++ potentials/albe.h | 12 ++++ potentials/albe_fast.c | 124 +++++++++++++++++++++++++++++++++++++++++ 5 files changed, 188 insertions(+) diff --git a/moldyn.c b/moldyn.c index e936b11..8d50a69 100644 --- 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 --- 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 --- 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 --- 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 --- a/potentials/albe_fast.c +++ b/potentials/albe_fast.c @@ -21,6 +21,10 @@ #include #endif +#ifdef PTHREAD +#include +#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 -- 2.20.1