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