added del_atoms, playing with Ef calc, some more TODO

diff --git a/TODO.txt b/TODO.txt
index 8712c9d..eb0de57 100644 (file)
--- 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 (executable)
--- 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 (file)
--- 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;i<moldyn->count;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 (file)
--- 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 (file)
--- 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;i<moldyn->count;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;