#include <omp.h>
#endif
+#ifdef PTHREAD
+#include <pthread.h>
+#endif
+
#include "../moldyn.h"
#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
+
int albe_potential_force_calc(t_moldyn *moldyn) {
int i,j,k,count;
#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];
#endif
u8 bc_ij,bc_ik;
int dnlc;
+#ifdef PTHREADS
+ int ret;
+ t_kdata kdata[27];
+ pthread_t kthread[27];
+#endif
// needed to work
t_atom *ai;
count=moldyn->count;
itom=moldyn->atom;
lc=&(moldyn->lc);
-#ifdef STATIC_LISTS
- atom=moldyn->atom;
-#endif
// optimized
params=moldyn->pot_params;
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
#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);
/* 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;
- 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);
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) */
#ifdef STATIC_LISTS
}
+#elif LOWMEM_LISTS
+ }
#else
} while(list_next_f(that)!=\
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 ... */
#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);
#ifdef STATIC_LISTS
}
+#elif LOWMEM_LISTS
+ }
#else
} while(list_next_f(that)!=\
L_NO_NEXT_ELEMENT);
#ifdef STATIC_LISTS
}
+#elif LOWMEM_LISTS
+ }
#else
} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
#endif