X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=b73cd3432b65971aabc055e7f126fc5efdb3deb4;hb=37d949bdf2aaff618e70e6fd973ef350f4dacf12;hp=48129029a1c6e47a6fc0a59893c9fef04caa8999;hpb=c2697c8b6f808b7edef2f9e2ceded0aa65374168;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 4812902..b73cd34 100644
--- 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;
@@ -872,7 +888,7 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 		scale*=2.0;
 	else
 		if(moldyn->pt_scale&T_SCALE_BERENDSEN)
-			scale=1.0+(scale-1.0)*moldyn->timestep/moldyn->t_tc;
+			scale=1.0+(scale-1.0)*moldyn->tau/moldyn->t_tc;
 	scale=sqrt(scale);
 
 	/* velocity scaling */
@@ -1182,8 +1198,8 @@ int scale_volume(t_moldyn *moldyn) {
 
 	/* scaling factor */
 	if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
-		scale=(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->timestep;
-		scale=pow(1.0-scale,ONE_THIRD);
+		scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau;
+		scale=pow(scale,ONE_THIRD);
 	}
 	else {
 		scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
@@ -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;
+}
+