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                         md->atom[4372].r.x=(-0.5+0.125+0.125)*ALBE_LC_SI;
  72                         md->atom[4372].r.y=(-0.5+0.125+0.125)*ALBE_LC_SI;
  73 #endif
  74 #ifdef INS_RAND
  75                         // random
  76                         r.x=(rand_get_double(&(moldyn->random))-0.5)*INS_LENX;
  77                         r.y=(rand_get_double(&(moldyn->random))-0.5)*INS_LENY;
  78                         r.z=(rand_get_double(&(moldyn->random))-0.5)*INS_LENZ;
  79 #endif
  80                         // offset
  81                         r.x+=INS_OFFSET;
  82                         r.y+=INS_OFFSET;
  83                         r.z+=INS_OFFSET;
  84                         /* assume valid coordinates */
  85                         run=0;
  86                         dmin=10000000000.0;             // for sure too high!
  87                         for(i=0;i<moldyn->count;i++) {
  88                                 atom=&(moldyn->atom[i]);
  89                                 v3_sub(&dist,&(atom->r),&r);
  90                                 check_per_bound(moldyn,&dist);
  91                                 d=v3_absolute_square(&dist);
  92                                 /* reject coordinates */
  93                                 if(d<INS_R_C) {
  94                                         //printf("atom %d - %f\n",i,d);
  95                                         run=1;
  96                                         break;
  97                                 }
  98                                 if(d<dmin)
  99                                         dmin=d;
 100                         }
 101                 }
 102                 add_atom(moldyn,INS_TYPE,INS_MASS,INS_BRAND,
 103                          ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|\
 104                          INS_ATTR,
 105                          &r,&v);
 106                 printf(" %02d: atom %d | %f %f %f | %f\n",
 107                        j,moldyn->count-1,r.x,r.y,r.z,dmin);
 108         }
 109
 110         return 0;
 111 }
 112
 113 int sic_hook(void *moldyn,void *hook_params) {
 114
 115         t_hp *hp;
 116         t_moldyn *md;
 117         int steps;
 118         double tau;
 119         double dt;
 120
 121         hp=hook_params;
 122         md=moldyn;
 123
 124         tau=1.0;
 125         steps=0;
 126
 127         /* switch on t scaling */
 128         if(md->schedule.count==0)
 129                 set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0);
 130
 131         /* my lousy state machine ! */
 132
 133         /* switch to insert state immediately */
 134         if(hp->state==STATE_PRERUN)
 135                 hp->state=STATE_INSERT;
 136
 137         /* act according to state */
 138         switch(hp->state) {
 139                 case STATE_INSERT:
 140                         /* assigne values */
 141                         steps=INS_RELAX;
 142                         tau=INS_TAU;
 143                         /* check temperature */
 144                         dt=md->t_avg-md->t_ref;
 145                         if(dt<0)
 146                                 dt=-dt;
 147                         if(dt>INS_DELTA_TC)
 148                                 break;
 149                         /* insert atoms */
 150                         hp->insert_count+=1;
 151                         printf("   ### insert atoms (%d/%d) ###\n",
 152                                hp->insert_count*INS_ATOMS,INS_RUNS*INS_ATOMS);
 153                         insert_atoms(md);
 154                         /* change state after last insertion */
 155                         if(hp->insert_count==INS_RUNS)
 156                                 hp->state=STATE_POSTRUN;
 157                         break;
 158                 case STATE_POSTRUN:
 159                         /* assigne values */
 160                         steps=POST_RELAX;
 161                         tau=POST_TAU;
 162                         /* check temperature */
 163                         dt=md->t_avg-md->t_ref;
 164                         if(dt<0)
 165                                 dt=-dt;
 166                         if(dt>POST_DELTA_TC)
 167                                 break;
 168                         /* decrease temperature */
 169                         hp->postrun_count+=1;
 170                         printf(" ### postrun (%d/%d) ###\n",
 171                                hp->postrun_count,POST_RUNS);
 172                         set_temperature(md,md->t_ref-POST_DT);
 173                         if(hp->postrun_count==POST_RUNS)
 174                                 return 0;
 175                         break;
 176                 default:
 177                         printf("[hook] FATAL (default case!?!)\n");
 178                         break;
 179         }
 180
 181         /* reset the average counters */
 182         average_reset(md);
 183
 184         /* add schedule */
 185         moldyn_add_schedule(md,steps,tau);
 186
 187         return 0;
 188 }
 189
 190 int main(int argc,char **argv) {
 191
 192         /* main moldyn structure */
 193         t_moldyn md;
 194
 195         /* hook parameter structure */
 196         t_hp hookparam;
 197
 198         /* potential parameters */
 199         t_tersoff_mult_params tp;
 200         t_albe_mult_params ap;
 201
 202         /* testing location & velocity vector */
 203         t_3dvec r,v;
 204         memset(&r,0,sizeof(t_3dvec));
 205         memset(&v,0,sizeof(t_3dvec));
 206
 207         /* initialize moldyn */
 208         moldyn_init(&md,argc,argv);
 209
 210         /* choose integration algorithm */
 211         set_int_alg(&md,MOLDYN_INTEGRATE_VERLET);
 212
 213         /* choose potential */
 214 #ifdef ALBE
 215         set_potential3b_j1(&md,albe_mult_3bp_j1);
 216         set_potential3b_k1(&md,albe_mult_3bp_k1);
 217         set_potential3b_j2(&md,albe_mult_3bp_j2);
 218         set_potential3b_k2(&md,albe_mult_3bp_k2);
 219 #else
 220         set_potential1b(&md,tersoff_mult_1bp);
 221         set_potential3b_j1(&md,tersoff_mult_3bp_j1);
 222         set_potential3b_k1(&md,tersoff_mult_3bp_k1);
 223         set_potential3b_j2(&md,tersoff_mult_3bp_j2);
 224         set_potential3b_k2(&md,tersoff_mult_3bp_k2);
 225 #endif
 226
 227 #ifdef ALBE
 228         set_potential_params(&md,&ap);
 229 #else
 230         set_potential_params(&md,&tp);
 231 #endif
 232
 233         /* cutoff radius & bondlen */
 234 #ifdef ALBE
 235         set_cutoff(&md,ALBE_S_SI);
 236         set_bondlen(&md,ALBE_S_SI,ALBE_S_C,ALBE_S_SIC);
 237         //set_cutoff(&md,ALBE_S_C);
 238 #else
 239         set_cutoff(&md,TM_S_SI);
 240         set_bondlen(&md,TM_S_SI,TM_S_C,-1.0);
 241         //set_cutoff(&md,TM_S_C);
 242 #endif
 243
 244         /*
 245          * potential parameters
 246          */
 247
 248         /*
 249          * tersoff mult potential parameters for SiC
 250          */
 251         tp.S[0]=TM_S_SI;
 252         tp.R[0]=TM_R_SI;
 253         tp.A[0]=TM_A_SI;
 254         tp.B[0]=TM_B_SI;
 255         tp.lambda[0]=TM_LAMBDA_SI;
 256         tp.mu[0]=TM_MU_SI;
 257         tp.beta[0]=TM_BETA_SI;
 258         tp.n[0]=TM_N_SI;
 259         tp.c[0]=TM_C_SI;
 260         tp.d[0]=TM_D_SI;
 261         tp.h[0]=TM_H_SI;
 262
 263         tp.S[1]=TM_S_C;
 264         tp.R[1]=TM_R_C;
 265         tp.A[1]=TM_A_C;
 266         tp.B[1]=TM_B_C;
 267         tp.lambda[1]=TM_LAMBDA_C;
 268         tp.mu[1]=TM_MU_C;
 269         tp.beta[1]=TM_BETA_C;
 270         tp.n[1]=TM_N_C;
 271         tp.c[1]=TM_C_C;
 272         tp.d[1]=TM_D_C;
 273         tp.h[1]=TM_H_C;
 274
 275         tp.chi=TM_CHI_SIC;
 276
 277         tersoff_mult_complete_params(&tp);
 278
 279         /*
 280          * albe mult potential parameters for SiC
 281          */
 282         ap.S[0]=ALBE_S_SI;
 283         ap.R[0]=ALBE_R_SI;
 284         ap.A[0]=ALBE_A_SI;
 285         ap.B[0]=ALBE_B_SI;
 286         ap.r0[0]=ALBE_R0_SI;
 287         ap.lambda[0]=ALBE_LAMBDA_SI;
 288         ap.mu[0]=ALBE_MU_SI;
 289         ap.gamma[0]=ALBE_GAMMA_SI;
 290         ap.c[0]=ALBE_C_SI;
 291         ap.d[0]=ALBE_D_SI;
 292         ap.h[0]=ALBE_H_SI;
 293
 294         ap.S[1]=ALBE_S_C;
 295         ap.R[1]=ALBE_R_C;
 296         ap.A[1]=ALBE_A_C;
 297         ap.B[1]=ALBE_B_C;
 298         ap.r0[1]=ALBE_R0_C;
 299         ap.lambda[1]=ALBE_LAMBDA_C;
 300         ap.mu[1]=ALBE_MU_C;
 301         ap.gamma[1]=ALBE_GAMMA_C;
 302         ap.c[1]=ALBE_C_C;
 303         ap.d[1]=ALBE_D_C;
 304         ap.h[1]=ALBE_H_C;
 305
 306         ap.Smixed=ALBE_S_SIC;
 307         ap.Rmixed=ALBE_R_SIC;
 308         ap.Amixed=ALBE_A_SIC;
 309         ap.Bmixed=ALBE_B_SIC;
 310         ap.r0_mixed=ALBE_R0_SIC;
 311         ap.lambda_m=ALBE_LAMBDA_SIC;
 312         ap.mu_m=ALBE_MU_SIC;
 313         ap.gamma_m=ALBE_GAMMA_SIC;
 314         ap.c_mixed=ALBE_C_SIC;
 315         ap.d_mixed=ALBE_D_SIC;
 316         ap.h_mixed=ALBE_H_SIC;
 317
 318         albe_mult_complete_params(&ap);
 319
 320         /* set (initial) dimensions of simulation volume */
 321 #ifdef ALBE
 322  #ifdef INIT_SI
 323         set_dim(&md,LCNTX*ALBE_LC_SI,LCNTY*ALBE_LC_SI,LCNTZ*ALBE_LC_SI,TRUE);
 324  #endif
 325  #ifdef INIT_C
 326         set_dim(&md,LCNTX*ALBE_LC_C,LCNTY*ALBE_LC_C,LCNTZ*ALBE_LC_C,TRUE);
 327  #endif
 328  #ifdef INIT_3CSIC
 329         set_dim(&md,LCNTX*ALBE_LC_SIC,LCNTY*ALBE_LC_SIC,LCNTZ*ALBE_LC_SIC,TRUE);
 330  #endif
 331 #else
 332  #ifdef INIT_SI
 333         set_dim(&md,LCNTX*LC_SI,LCNTY*LC_SI,LCNTZ*LC_SI,TRUE);
 334  #endif
 335  #ifdef INIT_C
 336         set_dim(&md,LCNTX*LC_C,LCNTY*LC_C,LCNTZ*LC_C,TRUE);
 337  #endif
 338  #ifdef INIT_3CSIC
 339         set_dim(&md,LCNTX*TM_LC_SIC,LCNTY*TM_LC_SIC,LCNTZ*TM_LC_SIC,TRUE);
 340  #endif
 341 #endif
 342
 343         /* set periodic boundary conditions in all directions */
 344         set_pbc(&md,TRUE,TRUE,TRUE);
 345
 346         /* create the lattice / place atoms */
 347
 348         // diamond
 349 #ifdef ALBE
 350  #ifdef INIT_SI
 351         create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,M_SI,
 352                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 353                        0,LCNTX,LCNTY,LCNTZ,NULL);
 354  #endif
 355  #ifdef INIT_C
 356         create_lattice(&md,DIAMOND,ALBE_LC_C,C,M_C,
 357                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 358                        1,LCNTX,LCNTY,LCNTZ,NULL);
 359  #endif
 360 #else
 361  #ifdef INIT_SI
 362         create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
 363                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 364                        0,LCNTX,LCNTY,LCNTZ,NULL);
 365  #endif
 366  #ifdef INIT_C
 367         create_lattice(&md,DIAMOND,LC_C,SI,M_SI,
 368                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 369                        1,LCNTX,LCNTY,LCNTZ,NULL);
 370  #endif
 371 #endif
 372
 373         // zinkblende
 374 #ifdef INIT_3CSIC
 375  #ifdef ALBE
 376         r.x=0.5*0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
 377         create_lattice(&md,FCC,ALBE_LC_SIC,SI,M_SI,
 378                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 379                        0,LCNTX,LCNTY,LCNTZ,&r);
 380         r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
 381         create_lattice(&md,FCC,ALBE_LC_SIC,C,M_C,
 382                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 383                        1,LCNTX,LCNTY,LCNTZ,&r);
 384  #else
 385         r.x=0.5*0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
 386         create_lattice(&md,FCC,TM_LC_SIC,SI,M_SI,
 387                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 388                        0,LCNTX,LCNTY,LCNTZ,&r);
 389         r.x+=0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
 390         create_lattice(&md,FCC,TM_LC_SIC,C,M_C,
 391                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 392                        1,LCNTX,LCNTY,LCNTZ,&r);
 393  #endif
 394 #endif
 395
 396         /* check for right atom placing */
 397         moldyn_bc_check(&md);
 398
 399         /* testing configuration */
 400         //r.x=0.27*sqrt(3.0)*LC_SI/2.0; v.x=0;
 401         //r.x=(TM_S_SI+TM_R_SI)/4.0;    v.x=0;
 402         //r.y=0;                v.y=0;
 403         //r.z=0;                v.z=0;
 404         //add_atom(&md,SI,M_SI,0,
 405         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 406         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
 407         //           &r,&v);
 408         //r.x=-r.x;     v.x=-v.x;
 409         //r.y=0;                v.y=0;
 410         //r.z=0;                v.z=0;
 411         //add_atom(&md,SI,M_SI,0,
 412         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 413         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
 414         //           &r,&v);
 415         //r.z=0.27*sqrt(3.0)*LC_SI/2.0; v.z=0;
 416         //r.x=(TM_S_SI+TM_R_SI)/4.0;    v.x=0;
 417         //r.y=0;                v.y=0;
 418         //r.x=0;                v.x=0;
 419         //add_atom(&md,SI,M_SI,0,
 420         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 421         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
 422         //           &r,&v);
 423         //r.z=-r.z;     v.z=-v.z;
 424         //r.y=0;                v.y=0;
 425         //r.x=0;                v.x=0;
 426         //add_atom(&md,SI,M_SI,0,
 427         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
 428         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
 429         //           &r,&v);
 430
 431         /* set temperature & pressure */
 432         set_temperature(&md,atof(argv[2])+273.0);
 433         set_pressure(&md,BAR);
 434
 435         /* set amount of steps to skip before average calc */
 436         set_avg_skip(&md,AVG_SKIP);
 437
 438         /* set p/t scaling */
 439         //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
 440         //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,
 441         //                 T_SCALE_BERENDSEN,100.0);
 442         //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
 443         //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,0,0);
 444
 445         /* initial thermal fluctuations of particles (in equilibrium) */
 446         thermal_init(&md,TRUE);
 447
 448         /* create the simulation schedule */
 449         moldyn_add_schedule(&md,PRERUN,PRE_TAU);
 450
 451         /* schedule hook function */
 452         memset(&hookparam,0,sizeof(t_hp));
 453         hookparam.argc=argc;
 454         hookparam.argv=argv;
 455         moldyn_set_schedule_hook(&md,&sic_hook,&hookparam);
 456         //moldyn_set_schedule_hook(&md,&hook_del_atom,&hookparam);
 457         //moldyn_add_schedule(&md,POSTRUN,1.0);
 458
 459         /* activate logging */
 460         moldyn_set_log_dir(&md,argv[1]);
 461         moldyn_set_report(&md,"Frank Zirkelbach",R_TITLE);
 462         moldyn_set_log(&md,LOG_TOTAL_ENERGY,LOG_E);
 463         moldyn_set_log(&md,LOG_TEMPERATURE,LOG_T);
 464         moldyn_set_log(&md,LOG_PRESSURE,LOG_P);
 465         moldyn_set_log(&md,VISUAL_STEP,LOG_V);
 466         moldyn_set_log(&md,SAVE_STEP,LOG_S);
 467         moldyn_set_log(&md,CREATE_REPORT,0);
 468
 469         /* next neighbour distance for critical checking */
 470         set_nn_dist(&md,0.25*ALBE_LC_SI*sqrt(3.0));
 471
 472         /*
 473          * let's do the actual md algorithm now
 474          *
 475          * integration of newtons equations
 476          */
 477         moldyn_integrate(&md);
 478 #ifdef dEBUG
 479 return 0;
 480 #endif
 481
 482         /*
 483          * post processing the data
 484          */
 485
 486         /* close */
 487         moldyn_shutdown(&md);
 488
 489         return 0;
 490 }
 491