first attempts: parallel, albe_fast. fpu_ctrl. cleaned first j loop.

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index d56d6f3..b73cd34 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -17,6 +17,12 @@
 #include <time.h>
 #include <math.h>
 
+#include <fpu_control.h>
+
+#ifdef PARALLEL
+#include <omp.h>
+#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;i<count;i++) {
                /* check whether fixed atom */
@@ -1829,6 +1850,9 @@ int potential_force_calc(t_moldyn *moldyn) {
        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;i<count;i++) {
 
                /* reset force */
@@ -1883,18 +1907,6 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                        jtom=&(atom[neighbour_i[j][p]]);
                                        p++;
-
-                                       if(jtom==&(itom[i]))
-                                               continue;
-
-                                       if((jtom->attr&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;i<count;i++) {
                /* calculate global virial */
                moldyn->gvir.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;
+}
+