From: hackbard Date: Tue, 28 Nov 2006 10:35:42 +0000 (+0000) Subject: more updates, now get the code running .... :/ X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=11bf064825efc474fa93c7270870856769a63e84;p=physik%2Fposic.git more updates, now get the code running .... :/ --- diff --git a/Makefile b/Makefile index 3c8fcd1..b037404 100644 --- a/Makefile +++ b/Makefile @@ -5,10 +5,10 @@ LDFLAGS=-lm OBJS=init/init.o visual/visual.o math/math.o random/random.o list/list.o OBJS+=moldyn.o -all: posic +all: sic -posic: $(OBJS) +sic: $(OBJS) .PHONY:clean clean: - rm -f *.o posic */*.o + rm -f *.o sic */*.o diff --git a/moldyn.c b/moldyn.c index fe54330..c0b99eb 100644 --- a/moldyn.c +++ b/moldyn.c @@ -23,206 +23,155 @@ #include "visual/visual.h" #include "list/list.h" -int moldyn_usage(char **argv) { - - printf("\n%s usage:\n\n",argv[0]); - printf("--- general options ---\n"); - printf("-E (log total energy)\n"); - printf("-M (log total momentum)\n"); - printf("-D (dump total information)\n"); - printf("-S (single save file)\n"); - printf("-V (rasmol file)\n"); - printf("--- physics options ---\n"); - printf("-T [K] (%f)\n",MOLDYN_TEMP); - printf("-t [s] (%.15f)\n",MOLDYN_TAU); - printf("-C [m] (%.15f)\n",MOLDYN_CUTOFF); - printf("-R (%d)\n",MOLDYN_RUNS); - printf(" -- integration algo --\n"); - printf(" -I (%d)\n",MOLDYN_INTEGRATE_DEFAULT); - printf(" 0: velocity verlet\n"); - printf(" -- potential --\n"); - printf(" -P \n"); - printf(" 0: harmonic oscillator\n"); - printf(" param1: spring constant\n"); - printf(" param2: equilibrium distance\n"); - printf(" 1: lennard jones\n"); - printf(" param1: epsilon\n"); - printf(" param2: sigma\n"); - printf("\n"); + +int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { + + //int ret; + + //ret=moldyn_parse_argv(moldyn,argc,argv); + //if(ret<0) return ret; + + memset(moldyn,0,sizeof(t_moldyn)); + + rand_init(&(moldyn->random),NULL,1); + moldyn->random.status|=RAND_STAT_VERBOSE; return 0; } -int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) { +int moldyn_shutdown(t_moldyn *moldyn) { - int i; + link_cell_shutdown(moldyn); + moldyn_log_shutdown(moldyn); + rand_close(&(moldyn->random)); + free(moldyn->atom); - memset(moldyn,0,sizeof(t_moldyn)); + return 0; +} - /* default values */ - moldyn->t=MOLDYN_TEMP; - moldyn->tau=MOLDYN_TAU; - moldyn->time_steps=MOLDYN_RUNS; - moldyn->integrate=velocity_verlet; - - /* parse argv */ - for(i=1;iewrite=atoi(argv[++i]); - strncpy(moldyn->efb,argv[++i],64); - break; - case 'M': - moldyn->mwrite=atoi(argv[++i]); - strncpy(moldyn->mfb,argv[++i],64); - break; - case 'S': - moldyn->swrite=atoi(argv[++i]); - strncpy(moldyn->sfb,argv[++i],64); - break; - case 'V': - moldyn->vwrite=atoi(argv[++i]); - strncpy(moldyn->vfb,argv[++i],64); - break; - case 'T': - moldyn->t=atof(argv[++i]); - break; - case 't': - moldyn->tau=atof(argv[++i]); - break; - case 'C': - moldyn->cutoff=atof(argv[++i]); - break; - case 'R': - moldyn->time_steps=atoi(argv[++i]); - break; - case 'I': - /* integration algorithm */ - switch(atoi(argv[++i])) { - case MOLDYN_INTEGRATE_VERLET: - moldyn->integrate=velocity_verlet; - break; - default: - printf("unknown integration algo %s\n",argv[i]); - moldyn_usage(argv); - return -1; - } +int set_int_alg(t_moldyn *moldyn,u8 algo) { - case 'P': - /* potential + params */ - switch(atoi(argv[++i])) { - case MOLDYN_POTENTIAL_HO: - hop.spring_constant=atof(argv[++i]); - hop.equilibrium_distance=atof(argv[++i]); - moldyn->pot_params=malloc(sizeof(t_ho_params)); - memcpy(moldyn->pot_params,&hop,sizeof(t_ho_params)); - moldyn->potential_force_function=harmonic_oscillator; - break; - case MOLDYN_POTENTIAL_LJ: - e=atof(argv[++i]); - s=atof(argv[++i]); - ljp.epsilon4=4*e; - ljp.sigma6=s*s*s*s*s*s; - ljp.sigma12=ljp.sigma6*ljp.sigma6; - moldyn->pot_params=malloc(sizeof(t_lj_params)); - memcpy(moldyn->pot_params,&ljp,sizeof(t_lj_params)); - moldyn->potential_force_function=lennard_jones; + switch(alg) { + case 'MOLDYN_INTEGRATE_VERLET': + moldyn->integrate=velocity_verlet; break; default: - printf("unknown potential %s\n",argv[i]); - moldyn_usage(argv); - return -1; - } - - default: - printf("unknown option %s\n",argv[i]); - moldyn_usage(argv); - return -1; - } - } else { - moldyn_usage(argv); + printf("unknown integration algorithm: %02x\",alg); return -1; - } } return 0; } -int moldyn_log_init(t_moldyn *moldyn) { +int set_cutoff(t_moldyn *moldyn,double cutoff) { - moldyn->lvstat=0; - t_visual *vis; + moldyn->cutoff=cutoff; - vis=&(moldyn->vis); + return 0; +} - if(moldyn->ewrite) { - moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC); - if(moldyn->efd<0) { - perror("[moldyn] efd open"); - return moldyn->efd; - } - dprintf(moldyn->efd,"# moldyn total energy logfile\n"); - moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E; - } +int set_temperature(t_moldyn *moldyn,double t) { + + moldyn->t=t; - if(moldyn->mwrite) { - moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC); - if(moldyn->mfd<0) { - perror("[moldyn] mfd open"); - return moldyn->mfd; - } - dprintf(moldyn->mfd,"# moldyn total momentum logfile\n"); - moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M; - } + return 0; +} - if(moldyn->swrite) - moldyn->lvstat|=MOLDYN_LVSTAT_SAVE; +int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) { - if((moldyn->vwrite)&&(vis)) { - moldyn->visual=vis; - visual_init(vis,moldyn->vfb); - moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL; - } + moldyn->dim.x=x; + moldyn->dim.y=y; + moldyn->dim.z=z; - moldyn->lvstat|=MOLDYN_LVSTAT_INITIALIZED; + if(visualize) { + moldyn->vis.dim.x=x; + moldyn->vis.dim.y=y; + moldyn->vis.dim.z=z; + } return 0; } -int moldyn_log_shutdown(t_moldyn *moldyn) { +int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) { - if(moldyn->efd) close(moldyn->efd); - if(moldyn->mfd) close(moldyn->efd); - if(moldyn->dfd) close(moldyn->efd); - if(moldyn->visual) visual_tini(moldyn->visual); + if(x) + moldyn->status|=MOLDYN_STAT_PBX; + + if(y) + moldyn->status|=MOLDYN_STAT_PBY; + + if(z) + moldyn->status|=MOLDYN_STAT_PBZ; return 0; } -int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { +int set_potential(t_moldyn *moldyn,u8 type,(int *)(func),void *params) { - int ret; + switch(type) { + case MOLDYN_1BP: + moldyn->pf_func1b=func; + moldyn->pot1b_params=params; + break; + case MOLDYN_2BP: + moldyn->pf_func2b=func; + moldyn->pot2b_params=params; + break; + case MOLDYN_3BP: + moldyn->pf_func3b=func; + moldyn->pot3b_params=params; + break; + default: + printf("unknown potential type: %02x\n",type); + return -1; + } - ret=moldyn_parse_argv(moldyn,argc,argv); - if(ret<0) return ret; + return 0; +} - ret=moldyn_log_init(moldyn); - if(ret<0) return ret; +int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer) { - rand_init(&(moldyn->random),NULL,1); - moldyn->random.status|=RAND_STAT_VERBOSE; - - moldyn->status=0; + switch(type) { + case LOG_TOTAL_ENERGY: + moldyn->ewrite=timer; + moldyn->efd=open(fb,O_WRONLY|O_CREAT|O_TRUNC); + if(moldyn->efd<0) { + perror("[moldyn] efd open"); + return moldyn->efd; + } + dprintf("# moldyn total energy log file\n"); + break; + case LOG_TOTAL_MOMENTUM: + moldyn->mwrite=timer; + moldyn->mfd=open(fb,O_WRONLY|O_CREAT|O_TRUNC); + if(moldyn->mfd<0) { + perror("[moldyn] mfd open"); + return moldyn->mfd; + } + dprintf("# moldyn total momentum log file\n"); + break; + case SAVE_STEP: + moldyn->swrite=timer; + strncpy(moldyn->sfb,fb,63); + break; + case VISUAL_STEP: + moldyn->mwrite=timer; + strncpy(moldyn->vfb,fb,63); + visual_init(&(moldyn->vis),fb); + break; + default: + printf("unknown log mechanism: %02x\n",type); + return -1; + } return 0; } -int moldyn_shutdown(t_moldyn *moldyn) { +int moldyn_log_shutdown(t_moldyn *moldyn) { - moldyn_log_shutdown(moldyn); - rand_close(&(moldyn->random)); - free(moldyn->atom); + if(moldyn->efd) close(moldyn->efd); + if(moldyn->mfd) close(moldyn->mfd); + if(moldyn->visual) visual_tini(moldyn->visual); return 0; } @@ -239,6 +188,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, atom=moldyn->atom; if(type==FCC) count*=4; + if(type==DIAMOND) count*=8; atom=malloc(count*sizeof(t_atom)); @@ -269,6 +219,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, return -1; } + moldyn->count=count; + while(count) { atom[count-1].element=element; atom[count-1].mass=mass; @@ -382,18 +334,18 @@ int scale_velocity(t_moldyn *moldyn) { return 0; } -double get_e_kin(t_atom *atom,int count) { +double get_e_kin(t_moldyn *moldyn) { int i; - double e; + t_atom *atom; - e=0.0; + atom=moldyn->atom; + moldyn->ekin=0.0; - for(i=0;icount;i++) + moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); - return e; + return moldyn->ekin; } double get_e_pot(t_moldyn *moldyn) { @@ -403,18 +355,16 @@ double get_e_pot(t_moldyn *moldyn) { double get_total_energy(t_moldyn *moldyn) { - double e; - - e=get_e_kin(moldyn->atom,moldyn->count); - e+=get_e_pot(moldyn); - - return e; + return(get_e_kin(moldyn)+get_e_pot(moldyn)); } -t_3dvec get_total_p(t_atom *atom, int count) { +t_3dvec get_total_p(t_moldyn *moldyn) { t_3dvec p,p_total; int i; + t_atom *atom; + + atom=moldyn->atom; v3_zero(&p_total); for(i=0;i