#include <time.h>
#include <math.h>
+#include <fpu_control.h>
+
+#ifdef PARALLEL
+#include <omp.h>
+#endif
+
#include "moldyn.h"
#include "report/report.h"
#include "potentials/tersoff.h"
#endif
+/* pse */
+#define PSE_NAME
+#define PSE_COL
+#include "pse.h"
+#undef PSE_NAME
+#undef PSE_COL
+
/*
* the moldyn functions
*/
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;
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 */
/* 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 {
#ifdef STATIC_LISTS
lc->subcell=malloc(lc->cells*sizeof(int*));
+#elif LOWMEM_LIST
+ lc->subcell=malloc(t_lowmem_list);
#else
lc->subcell=malloc(lc->cells*sizeof(t_list));
#endif
}
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
printf("[moldyn] initializing 'static' linked cells (%d)\n",
lc->cells);
+#elif LOWMEM_LIST
+ printf("[moldyn] initializing 'lowmem' linked cells (%d)\n",
+ lc->cells);
#else
printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
lc->cells);
i,lc->subcell[0],lc->subcell);
*/
}
+#elif LOWMEM_LIST
+ lc->subcell->head=malloc(lc->cells*sizeof(int));
+ if(lc->subcell->head==NULL) {
+ perror("[moldyn] head init (malloc)");
+ return -1;
+ }
+ lc->subcell->list=malloc(moldyn->count*sizeof(int));
+ if(lc->subcell->list==NULL) {
+ perror("[moldyn] list init (malloc)");
+ return -1;
+ }
#else
for(i=0;i<lc->cells;i++)
list_init_f(&(lc->subcell[i]));
for(i=0;i<lc->cells;i++)
#ifdef STATIC_LISTS
- memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
+ memset(lc->subcell[i],-1,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
+#elif LOWMEM_LIST
+ lc->subcell->head[i]=-1;
#else
list_destroy_f(&(lc->subcell[i]));
#endif
#ifdef STATIC_LISTS
p=0;
- while(lc->subcell[i+j*nx+k*nx*ny][p]!=0)
+ while(lc->subcell[i+j*nx+k*nx*ny][p]!=-1)
p++;
if(p>=MAX_ATOMS_PER_LIST) {
}
lc->subcell[i+j*nx+k*nx*ny][p]=count;
+#elif LOWMEM_LIST
+ lc->subcell->list[count]=list->subcell->head[i+j*nx+k*nx*ny];
+ list->subcell->head=count;
#else
list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
&(atom[count]));
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 */
memset(&(moldyn->gvir),0,sizeof(t_virial));
/* reset force, site energy and virial of every atom */
+#ifdef PARALLEL
+ i=omp_get_thread_num();
+ #pragma omp parallel for private(virial)
+#endif
for(i=0;i<count;i++) {
/* reset force */
#ifdef STATIC_LISTS
p=0;
- while(neighbour_i[j][p]!=0) {
+ while(neighbour_i[j][p]!=-1) {
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);
do {
jtom=this->current->data;
+#endif
if(jtom==&(itom[i]))
continue;
jtom,
bc_ij);
}
+#ifdef STATIC_LISTS
+ }
+#else
} while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
#endif
#ifdef STATIC_LISTS
p=0;
- while(neighbour_i[j][p]!=0) {
+ while(neighbour_i[j][p]!=-1) {
jtom=&(atom[neighbour_i[j][p]]);
p++;
#ifdef STATIC_LISTS
q=0;
- while(neighbour_i[j][q]!=0) {
+ while(neighbour_i[k][q]!=-1) {
ktom=&(atom[neighbour_i[k][q]]);
q++;
#ifdef STATIC_LISTS
q=0;
- while(neighbour_i[j][q]!=0) {
+ while(neighbour_i[k][q]!=-1) {
ktom=&(atom[neighbour_i[k][q]]);
q++;
#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;
/* 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);
}
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;
#ifdef STATIC_LISTS
p=0;
- while(neighbour[j][p]!=0) {
+ while(neighbour[j][p]!=-1) {
jtom=&(moldyn->atom[neighbour[j][p]]);
p++;
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;
+}
+