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

diff --git a/moldyn.c b/moldyn.c
index dd8a527..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;
@@ -485,6 +501,10 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 	count=moldyn->count;
 
 	/* how many atoms do we expect */
+	if(type==NONE) {
+		new*=1;
+		printf("[moldyn] WARNING: create 'none' lattice called");
+	}
 	if(type==CUBIC) new*=1;
 	if(type==FCC) new*=4;
 	if(type==DIAMOND) new*=8;
@@ -868,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->t_tc;
+			scale=1.0+(scale-1.0)*moldyn->tau/moldyn->t_tc;
 	scale=sqrt(scale);
 
 	/* velocity scaling */
@@ -1178,7 +1198,7 @@ int scale_volume(t_moldyn *moldyn) {
 
 	/* scaling factor */
 	if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
-		scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc;
+		scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau;
 		scale=pow(scale,ONE_THIRD);
 	}
 	else {
@@ -1298,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
@@ -1574,9 +1595,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n");
 	if(moldyn->cutoff>0.5*moldyn->dim.z)
 		printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n");
-	ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
-	if(ds>0.05*moldyn->nnd)
+	if(moldyn->count) {
+		ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
+		if(ds>0.05*moldyn->nnd)
 		printf("[moldyn] WARNING: forces too high / tau too small!\n");
+	}
 
 	/* zero absolute time */
 	// should have right values!
@@ -1698,6 +1721,7 @@ printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
        sched->count,i,moldyn->total_steps,
        moldyn->t,moldyn->t_avg,
        moldyn->p/BAR,moldyn->p_avg/BAR,
+       //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
        moldyn->volume,
        (int)(t2.tv_sec-t1.tv_sec));
 
@@ -1767,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 */
@@ -1822,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 */
@@ -1876,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);
@@ -1897,6 +1916,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 				do {
 					jtom=this->current->data;
+#endif
 
 					if(jtom==&(itom[i]))
 						continue;
@@ -1908,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
 
@@ -2107,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;
@@ -2118,7 +2144,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 		/* check forces regarding the given timestep */
 		if(v3_norm(&(itom[i].f))>\
-		   0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
+                   0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
 			printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
 			       i);
 	}
@@ -2321,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;
@@ -2793,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;
+}
+