X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Falbe_fast.c;h=46102303d364ca038496753053a16a056c0403eb;hb=452a348351ae8a2255809918a05683d2510655ce;hp=3b5f68100a7b0e50463e1a2d97fcd4dc96a5203e;hpb=cb2811d92cb79b1f24d11d4f4897426ce5f7dc3a;p=physik%2Fposic.git diff --git a/potentials/albe_fast.c b/potentials/albe_fast.c index 3b5f681..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; @@ -45,7 +159,9 @@ int albe_potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS int *neighbour_i[27]; int p,q; - t_atom *atom; +#elif LOWMEM_LISTS + int neighbour_i[27]; + int p,q; #else t_list neighbour_i[27]; t_list neighbour_i2[27]; @@ -90,13 +206,18 @@ int albe_potential_force_calc(t_moldyn *moldyn) { t_3dvec dcosdrj,dcosdrk; t_3dvec tmp; + // even more ... + double gamma_i; + double c_i; + double d_i; + double h_i; + double ci2; + double di2; + double ci2di2; count=moldyn->count; itom=moldyn->atom; lc=&(moldyn->lc); -#ifdef STATIC_LISTS - atom=moldyn->atom; -#endif // optimized params=moldyn->pot_params; @@ -149,7 +270,9 @@ int albe_potential_force_calc(t_moldyn *moldyn) { dnlc=lc->dnlc; /* copy the neighbour lists */ -#ifndef STATIC_LISTS +#ifdef STATIC_LISTS +#elif LOWMEM_LISTS +#else memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); #endif @@ -163,10 +286,17 @@ int albe_potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS p=0; - while(neighbour_i[j][p]!=0) { + while(neighbour_i[j][p]!=-1) { - jtom=&(atom[neighbour_i[j][p]]); + jtom=&(itom[neighbour_i[j][p]]); p++; +#elif LOWMEM_LISTS + p=neighbour_i[j]; + + while(p!=-1) { + + jtom=&(itom[p]); + p=lc->subcell->list[p]; #else this=&(neighbour_i[j]); list_reset_f(this); @@ -228,10 +358,17 @@ int albe_potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS q=0; - while(neighbour_i[j][q]!=0) { + while(neighbour_i[k][q]!=-1) { - ktom=&(atom[neighbour_i[k][q]]); + ktom=&(itom[neighbour_i[k][q]]); q++; +#elif LOWMEM_LISTS + q=neighbour_i[k]; + + while(q!=-1) { + + ktom=&(itom[q]); + q=lc->subcell->list[q]; #else that=&(neighbour_i2[k]); list_reset_f(that); @@ -252,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) */ @@ -267,24 +412,27 @@ int albe_potential_force_calc(t_moldyn *moldyn) { Sk=params->S[brand_i]; Sk2=params->S2[brand_i]; /* albe needs i,k depending c,d,h and gamma values */ - exchange->gamma_i=&(params->gamma[brand_i]); - exchange->c_i=&(params->c[brand_i]); - exchange->d_i=&(params->d[brand_i]); - exchange->h_i=&(params->h[brand_i]); + 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 */ - exchange->gamma_i=&(params->gamma_m); - exchange->c_i=&(params->c_mixed); - exchange->d_i=&(params->d_mixed); - exchange->h_i=&(params->h_mixed); + 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; } - 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,&(ktom->r),&(ai->r)); @@ -307,12 +455,19 @@ int albe_potential_force_calc(t_moldyn *moldyn) { /* cos theta */ cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); - /* g_ijk */ + /* 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) { @@ -339,8 +494,12 @@ int albe_potential_force_calc(t_moldyn *moldyn) { /* increase k counter */ kcount++; +#endif // PTHREADS + #ifdef STATIC_LISTS } +#elif LOWMEM_LISTS + } #else } while(list_next_f(that)!=\ L_NO_NEXT_ELEMENT); @@ -448,10 +607,17 @@ if(moldyn->time>DSTART&&moldyn->timesubcell->list[q]; #else that=&(neighbour_i2[k]); list_reset_f(that); @@ -579,6 +745,8 @@ if(moldyn->time>DSTART&&moldyn->timetime>DSTART&&moldyn->time