From a16579f010cb5d53b7bf44e854089019ab9b1df0 Mon Sep 17 00:00:00 2001 From: hackbard Date: Wed, 28 Jan 2009 23:09:25 +0100 Subject: [PATCH] added del_atoms, playing with Ef calc, some more TODO --- TODO.txt | 25 +++++++----------- calc_delta_e | 24 ++++++++++++++++- mdrun.c | 42 ++++++++++++++++++++++++++++++ mdrun.h | 6 +++++ moldyn.c | 73 ++++++++++++++++++++++++++++++++++++---------------- 5 files changed, 131 insertions(+), 39 deletions(-) diff --git a/TODO.txt b/TODO.txt index 8712c9d..eb0de57 100644 --- a/TODO.txt +++ b/TODO.txt @@ -1,7 +1,10 @@ -interstitial simulations -************************ +parcas +****** +more DOUBLECHECKs with PARCAS !!!! -- more DOUBLECHECKs with PARCAS !!!! + * E_kin + * x/y/z(t) of both + * video of hex -> 100db implement ********* @@ -14,7 +17,6 @@ implement - clean up mdrun - config sanity checks - - introduce/improve fill command (multiple fills) - anneal -> from current to T with rate R @@ -36,11 +38,8 @@ implement - potentials: - * tersoff fast *** - * modifief tersoff/albe *** - * stillinger weber *** - * eam - * edip + * tersoff fast + * meam - optimize code! @@ -65,22 +64,16 @@ implement todo: - listen ! estimate time - more in fast (also think about -> 'type') - - algo! (though i don't believe!) - inline - - openmp (num_thread copies of neighbour_i) - + pthreads (simultaneously!) simulation runs *************** -- cryst simulations: reasonable pctrl? -- even higher temperatures - tctrl only in outer regions - only 1 atom per timestep -- EXTENDED C-C cutoff -- different sized SiC prec in Si (4:5! diff temperatures) *** - melting exps (both, anneal + interface method) - interstitials: + - sic ints (compare!) - more interstitial combinations - characterize interstitials by PCF, then inc temperature diff --git a/calc_delta_e b/calc_delta_e index 455ce57..e0ee86f 100755 --- a/calc_delta_e +++ b/calc_delta_e @@ -1,9 +1,31 @@ #!/bin/bash +# bulk chemical potentials +musi=-4.63 +muc=-7.374 +music=-12.68 + +# formation enthalpy +dhf=0.68 +#dhf=-9.64 + +# first method to calculate formation energy file=`ls $1/atomic_conf_* | tail -1` atom_cnt=`grep '# \[P\]' $file | awk '{ print $3 }'` e0=`awk 'NR==2' $1/energy | awk '{ print $4 }'` e1=`tail -n 1 $1/energy | awk '{ print $4 }'` echo "$e0 $e1 $atom_cnt" | \ - awk '{ print ($2-$1)*$3" eV" }' + awk '{ print " "($2-$1)*$3" eV (simple)" }' + +# second method to calculate formation energy +si_cnt=`grep ^Si $file | wc -l` +c_cnt=`grep ^C $file | wc -l` +ed=`tail -n 1 $1/energy | awk '{ print $3 }'` + +echo "$si_cnt $c_cnt $musi $muc $ed $music $dhf" | \ + awk '{ print " "$5*($1+$2)-0.5*($1+$2)*$6-0.5*($1-$2)*($3-$4)-0.5*($1-$2)*$7 " eV (Albe, #Si=" $1 " #C=" $2 ") " }' +echo "$si_cnt $c_cnt $ed $e0" | \ + awk '{ print " " $3*($1+$2)-($1+$2)*$4 " eV (Gao, #Si=" $1 " #C=" $2 ")" }' + + diff --git a/mdrun.c b/mdrun.c index 34807e7..bd70d00 100644 --- a/mdrun.c +++ b/mdrun.c @@ -98,6 +98,9 @@ int add_stage(t_mdrun *mdrun,u8 type,void *params) { case STAGE_DISPLACE_ATOM: psize=sizeof(t_displace_atom_params); break; + case STAGE_DEL_ATOMS: + psize=sizeof(t_del_atoms_params); + break; case STAGE_INSERT_ATOMS: psize=sizeof(t_insert_atoms_params); break; @@ -175,6 +178,7 @@ int mdrun_parse_config(t_mdrun *mdrun) { t_set_temp_params stp; t_set_timestep_params stsp; t_fill_params fp; + t_del_atoms_params delp; /* open config file */ fd=open(mdrun->cfile,O_RDONLY); @@ -207,6 +211,7 @@ int mdrun_parse_config(t_mdrun *mdrun) { memset(&stp,0,sizeof(t_set_temp_params)); memset(&stsp,0,sizeof(t_set_timestep_params)); memset(&fp,0,sizeof(t_fill_params)); + memset(&delp,0,sizeof(t_del_atoms_params)); // get command + args wcnt=0; @@ -509,6 +514,13 @@ int mdrun_parse_config(t_mdrun *mdrun) { dap.dz=atof(word[5]); add_stage(mdrun,STAGE_DISPLACE_ATOM,&dap); } + else if(!strncmp(word[1],"del_atoms",9)) { + delp.o.x=atof(word[2]); + delp.o.y=atof(word[3]); + delp.o.z=atof(word[4]); + delp.r=atof(word[5]); + add_stage(mdrun,STAGE_DEL_ATOMS,&delp); + } else if(!strncmp(word[1],"ins_atoms",9)) { iap.ins_steps=atoi(word[2]); iap.ins_atoms=atoi(word[3]); @@ -732,6 +744,31 @@ int displace_atom(t_moldyn *moldyn,t_mdrun *mdrun) { return 0; } +int del_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { + + t_stage *stage; + t_del_atoms_params *delp; + int i; + t_3dvec dist; + + stage=mdrun->stage.current->data; + delp=stage->params; + + for(i=0;icount;i++) { + v3_sub(&dist,&(delp->o),&(moldyn->atom[i].r)); +//printf("%d ----> %f %f %f = %f | %f\n",i,dist.x,dist.y,dist.z,v3_absolute_square(&dist),delp->r*delp->r); + if(v3_absolute_square(&dist)<=(delp->r*delp->r)) { + del_atom(moldyn,moldyn->atom[i].tag); + printf("%s atom deleted: %d %d %d\n",ME, + moldyn->atom[i].tag,moldyn->atom[i].element, + moldyn->atom[i].brand); + } + } + + return 0; + +} + int insert_atoms(t_moldyn *moldyn,t_mdrun *mdrun) { t_insert_atoms_params *iap; @@ -1108,6 +1145,11 @@ int mdrun_hook(void *ptr1,void *ptr2) { displace_atom(moldyn,mdrun); change_stage=TRUE; break; + case STAGE_DEL_ATOMS: + stage_print(" -> del atoms\n\n"); + del_atoms(moldyn,mdrun); + change_stage=TRUE; + break; case STAGE_INSERT_ATOMS: stage_print(" -> insert atoms\n\n"); iap=stage->params; diff --git a/mdrun.h b/mdrun.h index 3f2c6ec..be09d5a 100644 --- a/mdrun.h +++ b/mdrun.h @@ -52,6 +52,7 @@ typedef struct s_stage { #define STAGE_SET_TIMESTEP 0x08 #define STAGE_FILL 0x09 #define STAGE_THERMAL_INIT 0x10 +#define STAGE_DEL_ATOMS 0x11 typedef struct s_mdrun { char cfile[128]; // config file @@ -108,6 +109,11 @@ typedef struct s_displace_atom_params { double dx,dy,dz; } t_displace_atom_params; +typedef struct s_del_atoms_params { + double r; + t_3dvec o; +} t_del_atoms_params; + typedef struct s_insert_atoms_params { u8 type; double x0,y0,z0,x1,y1,z1; diff --git a/moldyn.c b/moldyn.c index bc4fdff..7698c24 100644 --- a/moldyn.c +++ b/moldyn.c @@ -736,6 +736,9 @@ int del_atom(t_moldyn *moldyn,int tag) { t_atom *new,*old; int cnt; +#if defined LOWMEM_LISTS || defined PTHREADS + void *ptr; +#endif old=moldyn->atom; @@ -758,6 +761,28 @@ int del_atom(t_moldyn *moldyn,int tag) { free(old); +#ifdef LOWMEM_LISTS + ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int)); + if(!ptr) { + perror("[moldyn] list realloc (del atom)"); + return -1; + } + moldyn->lc.subcell->list=ptr; + // update lists + link_cell_update(moldyn); +#endif + +#ifdef PTHREADS + ptr=realloc(amutex,moldyn->count*sizeof(pthread_mutex_t)); + if(!ptr) { + perror("[moldyn] mutex realloc (add atom)"); + return -1; + } + amutex=ptr; + pthread_mutex_destroy(&(amutex[moldyn->count+1])); +#endif + + return 0; } @@ -2065,10 +2090,11 @@ int moldyn_integrate(t_moldyn *moldyn) { if(e) { if(!(moldyn->total_steps%e)) dprintf(moldyn->efd, - "%f %f %f %f\n", + "%f %f %f %f %f %f\n", moldyn->time,moldyn->ekin/energy_scale, moldyn->energy/energy_scale, - get_total_energy(moldyn)/energy_scale); + get_total_energy(moldyn)/energy_scale, + moldyn->ekin/EV,moldyn->energy/EV); } if(m) { if(!(moldyn->total_steps%m)) { @@ -3203,9 +3229,8 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, ba=data; /* increase total bond counter - * ... double counting! */ - ba->tcnt+=2; + ba->tcnt+=1; if(itom->brand==0) ba->acnt[jtom->tag]+=1; @@ -3222,10 +3247,11 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, int bond_analyze(t_moldyn *moldyn,double *quality) { - // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total - - int qcnt; - int ccnt,cset; + int qcnt4; + int qcnt3; + int ccnt4; + int ccnt3; + int bcnt; t_ba ba; int i; t_atom *atom; @@ -3245,9 +3271,9 @@ int bond_analyze(t_moldyn *moldyn,double *quality) { memset(ba.bcnt,0,moldyn->count*sizeof(int)); ba.tcnt=0; - qcnt=0; - ccnt=0; - cset=0; + qcnt4=0; qcnt3=0; + ccnt4=0; ccnt3=0; + bcnt=0; atom=moldyn->atom; @@ -3256,27 +3282,30 @@ int bond_analyze(t_moldyn *moldyn,double *quality) { for(i=0;icount;i++) { if(atom[i].brand==0) { if((ba.acnt[i]==0)&(ba.bcnt[i]==4)) - qcnt+=4; + qcnt4+=4; + if((ba.acnt[i]==0)&(ba.bcnt[i]==3)) + qcnt3+=3; } else { if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) { - qcnt+=4; - ccnt+=1; + qcnt4+=4; + ccnt4+=1; } - cset+=1; + if((ba.acnt[i]==3)&(ba.bcnt[i]==0)) { + qcnt3+=4; + ccnt3+=1; + } + bcnt+=1; } } - printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset); - printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt); - if(quality) { - quality[0]=1.0*ccnt/cset; - quality[1]=1.0*qcnt/ba.tcnt; + quality[0]=1.0*ccnt4/bcnt; + quality[1]=1.0*ccnt3/bcnt; } else { - printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt); - printf("[moldyn] bond analyze: tot_q=%f\n",1.0*qcnt/ba.tcnt); + printf("[moldyn] bond analyze: %f %f\n", + 1.0*ccnt4/bcnt,1.0*ccnt3/bcnt); } return 0; -- 2.20.1