X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=mdrun.c;h=29d3160acdc9375864ff56788a293704e8c3ccee;hb=959801961cfdd51e080e6500f51647b3397d336c;hp=6f7a2f02a0704491b5d8bf50762031e0149711ba;hpb=d7f67c88195ab155f2737e57cc5e81973d3feb0c;p=physik%2Fposic.git diff --git a/mdrun.c b/mdrun.c index 6f7a2f0..29d3160 100644 --- a/mdrun.c +++ b/mdrun.c @@ -62,12 +62,12 @@ int mdrun_parse_argv(t_mdrun *mdrun,int argc,char **argv) { return 0; } -del_stages(t_mdrun *mdrun) { +int del_stages(t_mdrun *mdrun) { t_list *sl; t_stage *stage; - sl=mdrun->stage; + sl=&(mdrun->stage); list_reset_f(sl); @@ -83,13 +83,16 @@ del_stages(t_mdrun *mdrun) { return 0; } -add_stage(t_mdrun *mdrun,u8 type,void *params) { +int add_stage(t_mdrun *mdrun,u8 type,void *params) { int psize; t_stage *stage; switch(type) { + case STAGE_DISPLACE_ATOM: + psize=sizeof(t_displace_atom_params); + break; case STAGE_INSERT_ATOMS: psize=sizeof(t_insert_atoms_params); break; @@ -106,7 +109,8 @@ add_stage(t_mdrun *mdrun,u8 type,void *params) { psize=sizeof(t_chsattr_params); break; default: - + printf("%s unknown stage type: %02x\n",ME,type); + return -1; } stage=malloc(sizeof(t_stage)); @@ -126,7 +130,7 @@ add_stage(t_mdrun *mdrun,u8 type,void *params) { memcpy(stage->params,params,psize); - list_add_immediate_f(mdrun->stage,stage); + list_add_immediate_f(&(mdrun->stage),stage); return 0; } @@ -137,10 +141,11 @@ int mdrun_parse_config(t_mdrun *mdrun) { char error[128]; char line[128]; char *wptr; - char word[16][32]; + char word[16][64]; int wcnt; - int i; + int i,o; + t_displace_atom_params dap; t_insert_atoms_params iap; t_continue_params cp; t_anneal_params ap; @@ -168,6 +173,13 @@ int mdrun_parse_config(t_mdrun *mdrun) { if((line[0]=='#')|(ret==1)) continue; + // reset + memset(&iap,0,sizeof(t_insert_atoms_params)); + memset(&cp,0,sizeof(t_continue_params)); + memset(&ap,0,sizeof(t_anneal_params)); + memset(&cap,0,sizeof(t_chaattr_params)); + memset(&csp,0,sizeof(t_chsattr_params)); + // get command + args wcnt=0; while(1) { @@ -177,7 +189,7 @@ int mdrun_parse_config(t_mdrun *mdrun) { wptr=strtok(line," \t"); if(wptr==NULL) break; - strncpy(word[wcnt],wptr,32); + strncpy(word[wcnt],wptr,64); wcnt+=1; } @@ -191,8 +203,12 @@ int mdrun_parse_config(t_mdrun *mdrun) { if(!strncmp(word[1],"lj",2)) mdrun->potential=MOLDYN_POTENTIAL_LJ; } + else if(!strncmp(word[0],"continue",8)) + strncpy(mdrun->continue_file,word[1],128); else if(!strncmp(word[0],"cutoff",6)) mdrun->cutoff=atof(word[1]); + else if(!strncmp(word[0],"nnd",3)) + mdrun->nnd=atof(word[1]); else if(!strncmp(word[0],"intalgo",7)) { if(!strncmp(word[1],"verlet",5)) mdrun->intalgo=MOLDYN_INTEGRATE_VERLET; @@ -234,51 +250,53 @@ int mdrun_parse_config(t_mdrun *mdrun) { } else if(!strncmp(word[0],"element1",8)) { mdrun->element1=atoi(word[1]); - switch(mdrun->element1) { - case SI: - mdrun->m1=M_SI; - break; - case C: - mdrun->m1=M_C; - break; - default: - printf("%s unknown element1: %s|%d\n", - ME,word[1],mdrun->element1); - return -1; - } + mdrun->m1=pse_mass[mdrun->element1]; } else if(!strncmp(word[0],"element2",8)) { mdrun->element2=atoi(word[1]); - switch(mdrun->element2) { - case SI: - mdrun->m2=M_SI; - break; - case C: - mdrun->m2=M_C; - break; - default: - printf("%s unknown element2: %s|%d\n", - ME,word[1],mdrun->element2); - return -1; - } + mdrun->m2=pse_mass[mdrun->element2]; } else if(!strncmp(word[0],"fill",6)) { // only lc mode by now mdrun->lx=atoi(word[2]); mdrun->ly=atoi(word[3]); mdrun->lz=atoi(word[4]); + mdrun->lc=atof(word[5]); } else if(!strncmp(word[0],"aattr",5)) { // for aatrib line we need a special stage // containing one schedule of 0 loops ... - if(!strncmp(word[1],"all",3)) { - cap.type= - - - HIER WEITER + for(i=0;ip_tau=atof(word[1]); - mdrun->t_tau=atof(word[2]); + else if(!strncmp(word[0],"sattr",5)) { + // for satrib line we need a special stage + // containing one schedule of 0 loops ... + csp.type=0; + for(i=1;iprerun=atoi(word[1]); @@ -328,10 +374,92 @@ int mdrun_parse_config(t_mdrun *mdrun) { mdrun->visualize=atoi(word[1]); else if(!strncmp(word[0],"stage",5)) { // for every stage line, add a stage - printf("%s stage: %s\n",ME,word[1]); + if(!strncmp(word[1],"displace",8)) { + dap.nr=atoi(word[2]); + dap.dx=atof(word[3]); + dap.dy=atof(word[4]); + dap.dz=atof(word[5]); + add_stage(mdrun,STAGE_DISPLACE_ATOM,&dap); + } + else if(!strncmp(word[1],"ins_atoms",9)) { + iap.ins_steps=atoi(word[2]); + iap.ins_atoms=atoi(word[3]); + iap.element=atoi(word[4]); + iap.brand=atoi(word[5]); + for(i=0;isattr&SATTR_PRELAX)) + if(!(mdrun->sattr&SATTR_PRELAX)) { return TRUE; + } delta=moldyn->p_avg-moldyn->p_ref; @@ -378,33 +507,67 @@ int check_temperature(t_moldyn *moldyn,t_mdrun *mdrun) { return TRUE; } +int displace_atom(t_moldyn *moldyn,t_mdrun *mdrun) { + + t_displace_atom_params *dap; + t_stage *stage; + t_atom *atom; + + stage=mdrun->stage.current->data; + dap=stage->params; + + atom=&(moldyn->atom[dap->nr]); + atom->r.x+=dap->dx; + atom->r.y+=dap->dy; + atom->r.z+=dap->dz; + + return 0; +} + int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { - insert_atoms_params *iap; + t_insert_atoms_params *iap; t_stage *stage; t_atom *atom; t_3dvec r,v,dist; - double d,dmin; + double d,dmin,o; int cnt,i; double x,y,z; double x0,y0,z0; u8 cr_check,run; - stage=mdrun->stage->current->data; + stage=mdrun->stage.current->data; iap=stage->params; cr_check=FALSE; v.x=0.0; v.y=0.0; v.z=0.0; + x=0; y=0; z=0; + o=0; dmin=0; + + switch(mdrun->lattice) { + case CUBIC: + o=0.5*mdrun->lc; + break; + case FCC: + o=0.25*mdrun->lc; + break; + case DIAMOND: + case ZINCBLENDE: + o=0.125*mdrun->lc; + break; + default: + break; + } switch(iap->type) { case INS_TOTAL: x=moldyn->dim.x; - x0=0.0; + x0=-x/2.0; y=moldyn->dim.y; - y0=0.0; + y0=-y/2.0; z=moldyn->dim.z; - z0=0.0; + z0=-z/2.0; cr_check=TRUE; break; case INS_REGION: @@ -416,8 +579,15 @@ int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { z0=iap->z0; cr_check=TRUE; break; + case INS_POS: + x0=iap->x0; + y0=iap->y0; + z0=iap->z0; + cr_check=FALSE; + break; default: - printf("%s unknown insertion mode\n"); + printf("%s unknown insertion mode: %02x\n", + ME,iap->type); return -1; } @@ -425,9 +595,25 @@ int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { while(cntins_atoms) { run=1; while(run) { - r.x=(rand_get_double(&(moldyn->random))-0.5)*x+x0; - r.y=(rand_get_double(&(moldyn->random))-0.5)*y+y0; - r.z=(rand_get_double(&(moldyn->random))-0.5)*z+z0; + if(iap->type!=INS_POS) { + r.x=rand_get_double(&(moldyn->random))*x; + r.y=rand_get_double(&(moldyn->random))*y; + r.z=rand_get_double(&(moldyn->random))*z; + } + else { + r.x=0.0; + r.y=0.0; + r.z=0.0; + } + r.x+=x0; + r.y+=y0; + r.z+=z0; + // offset + if(iap->type!=INS_TOTAL) { + r.x+=o; + r.y+=o; + r.z+=o; + } run=0; if(cr_check==TRUE) { dmin=1000; // for sure too high! @@ -447,22 +633,24 @@ int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { } add_atom(moldyn,iap->element,pse_mass[iap->element], iap->brand,iap->attr,&r,&v); - printf("%s atom inserted: %f %f %f | d squared = %f\n", - ME,r.x,r.y,r.z,dmin); + printf("%s atom inserted (%d/%d): %f %f %f\n", + ME,(iap->cnt_steps+1)*iap->ins_atoms, + iap->ins_steps*iap->ins_atoms,r.x,r.y,r.z); + printf(" -> d2 = %f/%f\n",dmin,iap->cr*iap->cr); cnt+=1; } return 0; } -int chaatrib(t_moldyn *moldyn,t_mdrun *mdrun) { +int chaatr(t_moldyn *moldyn,t_mdrun *mdrun) { t_stage *stage; t_chaattr_params *cap; t_atom *atom; - u8 assigne; + int i; - stage=mdrun->stage->current->data; + stage=mdrun->stage.current->data; cap=stage->params; for(i=0;icount;i++) { @@ -496,41 +684,54 @@ int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) { t_stage *stage; t_chsattr_params *csp; - stage=mdrun->stage->current->data; + stage=mdrun->stage.current->data; csp=stage->params; - switch(csp->type) { - case CHSATTR_PCTRL: - set_p_scale(moldyn,csp->ctrl_type,csp->tau); - break; - case CHSATTR_TCTRL: - set_t_scale(moldyn,csp->ctrl_type,csp->tau); - break; - case CHSATTR_PRELAX: - mdrun->sattr&=(^(SATTR_PRELAX)); - if(csp->ctrl==TRUE) - mdrun->sattr|=SATTR_PRELAX; - break; - case CHSATTR_TRELAX: - mdrun->sattr&=(^(SATTR_TRELAX)); - if(csp->ctrl==TRUE) - mdrun->sattr|=SATTR_TRELAX; - break; - case CHSATTR_AVGRST: - mdrun->sattr&=(^(SATTR_AVGRST)); - if(csp->ctrl==TRUE) - mdrun->sattr|=SATTR_AVGRST; - break; - default: - printf("%s unknown system attribute change\n",ME); - return -1; + if(csp->type&CHSATTR_PCTRL) { + if(csp->ptau>0) + set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau); + else + set_p_scale(moldyn,P_SCALE_BERENDSEN,csp->ptau); + } + if(csp->type&CHSATTR_TCTRL) { + if(csp->ttau>0) + set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau); + else + set_t_scale(moldyn,T_SCALE_BERENDSEN,csp->ttau); + } + if(csp->type&CHSATTR_PRELAX) { + if(csp->dp<0) + mdrun->sattr&=(~(SATTR_PRELAX)); + else + mdrun->sattr|=SATTR_PRELAX; + mdrun->dp=csp->dp; + } + if(csp->type&CHSATTR_TRELAX) { + if(csp->dt<0) + mdrun->sattr&=(~(SATTR_TRELAX)); + else + mdrun->sattr|=SATTR_TRELAX; + mdrun->dt=csp->dt; + } + if(csp->type&CHSATTR_AVGRST) { + if(csp->avgrst) + mdrun->sattr|=SATTR_AVGRST; + else + mdrun->sattr&=(~(SATTR_AVGRST)); + } + if(csp->type&CHSATTR_RSTEPS) { + mdrun->relax_steps=csp->rsteps; } return 0; } -int mdrun_hook(t_moldyn *moldyn,void *ptr) { +#define stage_print(m) if(!(stage->executed)) \ + printf("%s",m) +int mdrun_hook(void *ptr1,void *ptr2) { + + t_moldyn *moldyn; t_mdrun *mdrun; t_stage *stage; t_list *sl; @@ -542,8 +743,10 @@ int mdrun_hook(t_moldyn *moldyn,void *ptr) { t_continue_params *cp; t_anneal_params *ap; - mdrun=ptr; - sl=mdrun->stage; + moldyn=ptr1; + mdrun=ptr2; + + sl=&(mdrun->stage); change_stage=FALSE; @@ -552,18 +755,30 @@ int mdrun_hook(t_moldyn *moldyn,void *ptr) { return 0; /* get stage description */ - stage=mdrun->sl->current->data; + stage=sl->current->data; /* default steps and tau values */ steps=mdrun->relax_steps; tau=mdrun->timestep; /* check whether relaxation steps are necessary */ - if(!((check_pressure==FALSE)|(check_temperature==FALSE))) { + if((check_pressure(moldyn,mdrun)==TRUE)&\ + (check_temperature(moldyn,mdrun)==TRUE)) { + + /* be verbose */ + stage_print("\n###########################\n"); + stage_print("# [mdrun] executing stage #\n"); + stage_print("###########################\n\n"); /* stage specific stuff */ switch(stage->type) { + case STAGE_DISPLACE_ATOM: + stage_print(" -> displace atom\n\n"); + displace_atom(moldyn,mdrun); + change_stage=TRUE; + break; case STAGE_INSERT_ATOMS: + stage_print(" -> insert atoms\n\n"); iap=stage->params; if(iap->cnt_steps==iap->ins_steps) { change_stage=TRUE; @@ -571,30 +786,36 @@ int mdrun_hook(t_moldyn *moldyn,void *ptr) { } insert_atoms(moldyn,mdrun); iap->cnt_steps+=1; - break; + break; case STAGE_CONTINUE: + stage_print(" -> continue\n\n"); if(stage->executed==TRUE) { change_stage=TRUE; break; } cp=stage->params; steps=cp->runs; - break; + break; case STAGE_ANNEAL: + stage_print(" -> anneal\n\n"); ap=stage->params; if(ap->count==ap->runs) { change_stage=TRUE; break; } - set_temperature(moldyn,moldyn->t_ref+ap->dt); + if(moldyn->t_ref+ap->dt>=0.0) + set_temperature(moldyn, + moldyn->t_ref+ap->dt); ap->count+=1; break; case STAGE_CHAATTR: - chaatrib(moldyn,mdrun); + stage_print(" -> chaattr\n\n"); + chaatr(moldyn,mdrun); change_stage=TRUE; break; case STAGE_CHSATTR: - chsatrib(moldyn,mdrun); + stage_print(" -> chsattr\n\n"); + chsattr(moldyn,mdrun); change_stage=TRUE; break; default: @@ -637,8 +858,9 @@ int main(int argc,char **argv) { t_moldyn moldyn; t_3dvec o; - /* clear mdrun struct */ + /* clear structs */ memset(&mdrun,0,sizeof(t_mdrun)); + memset(&moldyn,0,sizeof(t_moldyn)); /* parse arguments */ if(mdrun_parse_argv(&mdrun,argc,argv)<0) @@ -656,9 +878,24 @@ int main(int argc,char **argv) { /* sanity check! */ /* prepare simulation */ - moldyn_init(&moldyn,argc,argv); + + if(mdrun.continue_file[0]) { + // read the save file + moldyn_read_save_file(&moldyn,mdrun.continue_file); + // manualaadjusting some stuff + moldyn.argc=argc; + moldyn.args=argv; + rand_init(&(moldyn.random),NULL,1); + moldyn.random.status|=RAND_STAT_VERBOSE; + } + else { + moldyn_init(&moldyn,argc,argv); + } + if(set_int_alg(&moldyn,mdrun.intalgo)<0) return -1; + + /* potential */ set_cutoff(&moldyn,mdrun.cutoff); if(set_potential(&moldyn,mdrun.potential)<0) return -1; @@ -684,27 +921,35 @@ int main(int argc,char **argv) { ME,mdrun.potential); return -1; } + + /* if it is a continue run, reset schedule and skip lattice init */ + if(mdrun.continue_file[0]) { + memset(&(moldyn.schedule),0,sizeof(t_moldyn_schedule)); + goto addsched; + } + + /* initial lattice and dimensions */ set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis); set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz); switch(mdrun.lattice) { case FCC: create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1, - mdrun.m1,0,0,mdrun.lx, + mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx, mdrun.ly,mdrun.lz,NULL); break; case DIAMOND: create_lattice(&moldyn,DIAMOND,mdrun.lc,mdrun.element1, - mdrun.m1,0,0,mdrun.lx, + mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx, mdrun.ly,mdrun.lz,NULL); break; case ZINCBLENDE: o.x=0.5*0.25*mdrun.lc; o.y=o.x; o.z=o.x; create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1, - mdrun.m1,0,0,mdrun.lx, + mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx, mdrun.ly,mdrun.lz,&o); o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x; create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2, - mdrun.m2,0,0,mdrun.lx, + mdrun.m2,DEFAULT_ATOM_ATTR,1,mdrun.lx, mdrun.ly,mdrun.lz,&o); break; default: @@ -713,10 +958,17 @@ int main(int argc,char **argv) { return -1; } moldyn_bc_check(&moldyn); + + /* temperature and pressure */ set_temperature(&moldyn,mdrun.temperature); set_pressure(&moldyn,mdrun.pressure); thermal_init(&moldyn,TRUE); + +addsched: + /* first schedule */ moldyn_add_schedule(&moldyn,mdrun.prerun,mdrun.timestep); + + /* log */ moldyn_set_log_dir(&moldyn,mdrun.sdir); moldyn_set_report(&moldyn,"CHANGE ME","CHANGE ME TOO"); if(mdrun.elog) @@ -734,12 +986,6 @@ int main(int argc,char **argv) { moldyn_set_log(&moldyn,CREATE_REPORT,0); set_avg_skip(&moldyn,mdrun.avgskip); - /* pt scaling */ - if(mdrun.p_tau!=0) - set_p_scale(&moldyn,P_SCALE_BERENDSEN,mdrun.p_tau); - if(mdrun.t_tau!=0) - set_t_scale(&moldyn,T_SCALE_BERENDSEN,mdrun.t_tau); - /* prepare the hook function */ moldyn_set_schedule_hook(&moldyn,&mdrun_hook,&mdrun);