#ifdef PTHREADS
#include <pthread.h>
-#define MAX_THREADS 4
+#define MAX_THREADS 2
#endif
#include "../moldyn.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
+#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
#ifndef PTHREADS
}
/* force contribution for atom i */
+#ifdef MATTONI
+ scale=-0.5*(f_c*(df_r-b*df_a)); // - in albe formalism
+#else
scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
+#endif
v3_scale(&force,&(dist_ij),scale);
v3_add(&(ai->f),&(ai->f),&force);
v3_add(&(jtom->f),&(jtom->f),&force);
/* virial */
- virial_calc(ai,&force,&(dist_ij));
+ albe_v_calc(ai,&force,&(dist_ij));
+ //virial_calc(ai,&force,&(dist_ij));
#ifdef DEBUG
if(moldyn->time>DSTART&&moldyn->time<DEND) {
#endif
/* virial */
- virial_calc(ai,&force,&dist_ij);
+ albe_v_calc(ai,&force,&dist_ij);
+ //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 */
+#ifdef MATTONI
+ v3_scale(&tmp,&dcosdrk,fcdg);
+ v3_scale(&force,&tmp,pre_dzeta);
+#else
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);
+#endif
/* force contribution */
v3_add(&(ktom->f),&(ktom->f),&force);
#endif
/* virial */
- virial_calc(ai,&force,&dist_ik);
+ albe_v_calc(ai,&force,&dist_ik);
+ //virial_calc(ai,&force,&dist_ik);
/* force contribution to atom i */
v3_scale(&force,&force,-1.0);
if(!(itom[i].attr&ATOM_ATTR_3BP))
return 0;
- link_cell_neighbour_index(moldyn,
+ // thread safe this way!
+ dnlc=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
while(p!=-1) {
jtom=&(itom[p]);
- p=lc->subcell->list[p];
+ p=lc->subcell->list[p]; // thread safe!
#else
this=&(neighbour_i[j]);
list_reset_f(this);
}
/* force contribution for atom i */
+#ifdef MATTONI
+ scale=-0.5*(f_c*(df_r-b*df_a)); // - in albe formalism
+#else
scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
+#endif
v3_scale(&force,&(dist_ij),scale);
- pthread_mutex_lock(&(amutex[ai->tag]));
+ if(pthread_mutex_lock(&(amutex[ai->tag])))
+ perror("[albe fast] mutex lock (1)\n");
v3_add(&(ai->f),&(ai->f),&force);
- pthread_mutex_unlock(&(amutex[ai->tag]));
+ if(pthread_mutex_unlock(&(amutex[ai->tag])))
+ perror("[albe fast] mutex unlock (1)\n");
/* force contribution for atom j */
v3_scale(&force,&force,-1.0); // dri rij = - drj rij
- pthread_mutex_lock(&(amutex[jtom->tag]));
+ if(pthread_mutex_lock(&(amutex[jtom->tag])))
+ perror("[albe fast] mutex lock (2)\n");
v3_add(&(jtom->f),&(jtom->f),&force);
- pthread_mutex_unlock(&(amutex[jtom->tag]));
+ if(pthread_mutex_unlock(&(amutex[jtom->tag])))
+ perror("[albe fast] mutex unlock (2)\n");
/* virial */
- pthread_mutex_lock(&(amutex[ai->tag]));
- virial_calc(ai,&force,&(dist_ij));
- pthread_mutex_unlock(&(amutex[ai->tag]));
+ if(pthread_mutex_lock(&(amutex[ai->tag])))
+ perror("[albe fast] mutex lock (3)\n");
+ albe_v_calc(ai,&force,&(dist_ij));
+ //virial_calc(ai,&force,&(dist_ij));
+ if(pthread_mutex_unlock(&(amutex[ai->tag])))
+ perror("[albe fast] mutex unlock (3)\n");
#ifdef DEBUG
if(moldyn->time>DSTART&&moldyn->time<DEND) {
/* energy contribution */
energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
- pthread_mutex_lock(&emutex);
+ if(pthread_mutex_lock(&emutex))
+ perror("[albe fast] mutex lock (energy)\n");
moldyn->energy+=energy;
- pthread_mutex_unlock(&emutex);
- pthread_mutex_lock(&(amutex[ai->tag]));
+ if(pthread_mutex_unlock(&emutex))
+ perror("[albe fast] mutex unlock (energy)\n");
+ if(pthread_mutex_lock(&(amutex[ai->tag])))
+ perror("[albe fast] mutex lock (4)\n");
ai->e+=energy;
- pthread_mutex_unlock(&(amutex[ai->tag]));
+ if(pthread_mutex_unlock(&(amutex[ai->tag])))
+ perror("[albe fast] mutex unlock (4)\n");
/* reset k counter for second k loop */
kcount=0;
v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
/* force contribution */
- pthread_mutex_lock(&(amutex[jtom->tag]));
+ if(pthread_mutex_lock(&(amutex[jtom->tag])))
+ perror("[albe fast] mutex lock (5)\n");
v3_add(&(jtom->f),&(jtom->f),&force);
- pthread_mutex_unlock(&(amutex[jtom->tag]));
+ if(pthread_mutex_unlock(&(amutex[jtom->tag])))
+ perror("[albe fast] mutex unlock (5)\n");
#ifdef DEBUG
if(moldyn->time>DSTART&&moldyn->time<DEND) {
#endif
/* virial */
- pthread_mutex_lock(&(amutex[ai->tag]));
- virial_calc(ai,&force,&dist_ij);
+ if(pthread_mutex_lock(&(amutex[ai->tag])))
+ perror("[albe fast] mutex lock (6)\n");
+ albe_v_calc(ai,&force,&dist_ij);
+ //virial_calc(ai,&force,&dist_ij);
/* force contribution to atom i */
v3_scale(&force,&force,-1.0);
v3_add(&(ai->f),&(ai->f),&force);
- pthread_mutex_unlock(&(amutex[ai->tag]));
+ if(pthread_mutex_unlock(&(amutex[ai->tag])))
+ perror("[albe fast] mutex unlock (6)\n");
/* derivative wrt k */
+#ifdef MATTONI
+ v3_scale(&tmp,&dcosdrk,fcdg);
+ v3_scale(&force,&tmp,pre_dzeta);
+#else
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);
+#endif
/* force contribution */
- pthread_mutex_lock(&(amutex[ktom->tag]));
+ if(pthread_mutex_lock(&(amutex[ktom->tag])))
+ perror("[albe fast] mutex lock (7)\n");
v3_add(&(ktom->f),&(ktom->f),&force);
- pthread_mutex_unlock(&(amutex[ktom->tag]));
+ if(pthread_mutex_unlock(&(amutex[ktom->tag])))
+ perror("[albe fast] mutex unlock (7)\n");
#ifdef DEBUG
if(moldyn->time>DSTART&&moldyn->time<DEND) {
#endif
/* virial */
- pthread_mutex_lock(&(amutex[ai->tag]));
- virial_calc(ai,&force,&dist_ik);
+ if(pthread_mutex_lock(&(amutex[ai->tag])))
+ perror("[albe fast] mutex lock (8)\n");
+ albe_v_calc(ai,&force,&dist_ik);
+ //virial_calc(ai,&force,&dist_ik);
/* force contribution to atom i */
v3_scale(&force,&force,-1.0);
v3_add(&(ai->f),&(ai->f),&force);
- pthread_mutex_unlock(&(amutex[ai->tag]));
+ if(pthread_mutex_unlock(&(amutex[ai->tag])))
+ perror("[albe fast] mutex unlock (8)\n");
/* increase k counter */
kcount++;
/* prepare thread data */
pft_data[j].moldyn=moldyn;
- pft_data[j].start=j*count/MAX_THREADS;
+ pft_data[j].start=j*(count/MAX_THREADS);
if(j==MAX_THREADS-1) {
pft_data[j].end=count;
}