X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=mdrun.c;h=29d3160acdc9375864ff56788a293704e8c3ccee;hb=52485a28c5f9812c0f7a2d4e0094cd8000167b24;hp=8c5cf71c974f8816c58ee97805a9189d69ea7daa;hpb=0d1dfb1e5027d215fced8ca306dd658f692a2a44;p=physik%2Fposic.git diff --git a/mdrun.c b/mdrun.c index 8c5cf71..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,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)) @@ -332,7 +338,7 @@ int mdrun_parse_config(t_mdrun *mdrun) { csp.type|=CHSATTR_TCTRL; } if(!strncmp(word[i],"prelax",6)) { - csp.dp=atof(word[++i]); + csp.dp=atof(word[++i])*BAR; csp.type|=CHSATTR_PRELAX; } if(!strncmp(word[i],"trelax",6)) { @@ -368,11 +374,17 @@ 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 - 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; @@ -677,8 +762,8 @@ int mdrun_hook(void *ptr1,void *ptr2) { 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"); @@ -687,6 +772,11 @@ int mdrun_hook(void *ptr1,void *ptr2) { /* 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; @@ -696,7 +786,7 @@ int mdrun_hook(void *ptr1,void *ptr2) { } insert_atoms(moldyn,mdrun); iap->cnt_steps+=1; - break; + break; case STAGE_CONTINUE: stage_print(" -> continue\n\n"); if(stage->executed==TRUE) { @@ -705,7 +795,7 @@ int mdrun_hook(void *ptr1,void *ptr2) { } cp=stage->params; steps=cp->runs; - break; + break; case STAGE_ANNEAL: stage_print(" -> anneal\n\n"); ap=stage->params; @@ -713,7 +803,9 @@ int mdrun_hook(void *ptr1,void *ptr2) { 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: @@ -786,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; @@ -814,6 +921,14 @@ 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) { @@ -843,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)