X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=mdrun.c;h=29d3160acdc9375864ff56788a293704e8c3ccee;hb=959801961cfdd51e080e6500f51647b3397d336c;hp=63899a6a60e9ed72dfd07591792c486ccb6ef2a6;hpb=2dddad627c54f294c4c5b5d736178c3ef3288159;p=physik%2Fposic.git diff --git a/mdrun.c b/mdrun.c index 63899a6..29d3160 100644 --- a/mdrun.c +++ b/mdrun.c @@ -90,6 +90,9 @@ int add_stage(t_mdrun *mdrun,u8 type,void *params) { 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; @@ -138,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; @@ -185,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; } @@ -199,6 +203,8 @@ 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)) @@ -255,6 +261,7 @@ int mdrun_parse_config(t_mdrun *mdrun) { 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 @@ -274,7 +281,7 @@ int mdrun_parse_config(t_mdrun *mdrun) { break; } } - i=1; + i=2; if(cap.type&CHAATTR_REGION) { cap.x0=atof(word[1]); cap.y0=atof(word[2]); @@ -288,9 +295,8 @@ int mdrun_parse_config(t_mdrun *mdrun) { cap.element=atoi(word[i]); i+=1; } - wptr=word[i]; - for(i=0;ivisualize=atoi(word[1]); else if(!strncmp(word[0],"stage",5)) { // for every stage line, add a stage - if(!strncmp(word[1],"ins_atoms",9)) { + 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.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; @@ -480,13 +507,30 @@ 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) { 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; @@ -498,15 +542,32 @@ int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { 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: @@ -518,6 +579,12 @@ 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: %02x\n", ME,iap->type); @@ -528,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! @@ -550,8 +633,10 @@ 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; } @@ -630,9 +715,9 @@ int chsattr(t_moldyn *moldyn,t_mdrun *mdrun) { } if(csp->type&CHSATTR_AVGRST) { if(csp->avgrst) - mdrun->sattr|=CHSATTR_AVGRST; + mdrun->sattr|=SATTR_AVGRST; else - mdrun->sattr&=(~(CHSATTR_AVGRST)); + mdrun->sattr&=(~(SATTR_AVGRST)); } if(csp->type&CHSATTR_RSTEPS) { mdrun->relax_steps=csp->rsteps; @@ -674,12 +759,11 @@ int mdrun_hook(void *ptr1,void *ptr2) { /* default steps and tau values */ steps=mdrun->relax_steps; -printf("-------> %d\n",mdrun->relax_steps); tau=mdrun->timestep; /* check whether relaxation steps are necessary */ - if(!((check_pressure(moldyn,mdrun)==FALSE)|\ - (check_temperature(moldyn,mdrun)==FALSE))) { + if((check_pressure(moldyn,mdrun)==TRUE)&\ + (check_temperature(moldyn,mdrun)==TRUE)) { /* be verbose */ stage_print("\n###########################\n"); @@ -688,6 +772,11 @@ printf("-------> %d\n",mdrun->relax_steps); /* 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; @@ -697,7 +786,7 @@ printf("-------> %d\n",mdrun->relax_steps); } insert_atoms(moldyn,mdrun); iap->cnt_steps+=1; - break; + break; case STAGE_CONTINUE: stage_print(" -> continue\n\n"); if(stage->executed==TRUE) { @@ -706,7 +795,7 @@ printf("-------> %d\n",mdrun->relax_steps); } cp=stage->params; steps=cp->runs; - break; + break; case STAGE_ANNEAL: stage_print(" -> anneal\n\n"); ap=stage->params; @@ -714,7 +803,9 @@ printf("-------> %d\n",mdrun->relax_steps); 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: @@ -767,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) @@ -786,11 +878,25 @@ 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); - set_nn_dist(&moldyn,mdrun.nnd); if(set_potential(&moldyn,mdrun.potential)<0) return -1; switch(mdrun.potential) { @@ -815,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,1,mdrun.lx, + mdrun.m2,DEFAULT_ATOM_ATTR,1,mdrun.lx, mdrun.ly,mdrun.lz,&o); break; default: @@ -844,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)