virial calc define + default amount threads -> 2

[physik/posic.git] / potentials / albe_fast.c

diff --git a/potentials/albe_fast.c b/potentials/albe_fast.c
index 952f1fe..0db5f66 100644 (file)
--- a/potentials/albe_fast.c
+++ b/potentials/albe_fast.c
@@ -23,7 +23,7 @@
 
 #ifdef PTHREADS
 #include <pthread.h>
-#define MAX_THREADS 4
+#define MAX_THREADS 2
 #endif
 
 #include "../moldyn.h"
@@ -39,12 +39,12 @@ extern pthread_mutex_t emutex;
  * 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
 
@@ -461,7 +461,8 @@ int albe_potential_force_calc(t_moldyn *moldyn) {
        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) {
@@ -595,7 +596,8 @@ 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);
@@ -623,7 +625,8 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 #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);
@@ -703,7 +706,7 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 
 typedef struct s_pft_data {
        t_moldyn *moldyn;
-       int i;
+       int start,end;
 } t_pft_data;
 
 void *potential_force_thread(void *ptr) {
@@ -785,14 +788,12 @@ void *potential_force_thread(void *ptr) {
        // optimized
        params=moldyn->pot_params;
 
+       /* get energy, force and virial for atoms */
 
-       /* get energy, force and virial of every atom */
-
-       /* first (and only) loop over atoms i */
-       for(i=0;i<count;i++) {
+       for(i=pft_data->start;i<pft_data->end;i++) {
 
                if(!(itom[i].attr&ATOM_ATTR_3BP))
-                       continue;
+                       return 0;
 
                link_cell_neighbour_index(moldyn,
                                          (itom[i].r.x+moldyn->dim.x/2)/lc->x,
@@ -1096,7 +1097,8 @@ void *potential_force_thread(void *ptr) {
 
        /* virial */
        pthread_mutex_lock(&(amutex[ai->tag])); 
-       virial_calc(ai,&force,&(dist_ij));
+       albe_v_calc(ai,&force,&(dist_ij));
+       //virial_calc(ai,&force,&(dist_ij));
        pthread_mutex_unlock(&(amutex[ai->tag]));       
 
 #ifdef DEBUG
@@ -1238,7 +1240,8 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 
        /* virial */
        pthread_mutex_lock(&(amutex[ai->tag])); 
-       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);
@@ -1270,7 +1273,8 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 
        /* virial */
        pthread_mutex_lock(&(amutex[ai->tag])); 
-       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);
@@ -1302,6 +1306,8 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
 #endif
                
                }
+
+       } // i loop
                
 #ifdef DEBUG
        //printf("\n\n");
@@ -1310,8 +1316,6 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
        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) {
@@ -1365,49 +1369,37 @@ int albe_potential_force_calc(t_moldyn *moldyn) {
 
        }
 
-       i=0;
-       while(i<count) {
-
-               /* start threads */
-               for(j=0;j<MAX_THREADS;j++) {
-
-                       /* break if all atoms are processed */
-                       if(j+i==count)
-                               break;
-
-                       /* prepare thread data */
-                       pft_data[j].moldyn=moldyn;
-                       pft_data[j].i=j+i;
+       /* start threads */
+       for(j=0;j<MAX_THREADS;j++) {
 
-                       ret=pthread_create(&(pft_thread[j]),NULL,
-                                           potential_force_thread,
-                                           &(pft_data[j]));
-                       if(ret)  {
-                               perror("[albe fast] pf thread create");
-                               return ret;
-                       }
+               /* prepare thread data */
+               pft_data[j].moldyn=moldyn;
+               pft_data[j].start=j*(count/MAX_THREADS);
+               if(j==MAX_THREADS-1) {
+                       pft_data[j].end=count;
                }
-
-               //printf("threads created! %d\n",j);
-
-               /* join threads */
-               for(j=0;j<MAX_THREADS;j++) {
-
-                       if(j+i==count)
-                               break;
-
-                       ret=pthread_join(pft_thread[j],NULL);
-                       if(ret) {
-                               perror("[albe fast] pf thread join");
-                               return ret;
-                       }
+               else {
+                       pft_data[j].end=pft_data[j].start;
+                       pft_data[j].end+=count/MAX_THREADS;
                }
 
-               /* increment counter */
-               i+=MAX_THREADS;
+               ret=pthread_create(&(pft_thread[j]),NULL,
+                                   potential_force_thread,
+                                   &(pft_data[j]));
+               if(ret)  {
+                       perror("[albe fast] pf thread create");
+                       return ret;
+               }
+       }
 
-               //printf("threads joined! -> %d\n",i);
+       /* join threads */
+       for(j=0;j<MAX_THREADS;j++) {
 
+               ret=pthread_join(pft_thread[j],NULL);
+               if(ret) {
+                       perror("[albe fast] pf thread join");
+                       return ret;
+               }
        }
 
        /* some postprocessing */
@@ -1427,8 +1419,6 @@ int albe_potential_force_calc(t_moldyn *moldyn) {
                               i);
        }
 
-       pthread_exit(NULL);
-
        return 0;
 }