From: hackbard Date: Wed, 17 Dec 2008 16:18:12 +0000 (+0100) Subject: cleaned partial lattice create + added defect feature, adabted md code + X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=dfb47bde6bf974d1b231048ec8bcda712404ccc1;p=physik%2Fposic.git cleaned partial lattice create + added defect feature, adabted md code + some minor changes: ask diff! :) - xmas checkin ... --- diff --git a/TODO.txt b/TODO.txt index 570ee54..8712c9d 100644 --- a/TODO.txt +++ b/TODO.txt @@ -1,19 +1,23 @@ interstitial simulations ************************ - DOUBLECHECK with PARCAS !!!! - - - visualize md.movie (converter tool - maybe good!) - - more checking ... +- more DOUBLECHECKs with PARCAS !!!! implement ********* - improve ins_m_atoms (merge to general ins_atoms function maybe) + + * modified first: try rules to better distribute Si and C *** + * general cleanup + * type 1 square, type 2 radius! + - clean up mdrun - config sanity checks - introduce/improve fill command (multiple fills) + - anneal -> from current to T with rate R + - check virial calc, where does the - come from? - angular distribution @@ -25,8 +29,18 @@ implement - improve diff calc - - make it parallel! (mpi/openmp) <- email to ralfu (asap) - * openmp: doch auf verlet listen (pointer problem!) + - pthreads: + + * threads nur einmal oeffnen + * verteilung auf laufende threads + + - potentials: + + * tersoff fast *** + * modifief tersoff/albe *** + * stillinger weber *** + * eam + * edip - optimize code! @@ -44,6 +58,9 @@ implement * -> 3.50 ^ + static lists * -> 3.37 source c, albe f, arch opts, c,d,h,gamma + lowmem lists (test for bigger + atoms!) + * -> 3.31(36) ^ + inline v_calc + + * -> 4.44 orig albe, lowmem, source c, arch opts todo: - listen ! estimate time @@ -61,7 +78,7 @@ simulation runs - tctrl only in outer regions - only 1 atom per timestep - EXTENDED C-C cutoff -- different sized SiC prec in Si (4:5! diff temperatures) +- different sized SiC prec in Si (4:5! diff temperatures) *** - melting exps (both, anneal + interface method) - interstitials: - more interstitial combinations diff --git a/calc_delta_e b/calc_delta_e index 4ff46fb..455ce57 100755 --- a/calc_delta_e +++ b/calc_delta_e @@ -2,7 +2,7 @@ file=`ls $1/atomic_conf_* | tail -1` atom_cnt=`grep '# \[P\]' $file | awk '{ print $3 }'` -e0=`grep ^0 $1/energy | awk '{ print $4 }'` +e0=`awk 'NR==2' $1/energy | awk '{ print $4 }'` e1=`tail -n 1 $1/energy | awk '{ print $4 }'` echo "$e0 $e1 $atom_cnt" | \ diff --git a/mdrun.c b/mdrun.c index 6effe1a..d063693 100644 --- a/mdrun.c +++ b/mdrun.c @@ -277,66 +277,114 @@ 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)) { - fp.lx=atoi(word[2]); - fp.ly=atoi(word[3]); - fp.lz=atoi(word[4]); - fp.lc=atof(word[5]); - mdrun->lc=fp.lc; - if(!strncmp(word[1],"lc",2)) { - if(wcnt==8) { - fp.fill_element=atoi(word[6]); - fp.fill_brand=atoi(word[7]); + // 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; + // parse fill command + i=1; + while(ilc=fp.lc; } - else { - fp.fill_element=mdrun->element1; - fp.fill_brand=0; + if(!strncmp(word[i],"eb",2)) { + fp.fill_element=atoi(word[++i]); + fp.fill_brand=atoi(word[++i]); } - } - else { - switch(word[6][0]) { + if(word[i][0]=='p') { + i+=1; + switch(word[i][0]) { case 'i': - if(word[6][1]=='r') - fp.p_type=PART_INSIDE_R; + if(word[i][1]=='r') + fp.p_params.type=PART_INSIDE_R; else - fp.p_type=PART_INSIDE_D; + fp.p_params.type=PART_INSIDE_D; break; case 'o': - if(word[6][1]=='r') - fp.p_type=PART_OUTSIDE_R; + if(word[i][1]=='r') + fp.p_params.type=PART_OUTSIDE_R; else - fp.p_type=PART_OUTSIDE_D; + 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; } - if((fp.p_type==PART_INSIDE_R)|| - (fp.p_type==PART_OUTSIDE_R)) { - fp.p_vals.r=atof(word[7]); - fp.p_vals.p.x=atof(word[8]); - fp.p_vals.p.y=atof(word[9]); - fp.p_vals.p.z=atof(word[10]); - } - if((fp.p_type==PART_INSIDE_D)|| - (fp.p_type==PART_OUTSIDE_D)) { - fp.p_vals.p.x=atof(word[7]); - fp.p_vals.p.y=atof(word[8]); - fp.p_vals.p.z=atof(word[9]); - fp.p_vals.d.x=atof(word[10]); - fp.p_vals.d.y=atof(word[11]); - fp.p_vals.d.z=atof(word[12]); - } - fp.lattice=mdrun->lattice; add_stage(mdrun,STAGE_FILL,&fp); } else if(!strncmp(word[0],"thermal_init",12)) { @@ -777,7 +825,7 @@ int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { } } } - 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, @@ -846,7 +894,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); @@ -880,7 +928,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); @@ -1103,39 +1151,44 @@ int mdrun_hook(void *ptr1,void *ptr2) { case STAGE_FILL: stage_print(" -> fill lattice\n\n"); fp=stage->params; - if(fp->lattice!=ZINCBLENDE) { - create_lattice(moldyn, - fp->lattice,fp->lc, - mdrun->element1, - mdrun->m1, - DEFAULT_ATOM_ATTR,0, - fp->lx,fp->ly,fp->lz, - NULL,fp->p_type, - &(fp->p_vals)); - } - else { + 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, - mdrun->m1, DEFAULT_ATOM_ATTR,0, fp->lx,fp->ly,fp->lz, - &o,fp->p_type, - &(fp->p_vals)); + &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, - mdrun->m2, DEFAULT_ATOM_ATTR,1, fp->lx,fp->ly,fp->lz, - &o,fp->p_type, - &(fp->p_vals)); + &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; diff --git a/mdrun.h b/mdrun.h index d6476e5..fad1f62 100644 --- a/mdrun.h +++ b/mdrun.h @@ -71,9 +71,7 @@ typedef struct s_mdrun { u8 pbcz; int element1; // element 1 - double m1; int element2; // element 2 - double m2; double lc; // lattice constant u8 lattice; // type of lattice @@ -199,8 +197,8 @@ typedef struct s_fill_params { u8 lattice; int fill_element; u8 fill_brand; - u8 p_type; - t_part_vals p_vals; + t_part_params p_params; + t_defect_params d_params; } t_fill_params; /* diff --git a/moldyn.c b/moldyn.c index a934194..bc4fdff 100644 --- a/moldyn.c +++ b/moldyn.c @@ -41,9 +41,11 @@ #endif /* pse */ +#define PSE_MASS #define PSE_NAME #define PSE_COL #include "pse.h" +#undef PSE_MASS #undef PSE_NAME #undef PSE_COL @@ -512,9 +514,9 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { * creating lattice functions */ -int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, +int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element, u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin, - u8 p_type,t_part_vals *p_vals) { + t_part_params *p_params,t_defect_params *d_params) { int new,count; int ret; @@ -538,6 +540,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, if(type==FCC) new*=4; if(type==DIAMOND) new*=8; + /* defects */ + if(d_params->type) { + switch(d_params->stype) { + case DEFECT_STYPE_DB_X: + case DEFECT_STYPE_DB_Y: + case DEFECT_STYPE_DB_Z: + case DEFECT_STYPE_DB_R: + new*=2; + break; + default: + printf("[moldyn] WARNING: cl unknown defect\n"); + break; + } + } + /* allocate space for atoms */ ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom)); if(!ptr) { @@ -572,21 +589,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, switch(type) { case CUBIC: set_nn_dist(moldyn,lc); - ret=cubic_init(a,b,c,lc,atom,&orig,p_type,p_vals); + ret=cubic_init(a,b,c,lc,atom,&orig,p_params,d_params); strcpy(name,"cubic"); break; case FCC: if(!origin) v3_scale(&orig,&orig,0.5); set_nn_dist(moldyn,0.5*sqrt(2.0)*lc); - ret=fcc_init(a,b,c,lc,atom,&orig,p_type,p_vals); + ret=fcc_init(a,b,c,lc,atom,&orig,p_params,d_params); strcpy(name,"fcc"); break; case DIAMOND: if(!origin) v3_scale(&orig,&orig,0.25); set_nn_dist(moldyn,0.25*sqrt(3.0)*lc); - ret=diamond_init(a,b,c,lc,atom,&orig,p_type,p_vals); + ret=diamond_init(a,b,c,lc,atom,&orig,p_params,d_params); strcpy(name,"diamond"); break; default: @@ -598,7 +615,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, if(ret!=new) { printf("[moldyn] creating %s lattice (lc=%f) incomplete\n", name,lc); - printf(" (ignore in case of partial lattice creation)\n"); + printf(" (ignore for partial lattice creation)\n"); printf(" amount of atoms\n"); printf(" - expected: %d\n",new); printf(" - created: %d\n",ret); @@ -610,7 +627,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, for(new=0;newtype) { + new+=1; + atom[new].element=d_params->element; + atom[new].mass=pse_mass[d_params->element]; + atom[new].attr=d_params->attr; + atom[new].brand=d_params->brand; + atom[new].tag=count+new; + check_per_bound(moldyn,&(atom[new].r)); + atom[new].r_0=atom[new].r; +#ifdef PTHREADS + pthread_mutex_init(&(mutex[new]),NULL); +#endif + } } /* fix allocation */ @@ -629,6 +659,8 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, } moldyn->atom=ptr; +// WHAT ABOUT AMUTEX !!!! + #ifdef LOWMEM_LISTS ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int)); if(!ptr) { @@ -644,7 +676,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, return ret; } -int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, +int add_atom(t_moldyn *moldyn,int element,u8 brand,u8 attr, t_3dvec *r,t_3dvec *v) { t_atom *atom; @@ -687,7 +719,7 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, atom[count].r=*r; atom[count].v=*v; atom[count].element=element; - atom[count].mass=mass; + atom[count].mass=pse_mass[element]; atom[count].brand=brand; atom[count].tag=count; atom[count].attr=attr; @@ -729,9 +761,44 @@ int del_atom(t_moldyn *moldyn,int tag) { return 0; } +#define set_atom_positions(pos) \ + if(d_params->type) {\ + d_o.x=0; d_o.y=0; d_o.z=0;\ + d_d.x=0; d_d.y=0; d_d.z=0;\ + switch(d_params->stype) {\ + case DEFECT_STYPE_DB_X:\ + d_o.x=d_params->od;\ + d_d.x=d_params->dd;\ + break;\ + case DEFECT_STYPE_DB_Y:\ + d_o.y=d_params->od;\ + d_d.y=d_params->dd;\ + break;\ + case DEFECT_STYPE_DB_Z:\ + d_o.z=d_params->od;\ + d_d.z=d_params->dd;\ + break;\ + case DEFECT_STYPE_DB_R:\ + break;\ + default:\ + printf("[moldyn] WARNING: unknown defect\n");\ + break;\ + }\ + v3_add(&dr,&pos,&d_o);\ + v3_copy(&(atom[count].r),&dr);\ + count+=1;\ + v3_add(&dr,&pos,&d_d);\ + v3_copy(&(atom[count].r),&dr);\ + count+=1;\ + }\ + else {\ + v3_copy(&(atom[count].r),&pos);\ + count+=1;\ + } + /* cubic init */ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, - u8 p_type,t_part_vals *p_vals) { + t_part_params *p_params,t_defect_params *d_params) { int count; t_3dvec r; @@ -739,6 +806,9 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, t_3dvec o; t_3dvec dist; t_3dvec p; + t_3dvec d_o; + t_3dvec d_d; + t_3dvec dr; p.x=0; p.y=0; p.z=0; @@ -749,10 +819,10 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, v3_zero(&o); /* shift partition values */ - if(p_type) { - p.x=p_vals->p.x+(a*lc)/2.0; - p.y=p_vals->p.y+(b*lc)/2.0; - p.z=p_vals->p.z+(c*lc)/2.0; + if(p_params->type) { + p.x=p_params->p.x+(a*lc)/2.0; + p.y=p_params->p.y+(b*lc)/2.0; + p.z=p_params->p.z+(c*lc)/2.0; } r.x=o.x; @@ -761,42 +831,39 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, for(j=0;jtype) { case PART_INSIDE_R: v3_sub(&dist,&r,&p); - if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&r); - count+=1; + if(v3_absolute_square(&dist)< + (p_params->r*p_params->r)) { + set_atom_positions(r); } break; case PART_OUTSIDE_R: v3_sub(&dist,&r,&p); - if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&r); - count+=1; + if(v3_absolute_square(&dist)>= + (p_params->r*p_params->r)) { + set_atom_positions(r); } break; case PART_INSIDE_D: v3_sub(&dist,&r,&p); - if((fabs(dist.x)d.x)&& - (fabs(dist.y)d.y)&& - (fabs(dist.z)d.z)) { - v3_copy(&(atom[count].r),&r); - count+=1; + if((fabs(dist.x)d.x)&& + (fabs(dist.y)d.y)&& + (fabs(dist.z)d.z)) { + set_atom_positions(r); } break; case PART_OUTSIDE_D: v3_sub(&dist,&r,&p); - if((fabs(dist.x)>=p_vals->d.x)|| - (fabs(dist.y)>=p_vals->d.y)|| - (fabs(dist.z)>=p_vals->d.z)) { - v3_copy(&(atom[count].r),&r); - count+=1; + if((fabs(dist.x)>=p_params->d.x)|| + (fabs(dist.y)>=p_params->d.y)|| + (fabs(dist.z)>=p_params->d.z)) { + set_atom_positions(r); } break; default: - v3_copy(&(atom[count].r),&r); - count+=1; + set_atom_positions(r); break; } r.z+=lc; @@ -817,7 +884,7 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, /* fcc lattice init */ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, - u8 p_type,t_part_vals *p_vals) { + t_part_params *p_params,t_defect_params *d_params) { int count; int i,j,k,l; @@ -825,6 +892,7 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, t_3dvec basis[3]; t_3dvec dist; t_3dvec p; + t_3dvec d_d,d_o,dr; p.x=0; p.y=0; p.z=0; @@ -844,10 +912,10 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, basis[2].z=0.5*lc; /* shift partition values */ - if(p_type) { - p.x=p_vals->p.x+(a*lc)/2.0; - p.y=p_vals->p.y+(b*lc)/2.0; - p.z=p_vals->p.z+(c*lc)/2.0; + if(p_params->type) { + p.x=p_params->p.x+(a*lc)/2.0; + p.y=p_params->p.y+(b*lc)/2.0; + p.z=p_params->p.z+(c*lc)/2.0; } /* fill up the room */ @@ -858,83 +926,77 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, r.z=o.z; for(k=0;ktype) { case PART_INSIDE_R: v3_sub(&dist,&r,&p); - if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&r); - count+=1; + if(v3_absolute_square(&dist)< + (p_params->r*p_params->r)) { + set_atom_positions(r); } break; case PART_OUTSIDE_R: v3_sub(&dist,&r,&p); - if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&r); - count+=1; + if(v3_absolute_square(&dist)>= + (p_params->r*p_params->r)) { + set_atom_positions(r); } break; case PART_INSIDE_D: v3_sub(&dist,&r,&p); - if((fabs(dist.x)d.x)&& - (fabs(dist.y)d.y)&& - (fabs(dist.z)d.z)) { - v3_copy(&(atom[count].r),&r); - count+=1; + if((fabs(dist.x)d.x)&& + (fabs(dist.y)d.y)&& + (fabs(dist.z)d.z)) { + set_atom_positions(r); } break; case PART_OUTSIDE_D: v3_sub(&dist,&r,&p); - if((fabs(dist.x)>=p_vals->d.x)|| - (fabs(dist.y)>=p_vals->d.y)|| - (fabs(dist.z)>=p_vals->d.z)) { - v3_copy(&(atom[count].r),&r); - count+=1; + if((fabs(dist.x)>=p_params->d.x)|| + (fabs(dist.y)>=p_params->d.y)|| + (fabs(dist.z)>=p_params->d.z)) { + set_atom_positions(r); } break; default: - v3_copy(&(atom[count].r),&r); - count+=1; + set_atom_positions(r); break; } /* the three face centered atoms */ for(l=0;l<3;l++) { v3_add(&n,&r,&basis[l]); - switch(p_type) { + switch(p_params->type) { case PART_INSIDE_R: v3_sub(&dist,&n,&p); - if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&n); - count+=1; + if(v3_absolute_square(&dist)< + (p_params->r*p_params->r)) { + set_atom_positions(n); } break; case PART_OUTSIDE_R: v3_sub(&dist,&n,&p); - if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) { - v3_copy(&(atom[count].r),&n); - count+=1; + if(v3_absolute_square(&dist)>= + (p_params->r*p_params->r)) { + set_atom_positions(n); } break; case PART_INSIDE_D: v3_sub(&dist,&n,&p); - if((fabs(dist.x)d.x)&& - (fabs(dist.y)d.y)&& - (fabs(dist.z)d.z)) { - v3_copy(&(atom[count].r),&n); - count+=1; + if((fabs(dist.x)d.x)&& + (fabs(dist.y)d.y)&& + (fabs(dist.z)d.z)) { + set_atom_positions(n); } break; case PART_OUTSIDE_D: v3_sub(&dist,&n,&p); - if((fabs(dist.x)>=p_vals->d.x)|| - (fabs(dist.y)>=p_vals->d.y)|| - (fabs(dist.z)>=p_vals->d.z)) { - v3_copy(&(atom[count].r),&n); - count+=1; + if((fabs(dist.x)>=p_params->d.x)|| + (fabs(dist.y)>=p_params->d.y)|| + (fabs(dist.z)>=p_params->d.z)) { + set_atom_positions(n); } break; default: - v3_copy(&(atom[count].r),&n); - count+=1; + set_atom_positions(n); break; } } @@ -956,12 +1018,12 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, } int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, - u8 p_type,t_part_vals *p_vals) { + t_part_params *p_params,t_defect_params *d_params) { int count; t_3dvec o; - count=fcc_init(a,b,c,lc,atom,origin,p_type,p_vals); + count=fcc_init(a,b,c,lc,atom,origin,p_params,d_params); o.x=0.25*lc; o.y=0.25*lc; @@ -969,7 +1031,7 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, if(origin) v3_add(&o,&o,origin); - count+=fcc_init(a,b,c,lc,&atom[count],&o,p_type,p_vals); + count+=fcc_init(a,b,c,lc,&atom[count],&o,p_params,d_params); return count; } @@ -1171,12 +1233,12 @@ double pressure_calc(t_moldyn *moldyn) { moldyn->p=2.0*moldyn->ekin+moldyn->virial; moldyn->p/=(3.0*moldyn->volume); - moldyn->px=2.0*moldyn->ekinx+moldyn->vir.xx; - moldyn->px/=moldyn->volume; - moldyn->py=2.0*moldyn->ekiny+moldyn->vir.yy; - moldyn->py/=moldyn->volume; - moldyn->pz=2.0*moldyn->ekinz+moldyn->vir.zz; - moldyn->pz/=moldyn->volume; + //moldyn->px=2.0*moldyn->ekinx+moldyn->vir.xx; + //moldyn->px/=moldyn->volume; + //moldyn->py=2.0*moldyn->ekiny+moldyn->vir.yy; + //moldyn->py/=moldyn->volume; + //moldyn->pz=2.0*moldyn->ekinz+moldyn->vir.zz; + //moldyn->pz/=moldyn->volume; /* pressure (absolute coordinates) */ //moldyn->gp=2.0*moldyn->ekin+moldyn->gv; @@ -1439,16 +1501,15 @@ int scale_dim_ind(t_moldyn *moldyn,double x,double y,double z) { int scale_volume(t_moldyn *moldyn) { t_3dvec *dim,*vdim; - //double scale; + double scale; t_linkcell *lc; - double sx,sy,sz; + //double sx,sy,sz; vdim=&(moldyn->vis.dim); dim=&(moldyn->dim); lc=&(moldyn->lc); /* scaling factor */ - /* if(moldyn->pt_scale&P_SCALE_BERENDSEN) { scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau; scale=pow(scale,ONE_THIRD); @@ -1456,20 +1517,22 @@ int scale_volume(t_moldyn *moldyn) { else { scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD); } - */ + + /* sx=1.0-(moldyn->p_ref-moldyn->px)*moldyn->p_tc*moldyn->tau; sy=1.0-(moldyn->p_ref-moldyn->py)*moldyn->p_tc*moldyn->tau; sz=1.0-(moldyn->p_ref-moldyn->pz)*moldyn->p_tc*moldyn->tau; sx=pow(sx,ONE_THIRD); sy=pow(sy,ONE_THIRD); sz=pow(sz,ONE_THIRD); + */ /* scale the atoms and dimensions */ - //scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); - //scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); - scale_atoms_ind(moldyn,sx,sy,sz); - scale_dim_ind(moldyn,sx,sy,sz); + scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE); + //scale_atoms_ind(moldyn,sx,sy,sz); + //scale_dim_ind(moldyn,sx,sy,sz); /* visualize dimensions */ if(vdim->x!=0) { @@ -1488,12 +1551,12 @@ int scale_volume(t_moldyn *moldyn) { link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); } else { - //lc->x*=scale; - //lc->y*=scale; - //lc->z*=scale; - lc->x*=sx; - lc->y*=sx; - lc->z*=sy; + lc->x*=scale; + lc->y*=scale; + lc->z*=scale; + //lc->x*=sx; + //lc->y*=sx; + //lc->z*=sy; } return 0; @@ -1507,16 +1570,16 @@ double e_kin_calc(t_moldyn *moldyn) { atom=moldyn->atom; moldyn->ekin=0.0; - moldyn->ekinx=0.0; - moldyn->ekiny=0.0; - moldyn->ekinz=0.0; + //moldyn->ekinx=0.0; + //moldyn->ekiny=0.0; + //moldyn->ekinz=0.0; for(i=0;icount;i++) { atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); moldyn->ekin+=atom[i].ekin; - moldyn->ekinx+=0.5*atom[i].mass*atom[i].v.x*atom[i].v.x; - moldyn->ekiny+=0.5*atom[i].mass*atom[i].v.y*atom[i].v.y; - moldyn->ekinz+=0.5*atom[i].mass*atom[i].v.z*atom[i].v.z; + //moldyn->ekinx+=0.5*atom[i].mass*atom[i].v.x*atom[i].v.x; + //moldyn->ekiny+=0.5*atom[i].mass*atom[i].v.y*atom[i].v.y; + //moldyn->ekinz+=0.5*atom[i].mass*atom[i].v.z*atom[i].v.z; } return moldyn->ekin; @@ -3264,6 +3327,7 @@ int visual_atoms(t_moldyn *moldyn) { t_visual *v; t_atom *atom; t_vb vb; + t_3dvec strain; #ifdef VISUAL_THREAD t_moldyn *moldyn; @@ -3294,14 +3358,18 @@ int visual_atoms(t_moldyn *moldyn) { moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help); // atomic configuration - for(i=0;icount;i++) + for(i=0;icount;i++) { + v3_sub(&strain,&(atom[i].r),&(atom[i].r_0)); + check_per_bound(moldyn,&strain); // atom type, positions, color and kinetic energy dprintf(vb.fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element], atom[i].r.x, atom[i].r.y, atom[i].r.z, pse_col[atom[i].element], - atom[i].ekin); + //atom[i].ekin); + sqrt(v3_absolute_square(&strain))); + } // bonds between atoms #ifndef VISUAL_THREAD diff --git a/moldyn.h b/moldyn.h index d065d43..538e0eb 100644 --- a/moldyn.h +++ b/moldyn.h @@ -258,17 +258,38 @@ typedef struct s_vb { int fd; } t_vb; -typedef struct s_part_vals { +typedef struct s_part_params { + u8 type; double r; t_3dvec p; t_3dvec d; -} t_part_vals; +} t_part_params; #define PART_INSIDE_R 1 #define PART_OUTSIDE_R 2 #define PART_INSIDE_D 3 #define PART_OUTSIDE_D 4 +typedef struct s_defect_params { + u8 type; + u8 stype; + double od; + double dd; + int element; + u8 brand; + u8 attr; +} t_defect_params; + +#define DEFECT_TYPE_0D 1 +#define DEFECT_TYPE_1D 2 +#define DEFECT_TYPE_2D 3 +#define DEFECT_TYPE_3D 4 + +#define DEFECT_STYPE_DB_X 1 +#define DEFECT_STYPE_DB_Y 2 +#define DEFECT_STYPE_DB_Z 3 +#define DEFECT_STYPE_DB_R 4 + /* * * defines @@ -415,18 +436,18 @@ int moldyn_set_report(t_moldyn *moldyn,char *author,char *title); int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer); int moldyn_log_shutdown(t_moldyn *moldyn); -int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, +int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element, u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin, - u8 p_type,t_part_vals *p_vals); -int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, + t_part_params *p_params,t_defect_params *d_params); +int add_atom(t_moldyn *moldyn,int element,u8 brand,u8 attr, t_3dvec *r,t_3dvec *v); int del_atom(t_moldyn *moldyn,int tag); int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, - u8 p_type,t_part_vals *p_vals); + t_part_params *p_params,t_defect_params *d_params); int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, - u8 p_type,t_part_vals *p_vals); + t_part_params *p_params,t_defect_params *d_params); int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, - u8 p_type,t_part_vals *p_vals); + t_part_params *p_params,t_defect_params *d_params); int destroy_atoms(t_moldyn *moldyn); int thermal_init(t_moldyn *moldyn,u8 equi_init); diff --git a/remember_me.txt b/remember_me.txt index 7247fbb..3ae541e 100644 --- a/remember_me.txt +++ b/remember_me.txt @@ -11,3 +11,13 @@ c-c 3.1: saves/c_in_si_prec_450_01/s-0566000.save 239458 240660 110 - 100 239703 241605 100 - 100 + + +concatenated, differently oriented 100 dumbbells +################################################ + +posic_new/saves/c_in_si_prec_450_tot_02/s-0272000.save b 3.08 0.01 + + atoms 238338/1 238631/1 - 3.078045 + atoms 238461/1 244278/1 - 3.085715 + diff --git a/sic.c b/sic.c index f769f6f..375c6d5 100644 --- a/sic.c +++ b/sic.c @@ -121,7 +121,7 @@ int insert_atoms(t_moldyn *moldyn) { dmin=d; } } - add_atom(moldyn,INS_TYPE,INS_MASS,INS_BRAND, + add_atom(moldyn,INS_TYPE,INS_BRAND, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|\ INS_ATTR, &r,&v); @@ -325,23 +325,23 @@ int main(int argc,char **argv) { // diamond #ifdef ALBE #ifdef INIT_SI - create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,M_SI, + create_lattice(&md,DIAMOND,ALBE_LC_SI,SI, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, 0,LCNTX,LCNTY,LCNTZ,NULL,0,NULL); #endif #ifdef INIT_C - create_lattice(&md,DIAMOND,ALBE_LC_C,C,M_C, + create_lattice(&md,DIAMOND,ALBE_LC_C,C, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, 1,LCNTX,LCNTY,LCNTZ,NULL,0,NULL); #endif #else #ifdef INIT_SI - create_lattice(&md,DIAMOND,LC_SI,SI,M_SI, + create_lattice(&md,DIAMOND,LC_SI,SI, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, 0,LCNTX,LCNTY,LCNTZ,NULL,0,NULL); #endif #ifdef INIT_C - create_lattice(&md,DIAMOND,LC_C,SI,M_SI, + create_lattice(&md,DIAMOND,LC_C,SI, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, 1,LCNTX,LCNTY,LCNTZ,NULL,0,NULL); #endif @@ -351,20 +351,20 @@ int main(int argc,char **argv) { #ifdef INIT_3CSIC #ifdef ALBE r.x=0.5*0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x; - create_lattice(&md,FCC,ALBE_LC_SIC,SI,M_SI, + create_lattice(&md,FCC,ALBE_LC_SIC,SI, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, 0,LCNTX,LCNTY,LCNTZ,&r,0,NULL); r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x; - create_lattice(&md,FCC,ALBE_LC_SIC,C,M_C, + create_lattice(&md,FCC,ALBE_LC_SIC,C, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB|ATOM_ATTR_VB, 1,LCNTX,LCNTY,LCNTZ,&r,0,NULL); #else r.x=0.5*0.25*TM_LC_SIC; r.y=r.x; r.z=r.x; - create_lattice(&md,FCC,TM_LC_SIC,SI,M_SI, + create_lattice(&md,FCC,TM_LC_SIC,SI, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, 0,LCNTX,LCNTY,LCNTZ,&r,0,NULL); r.x+=0.25*TM_LC_SIC; r.y=r.x; r.z=r.x; - create_lattice(&md,FCC,TM_LC_SIC,C,M_C, + create_lattice(&md,FCC,TM_LC_SIC,C, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, 1,LCNTX,LCNTY,LCNTZ,&r,0,NULL); #endif diff --git a/visualize b/visualize index 665cb67..0b8fc61 100755 --- a/visualize +++ b/visualize @@ -39,6 +39,7 @@ bx0=""; by0=""; bz0=""; bx1=""; by1=""; bz1=""; bcr=""; clx="0"; cly="0"; clz="0"; +extra=0 # parse argv @@ -58,6 +59,7 @@ while [ "$1" ]; do bx1=$5; by1=$6; bz1=$7; shift 7;; -B) bcr=$2; shift 2;; -C) lc=$2; shift 2;; + -e) extra=1; shift 1;; *) echo "options:" echo "########" @@ -159,18 +161,21 @@ EOF # atoms if [ -n "$x0" ]; then - export x0 y0 z0 x1 y1 z1 radius + export x0 y0 z0 x1 y1 z1 radius extra cat $file | grep -v '#' | awk '\ BEGIN { x0=ENVIRON["x0"]; y0=ENVIRON["y0"]; z0=ENVIRON["z0"]; x1=ENVIRON["x1"]; y1=ENVIRON["y1"]; z1=ENVIRON["z1"]; - radius=ENVIRON["radius"]; + radius=ENVIRON["radius"]; extra=ENVIRON["extra"]; } { if(($2>=x0)&&($3>=y0)&&($4>=z0)&&\ ($2<=x1)&&($3<=y1)&&($4<=z1)) { print "sphere { <"$2","$4","$3">, "radius" "; - print "texture { pigment { color "$5" } "; + if(extra) + print "texture { pigment { color rgb<"$6/4.4",0,"1-$6/4.4"> } "; + else + print "texture { pigment { color "$5" } "; print "finish { phong 1 metallic } } }"; } }' >> temp.pov @@ -194,29 +199,30 @@ EOF # boundaries if [ -z "$bx0" ]; then - if [ -z "$x0" ]; then + #if [ -z "$x0" ]; then cat $file | grep '# \[D\]' | while read foo bar x1 y1 z1 x2 y2 z2 ; do - draw_cyl $x1 $z1 $y1 $x2 $z2 $y2 0.05 + draw_cyl $x1 $y1 $z1 $x2 $y2 $z2 0.05 done - else + #else + # manually drawing the 3x4 boundaries ... - draw_cyl $x0 $y0 $z0 $x1 $y0 $z0 - draw_cyl $x0 $y0 $z0 $x0 $y1 $z0 - draw_cyl $x1 $y1 $z0 $x1 $y0 $z0 - draw_cyl $x0 $y1 $z0 $x1 $y1 $z0 - - draw_cyl $x0 $y0 $z1 $x1 $y0 $z1 - draw_cyl $x0 $y0 $z1 $x0 $y1 $z1 - draw_cyl $x1 $y1 $z1 $x1 $y0 $z1 - draw_cyl $x0 $y1 $z1 $x1 $y1 $z1 - - draw_cyl $x0 $y0 $z1 $x0 $y0 $z0 - draw_cyl $x0 $y1 $z1 $x0 $y1 $z0 - draw_cyl $x1 $y0 $z1 $x1 $y0 $z0 - draw_cyl $x1 $y1 $z1 $x1 $y1 $z0 - fi +# draw_cyl $x0 $y0 $z0 $x1 $y0 $z0 +# draw_cyl $x0 $y0 $z0 $x0 $y1 $z0 +# draw_cyl $x1 $y1 $z0 $x1 $y0 $z0 +# draw_cyl $x0 $y1 $z0 $x1 $y1 $z0 + +# draw_cyl $x0 $y0 $z1 $x1 $y0 $z1 +# draw_cyl $x0 $y0 $z1 $x0 $y1 $z1 +# draw_cyl $x1 $y1 $z1 $x1 $y0 $z1 +# draw_cyl $x0 $y1 $z1 $x1 $y1 $z1 + +# draw_cyl $x0 $y0 $z1 $x0 $y0 $z0 +# draw_cyl $x0 $y1 $z1 $x0 $y1 $z0 +# draw_cyl $x1 $y0 $z1 $x1 $y0 $z0 +# draw_cyl $x1 $y1 $z1 $x1 $y1 $z0 +# fi else @@ -287,7 +293,7 @@ EOF # mv png $POVRAY temp.pov > /dev/null 2>&1 - time=`echo $file | awk -F. '{ print $1 }' | awk -F_ '{ print $3 }'` + time=`basename $file | awk -F. '{ print $1 }' | awk -F_ '{ print $3 }'` convert $COPTS -draw "text 5,20 't = $time fs'" temp.png temp.png mv temp.png `echo $file | sed 's/\.xyz/\.png/'`