X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=moldyn.c;h=7a9709e107a69c0f735eff1f01006611bcbb30ec;hb=e3e1d8e2f7fc65748c8da0accb9b3e6b0b1eb049;hp=6fa9df25c008ac26c5e12e30a1ca66095b3a600f;hpb=2dddad627c54f294c4c5b5d736178c3ef3288159;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 6fa9df2..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,100 +36,12 @@ #include "potentials/tersoff.h" #endif - -/* - * global variables, pse and atom colors (only needed here) - */ - -static char *pse_name[]={ - "*", - "H", - "He", - "Li", - "Be", - "B", - "C", - "N", - "O", - "F", - "Ne", - "Na", - "Mg", - "Al", - "Si", - "P", - "S", - "Cl", - "Ar", -}; - -static char *pse_col[]={ - "*", - "White", - "He", - "Li", - "Be", - "B", - "Gray", - "N", - "Blue", - "F", - "Ne", - "Na", - "Mg", - "Al", - "Yellow", - "P", - "S", - "Cl", - "Ar", -}; - -/* -static double pse_mass[]={ - 0, - 0, - 0, - 0, - 0, - 0, - M_C, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - M_SI, - 0, - 0, - 0, - 0, -}; - -static double pse_lc[]={ - 0, - 0, - 0, - 0, - 0, - 0, - LC_C, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - LC_SI, - 0, - 0, - 0, - 0, -}; -*/ +/* pse */ +#define PSE_NAME +#define PSE_COL +#include "pse.h" +#undef PSE_NAME +#undef PSE_COL /* * the moldyn functions @@ -133,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; @@ -208,7 +129,7 @@ int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) { moldyn->pt_scale|=ptype; moldyn->p_tc=ptc; - printf("[moldyn] p/t scaling:\n"); + printf("[moldyn] p scaling:\n"); printf(" p: %s",ptype?"yes":"no "); if(ptype) @@ -224,7 +145,7 @@ int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) { moldyn->pt_scale|=ttype; moldyn->t_tc=ttc; - printf("[moldyn] p/t scaling:\n"); + printf("[moldyn] t scaling:\n"); printf(" t: %s",ttype?"yes":"no "); if(ttype) @@ -315,7 +236,7 @@ int set_potential(t_moldyn *moldyn,u8 type) { moldyn->func3b_k1=tersoff_mult_3bp_k1; moldyn->func3b_j2=tersoff_mult_3bp_j2; moldyn->func3b_k2=tersoff_mult_3bp_k2; - // missing: check 2b bond func + moldyn->check_2b_bond=tersoff_mult_check_2b_bond; break; case MOLDYN_POTENTIAL_AM: moldyn->func3b_j1=albe_mult_3bp_j1; @@ -574,11 +495,16 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, t_3dvec orig; void *ptr; t_atom *atom; + char name[16]; new=a*b*c; 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; @@ -608,18 +534,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, case CUBIC: set_nn_dist(moldyn,lc); ret=cubic_init(a,b,c,lc,atom,&orig); + strcpy(name,"cubic"); break; case FCC: if(!origin) v3_scale(&orig,&orig,0.5); set_nn_dist(moldyn,0.5*sqrt(2.0)*lc); ret=fcc_init(a,b,c,lc,atom,&orig); + strcpy(name,"fcc"); break; case DIAMOND: if(!origin) v3_scale(&orig,&orig,0.25); set_nn_dist(moldyn,0.25*sqrt(3.0)*lc); ret=diamond_init(a,b,c,lc,atom,&orig); + strcpy(name,"diamond"); break; default: printf("unknown lattice type (%02x)\n",type); @@ -636,7 +565,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, } moldyn->count+=new; - printf("[moldyn] created lattice with %d atoms\n",new); + printf("[moldyn] created %s lattice with %d atoms\n",name,new); for(ret=0;retpt_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 */ @@ -1269,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 { @@ -1379,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 @@ -1389,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); @@ -1418,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])); @@ -1447,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 @@ -1459,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) { @@ -1468,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])); @@ -1635,7 +1586,6 @@ int moldyn_integrate(t_moldyn *moldyn) { atom=moldyn->atom; /* initialize linked cell method */ -printf(" hier soll es sein\n"); link_cell_init(moldyn,VERBOSE); /* logging & visualization */ @@ -1666,13 +1616,16 @@ printf(" hier soll es sein\n"); 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 */ - moldyn->time=0.0; - moldyn->total_steps=0; + // should have right values! + //moldyn->time=0.0; + //moldyn->total_steps=0; /* debugging, ignore */ moldyn->debug=0; @@ -1704,7 +1657,11 @@ printf(" hier soll es sein\n"); temperature_calc(moldyn); virial_sum(moldyn); pressure_calc(moldyn); - //thermodynamic_pressure_calc(moldyn); + /* + thermodynamic_pressure_calc(moldyn); + printf("\n\nDEBUG: numeric pressure calc: %f\n\n", + moldyn->tp/BAR); + */ /* calculate fluctuations + averages */ average_and_fluctuation_calc(moldyn); @@ -1785,6 +1742,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)); @@ -1854,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); @@ -1984,6 +1938,7 @@ int potential_force_calc(t_moldyn *moldyn) { do { jtom=this->current->data; +#endif if(jtom==&(itom[i])) continue; @@ -1995,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 @@ -2020,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++; @@ -2064,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++; @@ -2119,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++; @@ -2194,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; @@ -2205,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); } @@ -2408,18 +2369,13 @@ 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); - - /* only init link cell if it doesn't exist! */ - HIER WEITER - link_cell_init(moldyn,VERBOSE); - itom=moldyn->atom; for(i=0;icount;i++) { @@ -2437,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++; @@ -2679,8 +2635,6 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE) return 0; - d=sqrt(d); - /* now count this bonding ... */ ba=data; @@ -2887,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; +} +