sic.c

   1 /*
   2  * sic.c - investigation of the sic precipitation process of silicon carbide
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #include <math.h>
   9
  10 #include "moldyn.h"
  11
  12 /* potential */
  13 #include "potentials/harmonic_oscillator.h"
  14 #include "potentials/lennard_jones.h"
  15 #include "potentials/albe.h"
  16 #ifdef TERSOFF_ORIG
  17 #include "potentials/tersoff_orig.h"
  18 #else
  19 #include "potentials/tersoff.h"
  20 #endif
  21
  22 typedef struct s_hp {
  23         int prerun_count;       /* prerun count */
  24         int insert_count;       /* insert count */
  25         int postrun_count;      /* post run count */
  26         unsigned char state;    /* current state */
  27         int argc;               /* arg count */
  28         char **argv;            /* args */
  29 } t_hp;
  30
  31 #define STATE_PRERUN    0x00
  32 #define STATE_INSERT    0x01
  33 #define STATE_POSTRUN   0x02
  34
  35 /* include the config file */
  36 #include "config.h"
  37
  38 int insert_atoms(t_moldyn *moldyn) {
  39
  40         int i,j;
  41         u8 run;
  42         t_3dvec r,v,dist;
  43         double d,dmin;
  44
  45         t_atom *atom;
  46
  47         atom=moldyn->atom;
  48
  49         v.x=0; v.y=0; v.z=0;
  50
  51         for(j=0;j<INS_ATOMS;j++) {
  52                 run=1;
  53                 while(run) {
  54 #ifdef INS_TETRA
  55                         // tetrahedral
  56                         r.x=0.0;
  57                         r.y=0.0;
  58                         r.z=0.0;
  59 #endif
  60 #ifdef INS_HEXA
  61                         // hexagonal
  62                         r.x=-1.0/8.0*ALBE_LC_SI;
  63                         r.y=-1.0/8.0*ALBE_LC_SI;
  64                         r.z=1.0/8.0*ALBE_LC_SI;
  65 #endif
  66 #ifdef INS_110DB
  67                         // 110 dumbbell
  68                         r.x=(-0.5+0.25+0.125)*ALBE_LC_SI;
  69                         r.y=(-0.5+0.25+0.125)*ALBE_LC_SI;
  70                         r.z=(-0.5+0.25)*ALBE_LC_SI;
  71                         moldyn->atom[4372].r.x=(-0.5+0.125+0.125)*ALBE_LC_SI;
  72                         moldyn->atom[4372].r.y=(-0.5+0.125+0.125)*ALBE_LC_SI;
  73 #endif
  74 #ifdef INS_001DB
  75                         // 001 dumbbell
  76                         r.x=(-0.5+0.25)*ALBE_LC_SI;
  77                         r.y=(-0.5+0.25)*ALBE_LC_SI;
  78                         r.z=(-0.1)*ALBE_LC_SI;
  79                         moldyn->atom[4372].r.z=(-0.4)*ALBE_LC_SI;
  80 #endif
  81 #ifdef INS_USER
  82                         // 001 dumbbell
  83                         r.x=INS_UX*ALBE_LC_SI;
  84                         r.y=INS_UY*ALBE_LC_SI;
  85                         r.z=INS_UZ*ALBE_LC_SI;
  86 #endif
  87 #ifdef INS_RAND
  88                         // random
  89 #ifdef INS_DYNAMIC_LEN
  90                         r.x=(rand_get_double(&(moldyn->random))-0.5)*\
  91                             moldyn->dim.x;
  92                         r.y=(rand_get_double(&(moldyn->random))-0.5)*\
  93                             moldyn->dim.y;
  94                         r.z=(rand_get_double(&(moldyn->random))-0.5)*\
  95                             moldyn->dim.z;
  96 #else
  97                         r.x=(rand_get_double(&(moldyn->random))-0.5)*INS_LENX;
  98                         r.y=(rand_get_double(&(moldyn->random))-0.5)*INS_LENY;
  99                         r.z=(rand_get_double(&(moldyn->random))-0.5)*INS_LENZ;
 100 #endif
 101 #endif
 102                         // offset
 103                         r.x+=INS_OFFSET;
 104                         r.y+=INS_OFFSET;
 105                         r.z+=INS_OFFSET;
 106                         /* assume valid coordinates */
 107                         run=0;
 108                         dmin=10000000000.0;             // for sure too high!
 109                         for(i=0;i<moldyn->count;i++) {
 110                                 atom=&(moldyn->atom[i]);
 111                                 v3_sub(&dist,&(atom->r),&r);
 112                                 check_per_bound(moldyn,&dist);
 113                                 d=v3_absolute_square(&dist);
 114                                 /* reject coordinates */
 115                                 if(d<INS_R_C) {
 116                                         //printf("atom %d - %f\n",i,d);
 117                                         run=1;
 118                                         break;
 119                                 }
 120                                 if(d<dmin)
 121                                         dmin=d;
 122                         }
 123                 }
 124                 add_atom(moldyn,INS_TYPE,INS_MASS,INS_BRAND,
 125                          ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|\
 126                          INS_ATTR,
 127                          &r,&v);
 128                 printf(" %02d: atom %d | %f %f %f | %f\n",
 129                        j,moldyn->count-1,r.x,r.y,r.z,dmin);
 130         }
 131
 132         return 0;
 133 }
 134
 135 int sic_hook(void *moldyn,void *hook_params) {
 136
 137         t_hp *hp;
 138         t_moldyn *md;
 139         int steps;
 140         double tau;
 141         double dt;
 142         double dp;
 143
 144         hp=hook_params;
 145         md=moldyn;
 146
 147         tau=1.0;
 148         steps=0;
 149
 150         /* switch on t scaling */
 151         if(md->schedule.count==0)
 152                 set_pt_scale(md,P_SCALE_BERENDSEN,P_SCALE_TAU,
 153                                 T_SCALE_BERENDSEN,T_SCALE_TAU);
 154
 155         /* my lousy state machine ! */
 156
 157         /* switch to insert state immediately */
 158         if(hp->state==STATE_PRERUN)
 159                 hp->state=STATE_INSERT;
 160
 161         switch(hp->state) {
 162                 case STATE_INSERT:
 163                         goto insert;
 164                         break;
 165                 case STATE_POSTRUN:
 166                         goto postrun;
 167                         break;
 168                 default:
 169                         printf("[sic hook] unknown state\n");
 170                         return -1;
 171         }
 172
 173         /* act according to state */
 174
 175 insert:
 176
 177         /* assigne values */
 178         steps=INS_RELAX;
 179         tau=INS_TAU;
 180
 181         /* check temperature */
 182         dt=md->t_avg-md->t_ref;
 183         dp=md->p_avg-md->p_ref;
 184         if(dt<0)
 185                 dt=-dt;
 186         if(dp<0)
 187                 dp=-dp;
 188         if((dt>INS_DELTA_TC)|(dp>INS_DELTA_PC))
 189                 goto addsched;
 190
 191         /* immediately go on if no job is to be done */
 192         if(hp->insert_count==INS_RUNS) {
 193                 printf("    --- leaving insert state ---\n");
 194                 hp->state=STATE_POSTRUN;
 195                 goto postrun;
 196         }
 197
 198         /* else -> insert atoms */
 199         hp->insert_count+=1;
 200         printf("   ### insert atoms (%d/%d) ###\n",
 201                hp->insert_count*INS_ATOMS,INS_RUNS*INS_ATOMS);
 202         insert_atoms(md);
 203         goto addsched;
 204
 205 postrun:
 206
 207         /* assigne values */
 208         steps=POST_RELAX;
 209         tau=POST_TAU;
 210
 211         /* check temperature */
 212         dt=md->t_avg-md->t_ref;
 213         dp=md->p_avg-md->p_ref;
 214         if(dt<0)
 215                 dt=-dt;
 216         if(dp<0)
 217                 dp=-dp;
 218         if((dt>POST_DELTA_TC)|(dp>POST_DELTA_PC))
 219                 goto addsched;
 220
 221         /* immediately return if no job is to be done */
 222         if(hp->postrun_count==POST_RUNS) {
 223                 printf("    --- leaving post run state ---\n");
 224                 return 0;
 225         }
 226
 227         /* postrun action */
 228         hp->postrun_count+=1;
 229         printf(" ### postrun (%d/%d) ###\n",
 230                hp->postrun_count,POST_RUNS);
 231         set_temperature(md,md->t_ref-POST_DT);
 232
 233 addsched:
 234
 235         /* reset the average counters */
 236         average_reset(md);
 237
 238         /* add schedule */
 239         moldyn_add_schedule(md,steps,tau);
 240
 241         return 0;
 242 }
 243
 244 int main(int argc,char **argv) {
 245
 246         /* main moldyn structure */
 247         t_moldyn md;
 248
 249         /* hook parameter structure */
 250         t_hp hookparam;
 251
 252         /* testing location & velocity vector */
 253         t_3dvec r,v;
 254         memset(&r,0,sizeof(t_3dvec));
 255         memset(&v,0,sizeof(t_3dvec));
 256
 257         /* initialize moldyn */
 258         moldyn_init(&md,argc,argv);
 259
 260         /* choose integration algorithm */
 261         set_int_alg(&md,MOLDYN_INTEGRATE_VERLET);
 262
 263         /* choose potential */
 264 #ifdef ALBE
 265         if(set_potential(&md,MOLDYN_POTENTIAL_AM)<0)
 266                 return -1;
 267 #else
 268         if(set_potential(&md,MOLDYN_POTENTIAL_TM)<0)
 269                 return -1;
 270 #endif
 271
 272         /* cutoff radius & bondlen */
 273 #ifdef ALBE
 274         set_cutoff(&md,ALBE_S_SI);
 275         set_bondlen(&md,ALBE_S_SI,ALBE_S_C,ALBE_S_SIC);
 276         //set_cutoff(&md,ALBE_S_C);
 277 #else
 278         set_cutoff(&md,TM_S_SI);
 279         set_bondlen(&md,TM_S_SI,TM_S_C,-1.0);
 280         //set_cutoff(&md,TM_S_C);
 281 #endif
 282
 283         /*
 284          * potential parameters
 285          */
 286
 287 #ifndef ALBE
 288         /*
 289          * tersoff mult potential parameters for SiC
 290          */
 291         tersoff_mult_set_params(&md,SI,C);
 292 #else
 293         /*
 294          * albe mult potential parameters for SiC
 295          */
 296         albe_mult_set_params(&md,SI,C);
 297 #endif
 298
 299         /* set (initial) dimensions of simulation volume */
 300 #ifdef ALBE
 301  #ifdef INIT_SI
 302         set_dim(&md,LCNTX*ALBE_LC_SI,LCNTY*ALBE_LC_SI,LCNTZ*ALBE_LC_SI,TRUE);
 303  #endif
 304  #ifdef INIT_C
 305         set_dim(&md,LCNTX*ALBE_LC_C,LCNTY*ALBE_LC_C,LCNTZ*ALBE_LC_C,TRUE);
 306  #endif
 307  #ifdef INIT_3CSIC
 308         set_dim(&md,LCNTX*ALBE_LC_SIC,LCNTY*ALBE_LC_SIC,LCNTZ*ALBE_LC_SIC,TRUE);
 309  #endif
 310 #else
 311  #ifdef INIT_SI
 312         set_dim(&md,LCNTX*LC_SI,LCNTY*LC_SI,LCNTZ*LC_SI,TRUE);
 313  #endif
 314  #ifdef INIT_C
 315         set_dim(&md,LCNTX*LC_C,LCNTY*LC_C,LCNTZ*LC_C,TRUE);
 316  #endif
 317  #ifdef INIT_3CSIC
 318         set_dim(&md,LCNTX*TM_LC_SIC,LCNTY*TM_LC_SIC,LCNTZ*TM_LC_SIC,TRUE);
 319  #endif
 320 #endif
 321
 322         /* set periodic boundary conditions in all directions */
 323         set_pbc(&md,TRUE,TRUE,TRUE);
 324
 325         /* create the lattice / place atoms */
 326
 327         // diamond
 328 #ifdef ALBE
 329  #ifdef INIT_SI
 330         create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,M_SI,
 331                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 332                        0,LCNTX,LCNTY,LCNTZ,NULL);
 333  #endif
 334  #ifdef INIT_C
 335         create_lattice(&md,DIAMOND,ALBE_LC_C,C,M_C,
 336                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 337                        1,LCNTX,LCNTY,LCNTZ,NULL);
 338  #endif
 339 #else
 340  #ifdef INIT_SI
 341         create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
 342                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 343                        0,LCNTX,LCNTY,LCNTZ,NULL);
 344  #endif
 345  #ifdef INIT_C
 346         create_lattice(&md,DIAMOND,LC_C,SI,M_SI,
 347                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 348                        1,LCNTX,LCNTY,LCNTZ,NULL);
 349  #endif
 350 #endif
 351
 352         // zinkblende
 353 #ifdef INIT_3CSIC
 354  #ifdef ALBE
 355         r.x=0.5*0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
 356         create_lattice(&md,FCC,ALBE_LC_SIC,SI,M_SI,
 357                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 358                        0,LCNTX,LCNTY,LCNTZ,&r);
 359         r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
 360         create_lattice(&md,FCC,ALBE_LC_SIC,C,M_C,
 361                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB|ATOM_ATTR_VB,
 362                        1,LCNTX,LCNTY,LCNTZ,&r);
 363  #else
 364         r.x=0.5*0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
 365         create_lattice(&md,FCC,TM_LC_SIC,SI,M_SI,
 366                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 367                        0,LCNTX,LCNTY,LCNTZ,&r);
 368         r.x+=0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
 369         create_lattice(&md,FCC,TM_LC_SIC,C,M_C,
 370                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 371                        1,LCNTX,LCNTY,LCNTZ,&r);
 372  #endif
 373 #endif
 374
 375         /* check for right atom placing */
 376         moldyn_bc_check(&md);
 377
 378         /* testing configuration */
 379         //r.x=0.27*sqrt(3.0)*LC_SI/2.0; v.x=0;
 380         //r.x=(TM_S_SI+TM_R_SI)/4.0;    v.x=0;
 381         //r.y=0;                v.y=0;
 382         //r.z=0;                v.z=0;
 383         //add_atom(&md,SI,M_SI,0,
 384         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 385         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
 386         //           &r,&v);
 387         //r.x=-r.x;     v.x=-v.x;
 388         //r.y=0;                v.y=0;
 389         //r.z=0;                v.z=0;
 390         //add_atom(&md,SI,M_SI,0,
 391         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 392         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
 393         //           &r,&v);
 394         //r.z=0.27*sqrt(3.0)*LC_SI/2.0; v.z=0;
 395         //r.x=(TM_S_SI+TM_R_SI)/4.0;    v.x=0;
 396         //r.y=0;                v.y=0;
 397         //r.x=0;                v.x=0;
 398         //add_atom(&md,SI,M_SI,0,
 399         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 400         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
 401         //           &r,&v);
 402         //r.z=-r.z;     v.z=-v.z;
 403         //r.y=0;                v.y=0;
 404         //r.x=0;                v.x=0;
 405         //add_atom(&md,SI,M_SI,0,
 406         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 407         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
 408         //           &r,&v);
 409
 410         /* set temperature & pressure */
 411         set_temperature(&md,atof(argv[2])+273.0);
 412         set_pressure(&md,0.0);
 413
 414         /* set amount of steps to skip before average calc */
 415         set_avg_skip(&md,AVG_SKIP);
 416
 417         /* set p/t scaling */
 418         //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
 419         //set_pt_scale(&md,P_SCALE_BERENDSEN,0.01/(100*GPA),
 420         //                 T_SCALE_BERENDSEN,100.0);
 421         //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
 422         //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,0,0);
 423
 424         /* initial thermal fluctuations of particles (in equilibrium) */
 425         thermal_init(&md,TRUE);
 426
 427         /* create the simulation schedule */
 428         moldyn_add_schedule(&md,PRERUN,PRE_TAU);
 429
 430         /* schedule hook function */
 431         memset(&hookparam,0,sizeof(t_hp));
 432         hookparam.argc=argc;
 433         hookparam.argv=argv;
 434         moldyn_set_schedule_hook(&md,&sic_hook,&hookparam);
 435         //moldyn_set_schedule_hook(&md,&hook_del_atom,&hookparam);
 436         //moldyn_add_schedule(&md,POSTRUN,1.0);
 437
 438         /* activate logging */
 439         moldyn_set_log_dir(&md,argv[1]);
 440         moldyn_set_report(&md,"Frank Zirkelbach",R_TITLE);
 441         moldyn_set_log(&md,LOG_TOTAL_ENERGY,LOG_E);
 442         moldyn_set_log(&md,LOG_TEMPERATURE,LOG_T);
 443         moldyn_set_log(&md,LOG_PRESSURE,LOG_P);
 444         moldyn_set_log(&md,LOG_VOLUME,LOG_V);
 445         moldyn_set_log(&md,VISUAL_STEP,LOG_A);
 446         moldyn_set_log(&md,SAVE_STEP,LOG_S);
 447         moldyn_set_log(&md,CREATE_REPORT,0);
 448
 449         /* next neighbour distance for critical checking */
 450         set_nn_dist(&md,0.25*ALBE_LC_SI*sqrt(3.0));
 451
 452         /*
 453          * let's do the actual md algorithm now
 454          *
 455          * integration of newtons equations
 456          */
 457         moldyn_integrate(&md);
 458 #ifdef dEBUG
 459 return 0;
 460 #endif
 461
 462         /*
 463          * post processing the data
 464          */
 465
 466         /* close */
 467         moldyn_shutdown(&md);
 468
 469         return 0;
 470 }
 471