#include "../math/math.h"
#include "albe.h"
-#ifdef PTHREADS
-typedef struct s_kdata {
- t_moldyn *moldyn;
- t_atom *itom,*jtom;
-} t_kdata;
-#endif
-
/*
* virial calculation
*/
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
+#ifndef PTHREADS
int albe_potential_force_calc(t_moldyn *moldyn) {
/* first loop over atoms k */
for(k=0;k<27;k++) {
-#ifdef PTHREADS
- // create threads
- kdata.moldyn=moldyn;
- kdata.jtom=jtom;
- kdata.itom=&(itom[i]);
- ret=pthread_create(&(kthread[k]),NULL,k1_thread,&(kdata[k]));
- if(ret) {
- perror("[albe fast] thread create");
- return ret;
- }
-#else
bc_ik=(k<dnlc)?0:1;
#ifdef STATIC_LISTS
q=0;
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) */
L_NO_NEXT_ELEMENT);
#endif
-#endif // PTHREADS
-
}
-#ifdef PTHREADS
- // join threads
- for(k=0;k<27;k++) {
- ret=pthread_join(kthread[k],NULL);
- if(ret) {
- perror("[albe fast] join thread");
- return ret;
- }
- }
-#endif
-
-
/* j2 func here ... */
return 0;
}
+
+#else // PTHREADS
+
+
+typedef struct s_pft_data {
+ t_moldyn *moldyn;
+ int i;
+} t_pft_data;
+
+void *potential_force_thread(void *ptr) {
+
+ t_pft_data *pft_data;
+ t_moldyn *moldyn;
+ t_albe_exchange ec;
+
+ int i,j,k,count;
+ t_atom *itom,*jtom,*ktom;
+ t_linkcell *lc;
+#ifdef STATIC_LISTS
+ int *neighbour_i[27];
+ int p,q;
+#elif LOWMEM_LISTS
+ int neighbour_i[27];
+ int p,q;
+#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;
+
+ // even more ...
+ double gamma_i;
+ double c_i;
+ double d_i;
+ double h_i;
+ double ci2;
+ double di2;
+ double ci2di2;
+
+ pft_data=ptr;
+ moldyn=pft_data->moldyn;
+ exchange=&ec;
+
+ count=moldyn->count;
+ itom=moldyn->atom;
+ lc=&(moldyn->lc);
+
+ // optimized
+ params=moldyn->pot_params;
+
+
+ /* get energy, force and virial of every atom */
+
+ /* first (and only) loop over atoms i */
+ for(i=0;i<count;i++) {
+
+ if(!(itom[i].attr&ATOM_ATTR_3BP))
+ continue;
+
+ link_cell_neighbour_index(moldyn,
+ (itom[i].r.x+moldyn->dim.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 */
+#ifdef STATIC_LISTS
+#elif LOWMEM_LISTS
+#else
+ 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=(j<dnlc)?0:1;
+#ifdef STATIC_LISTS
+ p=0;
+
+ while(neighbour_i[j][p]!=-1) {
+
+ 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);
+
+ if(this->start==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=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+ q=0;
+
+ while(neighbour_i[k][q]!=-1) {
+
+ 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);
+
+ if(that->start==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 */
+ 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++;
+
+#ifdef STATIC_LISTS
+ }
+#elif LOWMEM_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_ij<R) {
+ f_c=1.0;
+ df_c=0.0;
+ }
+ else {
+ s_r=S-R;
+ arg=M_PI*(d_ij-R)/s_r;
+ f_c=0.5+0.5*cos(arg);
+ df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
+ }
+
+ /* f_a, df_a */
+ f_a=-B*exp(-mu*(d_ij-r0));
+ df_a=mu*f_a/d_ij;
+
+ /* f_r, df_r */
+ f_r=A*exp(-lambda*(d_ij-r0));
+ df_r=lambda*f_r/d_ij;
+
+ /* b, db */
+ if(zeta_ij==0.0) {
+ b=1.0;
+ db=0.0;
+ }
+ else {
+ b=1.0/sqrt(1.0+zeta_ij);
+ db=-0.5*b/(1.0+zeta_ij);
+ }
+
+ /* force contribution for atom i */
+ scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
+ v3_scale(&force,&(dist_ij),scale);
+ v3_add(&(ai->f),&(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->time<DEND) {
+ if((ai==&(moldyn->atom[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=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+ q=0;
+
+ while(neighbour_i[k][q]!=-1) {
+
+ 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);
+
+ if(that->start==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->time<DEND) {
+ if(jtom==&(moldyn->atom[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->time<DEND) {
+ if(ktom==&(moldyn->atom[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
+ }
+#elif LOWMEM_LISTS
+ }
+#else
+ } while(list_next_f(that)!=\
+ L_NO_NEXT_ELEMENT);
+#endif
+
+ }
+
+#ifdef STATIC_LISTS
+ }
+#elif LOWMEM_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->time<DEND) {
+ printf("force:\n");
+ printf(" x: %0.40f\n",moldyn->atom[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;i<count;i++) {
+ /* calculate global virial */
+ moldyn->gvir.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;
+}
+
+int albe_potential_force_calc(t_moldyn *moldyn) {
+
+ int i,ret;
+ t_pft_data *pft_data;
+ int count;
+ pthread_t *pft_thread;
+ t_atom *itom;
+ t_virial *virial;
+
+ count=moldyn->count;
+ itom=moldyn->atom;
+
+ /* 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 */
+ for(i=0;i<count;i++) {
+
+ /* reset force */
+ v3_zero(&(itom[i].f));
+
+ /* reset virial */
+ virial=&(itom[i].virial);
+ virial->xx=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;
+
+ }
+
+ /* alloc thread memory */
+ pft_thread=malloc(count*sizeof(pthread_t));
+ if(pft_thread==NULL) {
+ perror("[albe fast] alloc thread mem");
+ return -1;
+ }
+ pft_data=malloc(count*sizeof(t_pft_data));
+ if(pft_data==NULL) {
+ perror("[albe fast] alloc thread mem");
+ return -1;
+ }
+
+ /* start threads */
+ for(i=0;i<count;i++) {
+
+ /* prepare thread data */
+ pft_data[i].moldyn=moldyn;
+ pft_data[i].i=i;
+
+ ret=pthread_create(&(pft_thread[i]),NULL,
+ potential_force_thread,&(pft_data[i]));
+ if(ret) {
+ perror("[albe fast] pf thread create");
+ return ret;
+ }
+ }
+
+ /* join threads */
+ for(i=0;i<count;i++) {
+ ret=pthread_join(pft_thread[i],NULL);
+ if(ret) {
+ perror("[albe fast] pf thread join");
+ return ret;
+ }
+ }
+
+ return 0;
+}
+
+#endif // PTHREADS
projects / physik / posic.git / commitdiff