From: hackbard Date: Wed, 6 Feb 2008 00:49:32 +0000 (+0100) Subject: float-store bugfix (unstatisfying!) + slightly new design of sic.c X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e1080fc0dd66b0cf5b7715c5e99e7a34ac04a8cf;p=physik%2Fposic.git float-store bugfix (unstatisfying!) + slightly new design of sic.c --- diff --git a/Makefile b/Makefile index 28fd2b1..83a14b6 100644 --- a/Makefile +++ b/Makefile @@ -13,7 +13,8 @@ CFLAGS += -DDSTART=400 -DDEND=600 #CFLAGS += -DVDEBUG #CFLAGS += -DTERSOFF_ORIG -LDFLAGS=-lm +LDFLAGS = -lm +#LDFLAGS += -lefence DEPS = moldyn.o random/random.o list/list.o DEPS += potentials/lennard_jones.o potentials/harmonic_oscillator.o @@ -23,7 +24,7 @@ all: posic sic fluctuation_calc postproc pair_correlation_calc posic: $(DEPS) -sic: $(DEPS) +sic: $(DEPS) config.h postproc: $(DEPS) diff --git a/moldyn.c b/moldyn.c index 5681779..edda007 100644 --- a/moldyn.c +++ b/moldyn.c @@ -956,10 +956,12 @@ double pressure_calc(t_moldyn *moldyn) { int average_and_fluctuation_calc(t_moldyn *moldyn) { + int denom; + if(moldyn->total_stepsavg_skip) return 0; - int denom=moldyn->total_steps+1-moldyn->avg_skip; + denom=moldyn->total_steps+1-moldyn->avg_skip; /* assume up to date energies, temperature, pressure etc */ @@ -1031,7 +1033,8 @@ printf(" --> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->coun double thermodynamic_pressure_calc(t_moldyn *moldyn) { - t_3dvec dim,*tp; + t_3dvec dim; + //t_3dvec *tp; double u_up,u_down,dv; double scale,p; t_atom *store; @@ -1436,7 +1439,6 @@ int link_cell_shutdown(t_moldyn *moldyn) { for(i=0;icells;i++) { #ifdef STATIC_LISTS - //printf(" ---> %d free %p\n",i,lc->subcell[i]); free(lc->subcell[i]); #else //printf(" ---> %d free %p\n",i,lc->subcell[i].start); @@ -1650,13 +1652,14 @@ int moldyn_integrate(t_moldyn *moldyn) { /* get current time */ gettimeofday(&t2,NULL); - printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)", - sched->count,i, - moldyn->t,moldyn->t_avg, - moldyn->p_avg/BAR,moldyn->gp_avg/BAR, - moldyn->volume, - (int)(t2.tv_sec-t1.tv_sec)); - fflush(stdout); +printf("\rsched:%d, steps:%d/%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)", + sched->count,i,moldyn->total_steps, + moldyn->t,moldyn->t_avg, + moldyn->p_avg/BAR,moldyn->gp_avg/BAR, + moldyn->volume, + (int)(t2.tv_sec-t1.tv_sec)); + + fflush(stdout); /* copy over time */ t1=t2; @@ -1670,8 +1673,8 @@ int moldyn_integrate(t_moldyn *moldyn) { /* check for hooks */ if(sched->hook) { - printf("\n ## schedule hook %d/%d start ##\n", - sched->count+1,sched->total_sched-1); + printf("\n ## schedule hook %d start ##\n", + sched->count); sched->hook(moldyn,sched->hook_params); printf(" ## schedule hook end ##\n"); } @@ -2182,10 +2185,14 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) { } size=sizeof(t_moldyn); - cnt=read(fd,moldyn,size); - if(cnt!=size) { - perror("[moldyn] load save file read (moldyn)"); - return cnt; + + while(size) { + cnt=read(fd,moldyn,size); + if(cnt<0) { + perror("[moldyn] load save file read (moldyn)"); + return cnt; + } + size-=cnt; } size=moldyn->count*sizeof(t_atom); @@ -2196,10 +2203,13 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) { return -1; } - cnt=read(fd,moldyn->atom,size); - if(cnt!=size) { - perror("[moldyn] load save file read (atoms)"); - return cnt; + while(size) { + cnt=read(fd,moldyn->atom,size); + if(cnt<0) { + perror("[moldyn] load save file read (atoms)"); + return cnt; + } + size-=cnt; } // hooks etc ... @@ -2258,13 +2268,14 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { t_list *this; unsigned char bc; t_3dvec dist; - double d,norm; + double d; + //double norm; int o,s; unsigned char ibrand; lc=&(moldyn->lc); - slots=(int)(moldyn->cutoff/dr); + slots=moldyn->cutoff/dr; o=2*slots; printf("[moldyn] pair correlation calc info:\n"); @@ -2323,12 +2334,9 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { jtom=this->current->data; #endif - - if(jtom==&(itom[i])) - continue; - - /* only count pairs once */ - if(itom[i].tag>jtom->tag) + /* only count pairs once, + * skip same atoms */ + if(itom[i].tag>=jtom->tag) continue; /* @@ -2340,14 +2348,18 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { if(bc) check_per_bound(moldyn,&dist); d=v3_absolute_square(&dist); - /* ignore if greater cutoff */ - if(d>moldyn->cutoff_square) + /* ignore if greater or equal cutoff */ + if(d>=moldyn->cutoff_square) continue; /* fill the slots */ d=sqrt(d); s=(int)(d/dr); + /* should never happen but it does 8) - + * related to -ffloat-store problem! */ + if(s>=slots) s=slots-1; + if(ibrand!=jtom->brand) { /* mixed */ stat[s]+=1; @@ -2360,7 +2372,6 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { /* type b - type b bonds */ stat[s+o]+=1; } - #ifdef STATIC_LISTS } #else diff --git a/moldyn.h b/moldyn.h index ee45228..2290743 100644 --- a/moldyn.h +++ b/moldyn.h @@ -70,7 +70,7 @@ typedef struct s_linkcell { int dnlc; /* direct neighbour lists counter */ } t_linkcell; -#define MAX_ATOMS_PER_LIST 10 +#define MAX_ATOMS_PER_LIST 20 /* moldyn schedule structure */ typedef struct s_moldyn_schedule { diff --git a/pair_correlation_calc.c b/pair_correlation_calc.c index f193614..f6a256f 100644 --- a/pair_correlation_calc.c +++ b/pair_correlation_calc.c @@ -47,8 +47,15 @@ int main(int argc,char **argv) { return ret; } + //moldyn.cutoff*=2; + //moldyn.cutoff_square*=4; + dr=atof(argv[2]); - slots=(int)(moldyn.cutoff/dr); + slots=moldyn.cutoff/dr; + printf("[pair corr calc]\n"); + printf(" slots: %d\n",slots); + printf(" cutoff: %f\n",moldyn.cutoff); + printf(" dr: %f\n",dr); stat=(double *)malloc(3*slots*sizeof(double)); if(stat==NULL) { diff --git a/posic.c b/posic.c index eb2db91..6b92a5e 100644 --- a/posic.c +++ b/posic.c @@ -34,6 +34,7 @@ int main(int argc,char **argv) { tp=NULL; ap=NULL; + memset(&md,0,sizeof(t_moldyn)); return 0; diff --git a/ppm2avi b/ppm2avi index 3680b80..a669290 100755 --- a/ppm2avi +++ b/ppm2avi @@ -5,4 +5,5 @@ OUTFILE=md.avi [ -n "$2" ] && OUTFILE=$2 -mencoder mf://$1/*.png -ovc copy -o $1/$OUTFILE >/dev/null 2>&1 +#mencoder mf://$1/*.png -ovc copy -o $1/$OUTFILE >/dev/null 2>&1 +mencoder mf://$1/*.png -ovc xvid -xvidencopts bitrate=1000 -o $1/$OUTFILE >/dev/null 2>&1 diff --git a/run b/run index 34b75bd..36b1b75 100755 --- a/run +++ b/run @@ -5,8 +5,10 @@ if [ "$?" == "0" ]; then #./perms if [ "$1" ] ; then #./visualize -w 640 -h 480 -d $1 - ./visualize -w 640 -h 480 -d $1 -nll -2.4 -2.4 -2.4 -fur 3.8 3.8 3.8 -b -2.03 -2.03 -2.03 3.39 3.39 3.39 -r 0.6 -c 1.5 -15.0 1.5 -B 0.1 + #./visualize -w 640 -h 480 -d $1 -nll -2.4 -2.4 -2.4 -fur 3.8 3.8 3.8 -b -2.03 -2.03 -2.03 3.39 3.39 3.39 -r 0.6 -c 1.5 -15.0 1.5 -B 0.1 #rasmol -32 -nodisplay < $1/visualize.scr > /dev/null 2>&1 - ./ppm2avi $1 + #./ppm2avi $1 + cp -v sic.c $1/sic.c + cp -v config.h $1/config.h fi fi diff --git a/sic.c b/sic.c index 05de728..e1103bd 100644 --- a/sic.c +++ b/sic.c @@ -13,101 +13,42 @@ #include "potentials/harmonic_oscillator.h" #include "potentials/lennard_jones.h" #include "potentials/albe.h" - #ifdef TERSOFF_ORIG #include "potentials/tersoff_orig.h" #else #include "potentials/tersoff.h" #endif -//#define INJECT 800 -#define INJECT 1 -#define NR_ATOMS 1 -#define R_C 1.5 -#define T_C 5.0 -//#define INJ_LENX (1*ALBE_LC_SIC) -//#define INJ_LENY (1*ALBE_LC_SIC) -//#define INJ_LENZ (1*ALBE_LC_SIC) -#define INJ_LENX (1*ALBE_LC_SI) -#define INJ_LENY (1*ALBE_LC_SI) -#define INJ_LENZ (1*ALBE_LC_SI) -#define INJ_TYPE_SILICON -//#define INJ_TYPE_CARBON -#define INJ_OFFSET (ALBE_LC_SI/8.0) -#define RELAX_S 20 - -#define LCNTX 9 -#define LCNTY 9 -#define LCNTZ 9 -#define PRERUN 40 -#define POSTRUN 3000 - -#define R_TITLE "Silicon self-interstitial" -#define LOG_E 10 -#define LOG_T 10 -#define LOG_P 10 -#define LOG_S 100 -#define LOG_V 20 - typedef struct s_hp { - int a_count; /* atom count */ - u8 quit; /* quit mark */ - int argc; /* arg count */ - char **argv; /* args */ + int prerun_count; /* prerun count */ + int insert_count; /* insert count */ + int postrun_count; /* post run count */ + unsigned char state; /* current state */ + int argc; /* arg count */ + char **argv; /* args */ } t_hp; -int hook_del_atom(void *moldyn,void *hook_params) { - - t_moldyn *md; - t_hp *hp; +#define STATE_PRERUN 0x00 +#define STATE_INSERT 0x01 +#define STATE_POSTRUN 0x02 - md=moldyn; - hp=hook_params; +/* include the config file */ +#include "config.h" - set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0); - del_atom(md,2); +int insert_atoms(t_moldyn *moldyn) { - return 0; -} - -int hook_add_atom(void *moldyn,void *hook_params) { - - t_moldyn *md; + int i,j; + u8 run; t_3dvec r,v,dist; double d; - unsigned char run; - int i,j; - t_atom *atom; - t_hp *hp; - - md=moldyn; - hp=hook_params; - /* quit */ - if(hp->quit) - return 0; - - /* switch on t scaling */ - if(md->schedule.count==0) - set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0); + t_atom *atom; - /* last schedule add if there is enough carbon inside */ - if(hp->a_count==(INJECT*NR_ATOMS)) { - hp->quit=1; - moldyn_add_schedule(md,POSTRUN,1.0); - return 0; - } + atom=moldyn->atom; - /* more relaxing time for too high temperatures */ - if(md->t-md->t_ref>T_C) { - moldyn_add_schedule(md,RELAX_S,1.0); - return 0; - } + v.x=0; v.y=0; v.z=0; - /* inject carbon atoms */ - printf("injecting another %d atoms ... (-> %d / %d)\n", - NR_ATOMS,hp->a_count+NR_ATOMS,INJECT*NR_ATOMS); - for(j=0;jatom[4372].r.y=(-0.5+0.125+0.125)*ALBE_LC_SI; */ // random - /* - r.x=(rand_get_double(&(md->random))-0.5)*INJ_LENX; - r.y=(rand_get_double(&(md->random))-0.5)*INJ_LENY; - r.z=(rand_get_double(&(md->random))-0.5)*INJ_LENZ; - */ + // + r.x=(rand_get_double(&(moldyn->random))-0.5)*INS_LENX; + r.y=(rand_get_double(&(moldyn->random))-0.5)*INS_LENY; + r.z=(rand_get_double(&(moldyn->random))-0.5)*INS_LENZ; + // // offset - r.x+=INJ_OFFSET; - r.y+=INJ_OFFSET; - r.z+=INJ_OFFSET; + r.x+=INS_OFFSET; + r.y+=INS_OFFSET; + r.z+=INS_OFFSET; /* assume valid coordinates */ run=0; - for(i=0;icount;i++) { - atom=&(md->atom[i]); + for(i=0;icount;i++) { + atom=&(moldyn->atom[i]); v3_sub(&dist,&(atom->r),&r); + check_per_bound(moldyn,&dist); d=v3_absolute_square(&dist); /* reject coordinates */ - if(da_count+=NR_ATOMS; - /* add schedule for simulating injected atoms ;) */ - moldyn_add_schedule(md,RELAX_S,1.0); + return 0; +} + +int sic_hook(void *moldyn,void *hook_params) { + + t_hp *hp; + t_moldyn *md; + int steps; + double tau; + + hp=hook_params; + md=moldyn; + + /* switch on t scaling */ + if(md->schedule.count==0) + set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0); + + /* my lousy state machine ! */ + + /* switch to insert state immediately */ + if(hp->state==STATE_PRERUN) + hp->state=STATE_INSERT; + + /* act according to state */ + switch(hp->state) { + case STATE_INSERT: + /* check temperature */ + if(md->t_avg-md->t_ref>INS_DELTA_TC) { + steps=INS_RELAX; + tau=INS_TAU; + break; + } + /* insert atoms */ + hp->insert_count+=1; + printf(" ### insert atoms (%d/%d) ###\n", + hp->insert_count*INS_ATOMS,INS_RUNS*INS_ATOMS); + insert_atoms(md); + /* change state after last insertion */ + if(hp->insert_count==INS_RUNS) + hp->state=STATE_POSTRUN; + break; + case STATE_POSTRUN: + /* settings */ + if(md->t-md->t_ref>POST_DELTA_TC) { + steps=POST_RELAX; + tau=POST_TAU; + } + /* decrease temperature */ + hp->postrun_count+=1; + printf(" ### postrun (%d/%d) ###\n", + hp->postrun_count,POST_RUNS); + set_temperature(md,md->t_ref-POST_DT); + if(hp->postrun_count==POST_RUNS) + return 0; + break; + default: + printf("[hook] FATAL (default case!?!)\n"); + break; + } + + /* add schedule */ + moldyn_add_schedule(md,steps,tau); return 0; } int main(int argc,char **argv) { - /* check argv */ - //if(argc!=3) { - // printf("[sic] usage: %s \n",argv[0]); - // return -1; - //} - /* main moldyn structure */ t_moldyn md; @@ -333,7 +325,6 @@ int main(int argc,char **argv) { // ATOM_ATTR_2BP|ATOM_ATTR_HB, 0,LCNTX,LCNTY,LCNTZ,NULL); // 1,LCNTX,LCNTY,LCNTZ,NULL); - // /* create zinkblende structure */ /* @@ -398,7 +389,7 @@ int main(int argc,char **argv) { set_pressure(&md,BAR); /* set amount of steps to skip before average calc */ - set_avg_skip(&md,(8.0/10.0*PRERUN)); + set_avg_skip(&md,AVG_SKIP); /* set p/t scaling */ //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0); @@ -411,13 +402,13 @@ int main(int argc,char **argv) { thermal_init(&md,TRUE); /* create the simulation schedule */ - moldyn_add_schedule(&md,PRERUN,1.0); + moldyn_add_schedule(&md,PRERUN,PRE_TAU); /* schedule hook function */ memset(&hookparam,0,sizeof(t_hp)); hookparam.argc=argc; hookparam.argv=argv; - moldyn_set_schedule_hook(&md,&hook_add_atom,&hookparam); + moldyn_set_schedule_hook(&md,&sic_hook,&hookparam); //moldyn_set_schedule_hook(&md,&hook_del_atom,&hookparam); //moldyn_add_schedule(&md,POSTRUN,1.0);