X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=moldyn.c;h=7a9709e107a69c0f735eff1f01006611bcbb30ec;hb=e3e1d8e2f7fc65748c8da0accb9b3e6b0b1eb049;hp=eb3c4b03dfbcdb6b427c66d03aace481a6fe91ff;hpb=67df6b3a722a44e36fd56f3f040c3c9726b3fc3f;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index eb3c4b0..7a9709e 100644 --- a/moldyn.c +++ b/moldyn.c @@ -17,6 +17,12 @@ #include #include +#include + +#ifdef PARALLEL +#include +#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->t_tc; + scale=1.0+(scale-1.0)*moldyn->tau/moldyn->t_tc; scale=sqrt(scale); /* velocity scaling */ @@ -1182,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 { @@ -1292,6 +1308,8 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) { #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 @@ -1302,12 +1320,16 @@ 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 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); @@ -1331,6 +1353,17 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) { 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;icells;i++) list_init_f(&(lc->subcell[i])); @@ -1360,7 +1393,9 @@ int link_cell_update(t_moldyn *moldyn) { for(i=0;icells;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 @@ -1372,7 +1407,7 @@ int link_cell_update(t_moldyn *moldyn) { #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) { @@ -1381,6 +1416,9 @@ int link_cell_update(t_moldyn *moldyn) { } 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])); @@ -1774,7 +1812,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;igvir),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;iattr&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 +1938,7 @@ int potential_force_calc(t_moldyn *moldyn) { do { jtom=this->current->data; +#endif if(jtom==&(itom[i])) continue; @@ -1915,6 +1950,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 @@ -1940,7 +1978,7 @@ int potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS p=0; - while(neighbour_i[j][p]!=0) { + while(neighbour_i[j][p]!=-1) { jtom=&(atom[neighbour_i[j][p]]); p++; @@ -1984,7 +2022,7 @@ int potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS q=0; - while(neighbour_i[j][q]!=0) { + while(neighbour_i[k][q]!=-1) { ktom=&(atom[neighbour_i[k][q]]); q++; @@ -2039,7 +2077,7 @@ int potential_force_calc(t_moldyn *moldyn) { #ifdef STATIC_LISTS q=0; - while(neighbour_i[j][q]!=0) { + while(neighbour_i[k][q]!=-1) { ktom=&(atom[neighbour_i[k][q]]); q++; @@ -2114,6 +2152,9 @@ int potential_force_calc(t_moldyn *moldyn) { #endif /* some postprocessing */ +#ifdef PARALLEL + #pragma omp parallel for +#endif for(i=0;igvir.xx+=itom[i].r.x*itom[i].f.x; @@ -2125,7 +2166,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); } @@ -2328,11 +2369,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; @@ -2352,7 +2393,7 @@ int process_2b_bonds(t_moldyn *moldyn,void *data, #ifdef STATIC_LISTS p=0; - while(neighbour[j][p]!=0) { + while(neighbour[j][p]!=-1) { jtom=&(moldyn->atom[neighbour[j][p]]); p++; @@ -2800,3 +2841,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; +} +