X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=mdrun.c;h=34807e76bd76f646f5d587705d09ce0e05840614;hb=1965279af348fbb9b69459d01516bbcd52b6f240;hp=1368880b7c62cde2eb66761bf726c962751f6cce;hpb=11fafa5d1eefeebd32e0996d2544adc39b16d8d1;p=physik%2Fposic.git diff --git a/mdrun.c b/mdrun.c index 1368880..34807e7 100644 --- a/mdrun.c +++ b/mdrun.c @@ -17,6 +17,11 @@ #include "potentials/tersoff.h" #endif +/* pse */ +#define PSE_MASS +#include "pse.h" +#undef PSE_MASS + #define ME "[mdrun]" /* @@ -111,6 +116,18 @@ int add_stage(t_mdrun *mdrun,u8 type,void *params) { case STAGE_CHSATTR: psize=sizeof(t_chsattr_params); break; + case STAGE_SET_TEMP: + psize=sizeof(t_set_temp_params); + break; + case STAGE_SET_TIMESTEP: + psize=sizeof(t_set_timestep_params); + break; + case STAGE_FILL: + psize=sizeof(t_fill_params); + break; + case STAGE_THERMAL_INIT: + psize=0; + break; default: printf("%s unknown stage type: %02x\n",ME,type); return -1; @@ -144,7 +161,7 @@ int mdrun_parse_config(t_mdrun *mdrun) { char error[128]; char line[128]; char *wptr; - char word[16][64]; + char word[32][64]; int wcnt; int i,o; @@ -155,6 +172,9 @@ int mdrun_parse_config(t_mdrun *mdrun) { t_anneal_params ap; t_chaattr_params cap; t_chsattr_params csp; + t_set_temp_params stp; + t_set_timestep_params stsp; + t_fill_params fp; /* open config file */ fd=open(mdrun->cfile,O_RDONLY); @@ -184,6 +204,9 @@ int mdrun_parse_config(t_mdrun *mdrun) { memset(&ap,0,sizeof(t_anneal_params)); memset(&cap,0,sizeof(t_chaattr_params)); memset(&csp,0,sizeof(t_chsattr_params)); + memset(&stp,0,sizeof(t_set_temp_params)); + memset(&stsp,0,sizeof(t_set_timestep_params)); + memset(&fp,0,sizeof(t_fill_params)); // get command + args wcnt=0; @@ -254,21 +277,119 @@ int mdrun_parse_config(t_mdrun *mdrun) { mdrun->lattice=DIAMOND; if(!strncmp(word[1],"none",4)) mdrun->lattice=NONE; + if(wcnt==3) + mdrun->lc=atof(word[2]); } else if(!strncmp(word[0],"element1",8)) { mdrun->element1=atoi(word[1]); - mdrun->m1=pse_mass[mdrun->element1]; } else if(!strncmp(word[0],"element2",8)) { mdrun->element2=atoi(word[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]); + // default values + fp.fill_element=mdrun->element1; + fp.fill_brand=0; + fp.lattice=mdrun->lattice; + fp.p_params.type=0; + fp.d_params.type=0; + fp.d_params.stype=0; + // parse fill command + i=1; + while(ilc=fp.lc; + } + if(!strncmp(word[i],"eb",2)) { + fp.fill_element=atoi(word[++i]); + fp.fill_brand=atoi(word[++i]); + } + if(word[i][0]=='p') { + i+=1; + switch(word[i][0]) { + case 'i': + if(word[i][1]=='r') + fp.p_params.type=PART_INSIDE_R; + else + fp.p_params.type=PART_INSIDE_D; + break; + case 'o': + if(word[i][1]=='r') + fp.p_params.type=PART_OUTSIDE_R; + else + fp.p_params.type=PART_OUTSIDE_D; + break; + default: + break; + } + if((fp.p_params.type==PART_INSIDE_R)|| + (fp.p_params.type==PART_OUTSIDE_R)) { + fp.p_params.r=atof(word[++i]); + fp.p_params.p.x=atof(word[++i]); + fp.p_params.p.y=atof(word[++i]); + fp.p_params.p.z=atof(word[++i]); + } + if((fp.p_params.type==PART_INSIDE_D)|| + (fp.p_params.type==PART_OUTSIDE_D)) { + fp.p_params.p.x=atof(word[++i]); + fp.p_params.p.y=atof(word[++i]); + fp.p_params.p.z=atof(word[++i]); + fp.p_params.d.x=atof(word[++i]); + fp.p_params.d.y=atof(word[++i]); + fp.p_params.d.z=atof(word[++i]); + } + } + if(word[i][0]=='d') { + switch(word[++i][0]) { + case '0': + + fp.d_params.type=DEFECT_TYPE_0D; + if(!strncmp(word[i+1],"dbx",3)) { + fp.d_params.stype=DEFECT_STYPE_DB_X; + } + if(!strncmp(word[i+1],"dby",3)) { + fp.d_params.stype=DEFECT_STYPE_DB_Y; + } + if(!strncmp(word[i+1],"dbz",3)) { + fp.d_params.stype=DEFECT_STYPE_DB_Z; + } + if(!strncmp(word[i+1],"dbr",3)) { + fp.d_params.stype=DEFECT_STYPE_DB_R; + } + i+=1; + fp.d_params.od=atof(word[++i]); + fp.d_params.dd=atof(word[++i]); + fp.d_params.element=atoi(word[++i]); + fp.d_params.brand=atoi(word[++i]); + // parsed in future + fp.d_params.attr=ATOM_ATTR_HB|ATOM_ATTR_VA; + fp.d_params.attr|=ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP; + break; + + case '1': + fp.d_params.type=DEFECT_TYPE_1D; + break; + case '2': + fp.d_params.type=DEFECT_TYPE_2D; + break; + case '3': + fp.d_params.type=DEFECT_TYPE_3D; + break; + default: + break; + } + + } + i+=1; + } + add_stage(mdrun,STAGE_FILL,&fp); + } + else if(!strncmp(word[0],"thermal_init",12)) { + add_stage(mdrun,STAGE_THERMAL_INIT,NULL); } else if(!strncmp(word[0],"aattr",5)) { // for aatrib line we need a special stage @@ -429,19 +550,32 @@ int mdrun_parse_config(t_mdrun *mdrun) { iap.z0=atof(word[10]); break; case 'r': - if(word[8][0]=='t') { + switch(word[8][0]) { + + case 't': iap.type=INS_TOTAL; iap.cr=atof(word[9]); - } - else { - iap.type=INS_REGION; - iap.x0=atof(word[8]); - iap.y0=atof(word[9]); - iap.z0=atof(word[10]); - iap.x1=atof(word[11]); - iap.y1=atof(word[12]); - iap.z1=atof(word[13]); - iap.cr=atof(word[14]); + break; + case 'r': + iap.type=INS_RECT; + iap.x0=atof(word[9]); + iap.y0=atof(word[10]); + iap.z0=atof(word[11]); + iap.x1=atof(word[12]); + iap.y1=atof(word[13]); + iap.z1=atof(word[14]); + iap.cr=atof(word[15]); + break; + case 's': + iap.type=INS_SPHERE; + iap.x0=atof(word[9]); + iap.y0=atof(word[10]); + iap.z0=atof(word[11]); + iap.x1=atof(word[12]); + iap.cr=atof(word[13]); + break; + default: + break; } default: break; @@ -507,8 +641,24 @@ int mdrun_parse_config(t_mdrun *mdrun) { ap.count=0; ap.runs=atoi(word[2]); ap.dt=atof(word[3]); + ap.interval=atoi(word[4]); add_stage(mdrun,STAGE_ANNEAL,&ap); } + else if(!strncmp(word[1],"set_temp",8)) { + if(word[2][0]=='c') { + stp.type=SET_TEMP_CURRENT; + stp.val=0.0; + } + else { + stp.type=SET_TEMP_VALUE; + stp.val=atof(word[2]); + } + add_stage(mdrun,STAGE_SET_TEMP,&stp); + } + else if(!strncmp(word[1],"set_timestep",12)) { + stsp.tau=atof(word[2]); + add_stage(mdrun,STAGE_SET_TIMESTEP,&stsp); + } else { printf("%s unknown stage type: %s\n", ME,word[1]); @@ -628,7 +778,7 @@ int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { z0=-z/2.0; cr_check=TRUE; break; - case INS_REGION: + case INS_RECT: x=iap->x1-iap->x0; x0=iap->x0; y=iap->y1-iap->y0; @@ -637,6 +787,15 @@ int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { z0=iap->z0; cr_check=TRUE; break; + case INS_SPHERE: + x=2.0*iap->x1; + x0=iap->x0-iap->x1; + y=x; + y0=iap->y0-iap->x1; + z=x; + z0=iap->z0-iap->x1; + cr_check=TRUE; + break; case INS_POS: x0=iap->x0; y0=iap->y0; @@ -688,8 +847,16 @@ int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { dmin=d; } } + if(iap->type==INS_SPHERE) { + if((r.x-iap->x0)*(r.x-iap->x0)+ + (r.y-iap->y0)*(r.y-iap->y0)+ + (r.z-iap->z0)*(r.z-iap->z0)> + (iap->x1*iap->x1)) { + run=1; + } + } } - add_atom(moldyn,iap->element,pse_mass[iap->element], + add_atom(moldyn,iap->element, iap->brand,iap->attr,&r,&v); printf("%s atom inserted (%d/%d): %f %f %f\n", ME,(iap->cnt_steps+1)*iap->ins_atoms, @@ -758,7 +925,7 @@ int insert_mixed_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { if(dmin>cmax) retry=1; } - add_atom(moldyn,imp->element1,pse_mass[imp->element1], + add_atom(moldyn,imp->element1, imp->brand1,imp->attr1,&r,&v); printf("%s (mixed) atom inserted (%d): %f %f %f\n", ME,imp->amount1,r.x,r.y,r.z); @@ -792,7 +959,7 @@ int insert_mixed_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { if(dmin>cmax) retry=1; } - add_atom(moldyn,imp->element2,pse_mass[imp->element2], + add_atom(moldyn,imp->element2, imp->brand2,imp->attr2,&r,&v); printf("%s (mixed) atom inserted (%d): %f %f %f\n", ME,imp->amount2,r.x,r.y,r.z); @@ -897,13 +1064,16 @@ int mdrun_hook(void *ptr1,void *ptr2) { t_stage *stage; t_list *sl; int steps; - double tau; u8 change_stage; + t_3dvec o; t_insert_atoms_params *iap; t_insert_mixed_atoms_params *imp; t_continue_params *cp; t_anneal_params *ap; + t_set_temp_params *stp; + t_set_timestep_params *stsp; + t_fill_params *fp; moldyn=ptr1; mdrun=ptr2; @@ -919,9 +1089,8 @@ int mdrun_hook(void *ptr1,void *ptr2) { /* get stage description */ stage=sl->current->data; - /* default steps and tau values */ + /* steps in next schedule */ steps=mdrun->relax_steps; - tau=mdrun->timestep; /* check whether relaxation steps are necessary */ if((check_pressure(moldyn,mdrun)==TRUE)&\ @@ -981,17 +1150,85 @@ int mdrun_hook(void *ptr1,void *ptr2) { set_temperature(moldyn, moldyn->t_ref+ap->dt); ap->count+=1; + steps=ap->interval; break; case STAGE_CHAATTR: - stage_print(" -> chaattr\n\n"); + stage_print(" -> change atom attributes\n\n"); chaatr(moldyn,mdrun); change_stage=TRUE; break; case STAGE_CHSATTR: - stage_print(" -> chsattr\n\n"); + stage_print(" -> change sys attributes\n\n"); chsattr(moldyn,mdrun); change_stage=TRUE; break; + case STAGE_SET_TEMP: + stage_print(" -> set temperature\n\n"); + stp=stage->params; + if(stp->type&SET_TEMP_CURRENT) { + set_temperature(moldyn,moldyn->t_avg); + } + else { + set_temperature(moldyn,stp->val); + } + change_stage=TRUE; + break; + case STAGE_SET_TIMESTEP: + stage_print(" -> set timestep\n\n"); + stsp=stage->params; + mdrun->timestep=stsp->tau; + change_stage=TRUE; + break; + case STAGE_FILL: + stage_print(" -> fill lattice\n\n"); + fp=stage->params; + switch(fp->lattice) { + case ZINCBLENDE: + + o.x=0.5*0.25*fp->lc; + o.y=o.x; + o.z=o.x; + create_lattice(moldyn, + FCC,fp->lc, + mdrun->element1, + DEFAULT_ATOM_ATTR,0, + fp->lx,fp->ly,fp->lz, + &o, + &(fp->p_params), + &(fp->d_params)); + o.x+=0.25*fp->lc; + o.y=o.x; + o.z=o.x; + create_lattice(moldyn, + FCC,fp->lc, + mdrun->element2, + DEFAULT_ATOM_ATTR,1, + fp->lx,fp->ly,fp->lz, + &o, + &(fp->p_params), + &(fp->d_params)); + break; + + default: + + create_lattice(moldyn, + fp->lattice,fp->lc, + mdrun->element1, + DEFAULT_ATOM_ATTR,0, + fp->lx,fp->ly,fp->lz, + NULL, + &(fp->p_params), + &(fp->d_params)); + break; + } + moldyn_bc_check(moldyn); + change_stage=TRUE; + break; + case STAGE_THERMAL_INIT: + stage_print(" -> thermal init\n\n"); + thermal_init(moldyn,TRUE); + change_stage=TRUE; + break; default: printf("%s unknwon stage type\n",ME); break; @@ -1008,7 +1245,6 @@ int mdrun_hook(void *ptr1,void *ptr2) { return 0; } steps=0; - tau=mdrun->timestep; } } @@ -1021,7 +1257,7 @@ int mdrun_hook(void *ptr1,void *ptr2) { } /* continue simulation */ - moldyn_add_schedule(moldyn,steps,tau); + moldyn_add_schedule(moldyn,steps,mdrun->timestep); return 0; } @@ -1030,7 +1266,7 @@ int main(int argc,char **argv) { t_mdrun mdrun; t_moldyn moldyn; - t_3dvec o; + //t_3dvec o; /* clear structs */ memset(&mdrun,0,sizeof(t_mdrun)); @@ -1105,28 +1341,33 @@ int main(int argc,char **argv) { /* 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); + /* replaced by fill stage !! switch(mdrun.lattice) { case FCC: - create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element1, - mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx, - mdrun.ly,mdrun.lz,NULL); + create_lattice(&moldyn,FCC,mdrun.lc,mdrun.fill_element, + mdrun.m1,DEFAULT_ATOM_ATTR, + mdrun.fill_brand,mdrun.lx, + mdrun.ly,mdrun.lz,NULL,0,NULL); break; case DIAMOND: - create_lattice(&moldyn,DIAMOND,mdrun.lc,mdrun.element1, - mdrun.m1,DEFAULT_ATOM_ATTR,0,mdrun.lx, - mdrun.ly,mdrun.lz,NULL); + create_lattice(&moldyn,DIAMOND,mdrun.lc, + mdrun.fill_element, + mdrun.m1,DEFAULT_ATOM_ATTR, + mdrun.fill_brand,mdrun.lx, + mdrun.ly,mdrun.lz,NULL,0,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,DEFAULT_ATOM_ATTR,0,mdrun.lx, - mdrun.ly,mdrun.lz,&o); + mdrun.ly,mdrun.lz,&o,0,NULL); o.x+=0.25*mdrun.lc; o.y=o.x; o.z=o.x; create_lattice(&moldyn,FCC,mdrun.lc,mdrun.element2, mdrun.m2,DEFAULT_ATOM_ATTR,1,mdrun.lx, - mdrun.ly,mdrun.lz,&o); + mdrun.ly,mdrun.lz,&o,0,NULL); break; case NONE: + set_nn_dist(&moldyn,mdrun.nnd); break; default: printf("%s unknown lattice type: %02x\n", @@ -1134,11 +1375,14 @@ int main(int argc,char **argv) { return -1; } moldyn_bc_check(&moldyn); + replaced by fill stage !! */ /* temperature and pressure */ set_temperature(&moldyn,mdrun.temperature); set_pressure(&moldyn,mdrun.pressure); + /* replaced thermal_init stage thermal_init(&moldyn,TRUE); + */ addsched: /* first schedule */