From: hackbard Date: Sun, 31 Aug 2008 16:02:10 +0000 (+0200) Subject: first attempts: parallel, albe_fast. fpu_ctrl. cleaned first j loop. X-Git-Url: https://www.hackdaworld.org/gitweb/?a=commitdiff_plain;h=cb2811d92cb79b1f24d11d4f4897426ce5f7dc3a;hp=21a29b60a222f81c311d0a3b85ccafa92c8b1279;p=physik%2Fposic.git first attempts: parallel, albe_fast. fpu_ctrl. cleaned first j loop. more ... --- diff --git a/moldyn.c b/moldyn.c index d56d6f3..b73cd34 100644 --- a/moldyn.c +++ b/moldyn.c @@ -17,6 +17,12 @@ #include #include +#include + +#ifdef PARALLEL +#include +#endif + #include "moldyn.h" #include "report/report.h" @@ -30,6 +36,13 @@ #include "potentials/tersoff.h" #endif +/* pse */ +#define PSE_NAME +#define PSE_COL +#include "pse.h" +#undef PSE_NAME +#undef PSE_COL + /* * the moldyn functions */ @@ -38,6 +51,9 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { printf("[moldyn] init\n"); + /* only needed if compiled without -msse2 (float-store prob!) */ + //fpu_set_rtd(); + memset(moldyn,0,sizeof(t_moldyn)); moldyn->argc=argc; @@ -1302,7 +1318,8 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) { } if(lc->cells<27) - printf("[moldyn] FATAL: less then 27 subcells!\n"); + printf("[moldyn] FATAL: less then 27 subcells! (%d)\n", + lc->cells); if(vol) { #ifdef STATIC_LISTS @@ -1774,7 +1791,11 @@ int velocity_verlet(t_moldyn *moldyn) { link_cell_update(moldyn); /* forces depending on chosen potential */ +#ifndef ALBE_FAST potential_force_calc(moldyn); +#else + albe_potential_force_calc(moldyn); +#endif for(i=0;igvir),0,sizeof(t_virial)); /* reset force, site energy and virial of every atom */ +#ifdef PARALLEL + #pragma omp parallel for private(virial) +#endif for(i=0;iattr&ATOM_ATTR_2BP)& - (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); - } - } #else this=&(neighbour_i[j]); list_reset_f(this); @@ -1904,6 +1916,7 @@ int potential_force_calc(t_moldyn *moldyn) { do { jtom=this->current->data; +#endif if(jtom==&(itom[i])) continue; @@ -1915,6 +1928,9 @@ int potential_force_calc(t_moldyn *moldyn) { jtom, bc_ij); } +#ifdef STATIC_LISTS + } +#else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); #endif @@ -2114,6 +2130,9 @@ int potential_force_calc(t_moldyn *moldyn) { #endif /* some postprocessing */ +#ifdef PARALLEL + #pragma omp parallel for +#endif for(i=0;igvir.xx+=itom[i].r.x*itom[i].f.x; @@ -2328,11 +2347,11 @@ int process_2b_bonds(t_moldyn *moldyn,void *data, int p; #else t_list neighbour[27]; + t_list *this; #endif u8 bc; t_atom *itom,*jtom; int i,j; - t_list *this; lc=&(moldyn->lc); itom=moldyn->atom; @@ -2800,3 +2819,22 @@ int visual_atoms(t_moldyn *moldyn) { return 0; } +/* + * fpu cntrol functions + */ + +// set rounding to double (eliminates -ffloat-store!) +int fpu_set_rtd(void) { + + fpu_control_t ctrl; + + _FPU_GETCW(ctrl); + + ctrl&=~_FPU_EXTENDED; + ctrl|=_FPU_DOUBLE; + + _FPU_SETCW(ctrl); + + return 0; +} + diff --git a/potentials/albe_fast.c b/potentials/albe_fast.c new file mode 100644 index 0000000..3b5f681 --- /dev/null +++ b/potentials/albe_fast.c @@ -0,0 +1,638 @@ +/* + * test: albe_new.c + * + * author: Frank Zirkelbach + * + */ + +#define _GNU_SOURCE +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef PARALLEL +#include +#endif + +#include "../moldyn.h" +#include "../math/math.h" +#include "albe.h" + +/* + * virial calculation + */ + +#define albe_v_calc(a,f,d) a->virial.xx+=f->x*d->x; \ + a->virial.yy+=f->y*d->y; \ + a->virial.zz+=f->z*d->z; \ + a->virial.xy+=f->x*d->y; \ + a->virial.xz+=f->x*d->z; \ + a->virial.yz+=f->y*d->z + +int albe_potential_force_calc(t_moldyn *moldyn) { + + int i,j,k,count; + t_atom *itom,*jtom,*ktom; + t_virial *virial; + t_linkcell *lc; +#ifdef STATIC_LISTS + int *neighbour_i[27]; + int p,q; + t_atom *atom; +#else + t_list neighbour_i[27]; + t_list neighbour_i2[27]; + t_list *this,*that; +#endif + u8 bc_ij,bc_ik; + int dnlc; + + // needed to work + t_atom *ai; + + // optimized + t_albe_mult_params *params; + t_albe_exchange *exchange; + t_3dvec dist_ij; + double d_ij2; + double d_ij; + u8 brand_i; + double S2; + int kcount; + double zeta_ij; + double pre_dzeta; + + // more ... + double Rk,Sk,Sk2; + t_3dvec dist_ik; + double d_ik2,d_ik; + double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg; + double f_c_ik,df_c_ik; + + t_3dvec force; + double f_a,df_a,b,db,f_c,df_c; + double f_r,df_r; + double scale; + double mu,B; + double lambda,A; + double r0; + double S,R; + double energy; + + double dijdik_inv,fcdg,dfcg; + t_3dvec dcosdrj,dcosdrk; + t_3dvec tmp; + + + count=moldyn->count; + itom=moldyn->atom; + lc=&(moldyn->lc); +#ifdef STATIC_LISTS + atom=moldyn->atom; +#endif + + // optimized + params=moldyn->pot_params; + exchange=&(params->exchange); + + + /* reset energy */ + moldyn->energy=0.0; + + /* reset global virial */ + memset(&(moldyn->gvir),0,sizeof(t_virial)); + + /* reset force, site energy and virial of every atom */ +#ifdef PARALLEL + #pragma omp parallel for private(virial) +#endif + for(i=0;ixx=0.0; + virial->yy=0.0; + virial->zz=0.0; + virial->xy=0.0; + virial->xz=0.0; + virial->yz=0.0; + + /* reset site energy */ + itom[i].e=0.0; + + } + + /* get energy, force and virial of every atom */ + + /* first (and only) loop over atoms i */ + for(i=0;idim.x/2)/lc->x, + (itom[i].r.y+moldyn->dim.y/2)/lc->y, + (itom[i].r.z+moldyn->dim.z/2)/lc->z, + neighbour_i); + + dnlc=lc->dnlc; + + /* copy the neighbour lists */ +#ifndef STATIC_LISTS + memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); +#endif + + ai=&(itom[i]); + brand_i=ai->brand; + + /* loop over atoms j */ + for(j=0;j<27;j++) { + + bc_ij=(jstart==NULL) + continue; + + do { + + jtom=this->current->data; +#endif + + if(jtom==&(itom[i])) + continue; + + if(!(jtom->attr&ATOM_ATTR_3BP)) + continue; + + /* reset 3bp run */ + moldyn->run3bp=1; + + +/* j1 func here ... */ +/* albe 3 body potential function (first ij loop) */ + + /* reset zeta sum */ + zeta_ij=0.0; + + /* + * set ij depending values + */ + + if(brand_i==jtom->brand) { + S2=params->S2[brand_i]; + } + else { + S2=params->S2mixed; + } + + /* dist_ij, d_ij2 */ + v3_sub(&dist_ij,&(jtom->r),&(ai->r)); + if(bc_ij) check_per_bound(moldyn,&dist_ij); + d_ij2=v3_absolute_square(&dist_ij); + + /* if d_ij2 > S2 => no force & potential energy contribution */ + if(d_ij2>S2) + continue; + + /* d_ij */ + d_ij=sqrt(d_ij2); + + /* reset k counter for first k loop */ + kcount=0; + + /* first loop over atoms k */ + for(k=0;k<27;k++) { + + bc_ik=(kstart==NULL) + continue; + + do { + ktom=that->current->data; +#endif + + if(!(ktom->attr&ATOM_ATTR_3BP)) + continue; + + if(ktom==jtom) + continue; + + if(ktom==&(itom[i])) + continue; + + +/* k1 func here ... */ +/* albe 3 body potential function (first k loop) */ + + 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 */ + 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]); + } + 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); + } + 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)); + 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.. + + /* 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++; + +#ifdef STATIC_LISTS + } +#else + } while(list_next_f(that)!=\ + L_NO_NEXT_ELEMENT); +#endif + + } + + +/* j2 func here ... */ + + + if(brand_i==jtom->brand) { + S=params->S[brand_i]; + R=params->R[brand_i]; + B=params->B[brand_i]; + A=params->A[brand_i]; + r0=params->r0[brand_i]; + mu=params->mu[brand_i]; + lambda=params->lambda[brand_i]; + } + else { + S=params->Smixed; + R=params->Rmixed; + B=params->Bmixed; + A=params->Amixed; + r0=params->r0_mixed; + mu=params->mu_m; + lambda=params->lambda_m; + } + + /* f_c, df_c */ + if(d_ijf),&(ai->f),&force); + + /* force contribution for atom j */ + v3_scale(&force,&force,-1.0); // dri rij = - drj rij + v3_add(&(jtom->f),&(jtom->f),&force); + + /* virial */ + virial_calc(ai,&force,&(dist_ij)); + +#ifdef DEBUG +if(moldyn->time>DSTART&&moldyn->timeatom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) { + printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag); + printf(" adding %f %f %f\n",force.x,force.y,force.z); + if(ai==&(moldyn->atom[0])) + printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); + if(jtom==&(moldyn->atom[0])) + printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z); + printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r), + f_c,b,f_a,f_r); + printf(" %f %f %f\n",zeta_ij,.0,.0); + } +} +#endif + + /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */ + pre_dzeta=0.5*f_a*f_c*db; + + /* energy contribution */ + energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism + moldyn->energy+=energy; + ai->e+=energy; + + /* reset k counter for second k loop */ + kcount=0; + + + /* second loop over atoms k */ + for(k=0;k<27;k++) { + + bc_ik=(kstart==NULL) + continue; + + do { + ktom=that->current->data; +#endif + + if(!(ktom->attr&ATOM_ATTR_3BP)) + continue; + + if(ktom==jtom) + continue; + + if(ktom==&(itom[i])) + continue; + + +/* k2 func here ... */ +/* albe 3 body potential function (second k loop) */ + + if(kcount>ALBE_MAXN) + printf("FATAL: neighbours!\n"); + + /* d_ik2 */ + d_ik2=exchange->d_ik2[kcount]; + + if(brand_i==ktom->brand) + Sk2=params->S2[brand_i]; + else + Sk2=params->S2mixed; + + /* return if d_ik > S */ + if(d_ik2>Sk2) { + kcount++; + continue; + } + + /* dist_ik, d_ik */ + dist_ik=exchange->dist_ik[kcount]; + d_ik=exchange->d_ik[kcount]; + + /* f_c_ik, df_c_ik */ + f_c_ik=exchange->f_c_ik[kcount]; + df_c_ik=exchange->df_c_ik[kcount]; + + /* g, dg, cos_theta */ + g=exchange->g[kcount]; + dg=exchange->dg[kcount]; + cos_theta=exchange->cos_theta[kcount]; + + /* cos_theta derivatives wrt j,k */ + dijdik_inv=1.0/(d_ij*d_ik); + v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j + v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2); + v3_add(&dcosdrj,&dcosdrj,&tmp); + v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k + v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2); + v3_add(&dcosdrk,&dcosdrk,&tmp); + + /* f_c_ik * dg, df_c_ik * g */ + fcdg=f_c_ik*dg; + dfcg=df_c_ik*g; + + /* derivative wrt j */ + v3_scale(&force,&dcosdrj,fcdg*pre_dzeta); + + /* force contribution */ + v3_add(&(jtom->f),&(jtom->f),&force); + +#ifdef DEBUG +if(moldyn->time>DSTART&&moldyn->timeatom[DATOM])) { + printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag); + printf(" adding %f %f %f\n",force.x,force.y,force.z); + printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z); + printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI)); + printf(" d ij ik = %f %f\n",d_ij,d_ik); + } +} +#endif + + /* virial */ + virial_calc(ai,&force,&dist_ij); + + /* force contribution to atom i */ + v3_scale(&force,&force,-1.0); + v3_add(&(ai->f),&(ai->f),&force); + + /* derivative wrt k */ + v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik + v3_scale(&tmp,&dcosdrk,fcdg); + v3_add(&force,&force,&tmp); + v3_scale(&force,&force,pre_dzeta); + + /* force contribution */ + v3_add(&(ktom->f),&(ktom->f),&force); + +#ifdef DEBUG +if(moldyn->time>DSTART&&moldyn->timeatom[DATOM])) { + printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag); + printf(" adding %f %f %f\n",force.x,force.y,force.z); + printf(" total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z); + printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI)); + printf(" d ij ik = %f %f\n",d_ij,d_ik); + } +} +#endif + + /* virial */ + virial_calc(ai,&force,&dist_ik); + + /* force contribution to atom i */ + v3_scale(&force,&force,-1.0); + v3_add(&(ai->f),&(ai->f),&force); + + /* increase k counter */ + kcount++; + + + +#ifdef STATIC_LISTS + } +#else + } while(list_next_f(that)!=\ + L_NO_NEXT_ELEMENT); +#endif + + } + +#ifdef STATIC_LISTS + } +#else + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif + + } + +#ifdef DEBUG + //printf("\n\n"); +#endif +#ifdef VDEBUG + printf("\n\n"); +#endif + + } + +#ifdef DEBUG + //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z); + if(moldyn->time>DSTART&&moldyn->timeatom[DATOM].f.x); + printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y); + printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z); + } +#endif + + /* some postprocessing */ +#ifdef PARALLEL + #pragma omp parallel for +#endif + for(i=0;igvir.xx+=itom[i].r.x*itom[i].f.x; + moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y; + moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z; + moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x; + moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x; + moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y; + + /* check forces regarding the given timestep */ + if(v3_norm(&(itom[i].f))>\ + 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square) + printf("[moldyn] WARNING: pfc (high force: atom %d)\n", + i); + } + + return 0; +} +