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