2 * sic.c - investigation of the sic precipitation process of silicon carbide
4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
13 #include "potentials/harmonic_oscillator.h"
14 #include "potentials/lennard_jones.h"
15 #include "potentials/albe.h"
17 #include "potentials/tersoff_orig.h"
19 #include "potentials/tersoff.h"
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 */
31 #define STATE_PRERUN 0x00
32 #define STATE_INSERT 0x01
33 #define STATE_POSTRUN 0x02
35 /* include the config file */
38 int insert_atoms(t_moldyn *moldyn) {
51 for(j=0;j<INS_ATOMS;j++) {
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;
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;
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;
84 /* assume valid coordinates */
86 for(i=0;i<moldyn->count;i++) {
87 atom=&(moldyn->atom[i]);
88 v3_sub(&dist,&(atom->r),&r);
89 check_per_bound(moldyn,&dist);
90 d=v3_absolute_square(&dist);
91 /* reject coordinates */
93 //printf("atom %d - %f\n",i,d);
99 add_atom(moldyn,INS_TYPE,INS_MASS,INS_BRAND,
100 ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|\
101 //ATOM_ATTR_HB|ATOM_ATTR_VB,
109 int sic_hook(void *moldyn,void *hook_params) {
119 /* switch on t scaling */
120 if(md->schedule.count==0)
121 set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0);
123 /* my lousy state machine ! */
125 /* switch to insert state immediately */
126 if(hp->state==STATE_PRERUN)
127 hp->state=STATE_INSERT;
129 /* act according to state */
132 /* check temperature */
133 if(md->t_avg-md->t_ref>INS_DELTA_TC) {
140 printf(" ### insert atoms (%d/%d) ###\n",
141 hp->insert_count*INS_ATOMS,INS_RUNS*INS_ATOMS);
143 /* change state after last insertion */
144 if(hp->insert_count==INS_RUNS)
145 hp->state=STATE_POSTRUN;
149 if(md->t-md->t_ref>POST_DELTA_TC) {
153 /* decrease temperature */
154 hp->postrun_count+=1;
155 printf(" ### postrun (%d/%d) ###\n",
156 hp->postrun_count,POST_RUNS);
157 set_temperature(md,md->t_ref-POST_DT);
158 if(hp->postrun_count==POST_RUNS)
162 printf("[hook] FATAL (default case!?!)\n");
167 moldyn_add_schedule(md,steps,tau);
172 int main(int argc,char **argv) {
174 /* main moldyn structure */
177 /* hook parameter structure */
180 /* potential parameters */
181 t_tersoff_mult_params tp;
182 t_albe_mult_params ap;
184 /* testing location & velocity vector */
186 memset(&r,0,sizeof(t_3dvec));
187 memset(&v,0,sizeof(t_3dvec));
189 /* initialize moldyn */
190 moldyn_init(&md,argc,argv);
192 /* choose integration algorithm */
193 set_int_alg(&md,MOLDYN_INTEGRATE_VERLET);
195 /* choose potential */
197 set_potential3b_j1(&md,albe_mult_3bp_j1);
198 set_potential3b_k1(&md,albe_mult_3bp_k1);
199 set_potential3b_j2(&md,albe_mult_3bp_j2);
200 set_potential3b_k2(&md,albe_mult_3bp_k2);
202 set_potential1b(&md,tersoff_mult_1bp);
203 set_potential3b_j1(&md,tersoff_mult_3bp_j1);
204 set_potential3b_k1(&md,tersoff_mult_3bp_k1);
205 set_potential3b_j2(&md,tersoff_mult_3bp_j2);
206 set_potential3b_k2(&md,tersoff_mult_3bp_k2);
210 set_potential_params(&md,&ap);
212 set_potential_params(&md,&tp);
215 /* cutoff radius & bondlen */
217 set_cutoff(&md,ALBE_S_SI);
218 set_bondlen(&md,ALBE_S_SI,ALBE_S_C,ALBE_S_SIC);
219 //set_cutoff(&md,ALBE_S_C);
221 set_cutoff(&md,TM_S_SI);
222 set_bondlen(&md,TM_S_SI,TM_S_C,-1.0);
223 //set_cutoff(&md,TM_S_C);
227 * potential parameters
231 * tersoff mult potential parameters for SiC
237 tp.lambda[0]=TM_LAMBDA_SI;
239 tp.beta[0]=TM_BETA_SI;
249 tp.lambda[1]=TM_LAMBDA_C;
251 tp.beta[1]=TM_BETA_C;
259 tersoff_mult_complete_params(&tp);
262 * albe mult potential parameters for SiC
269 ap.lambda[0]=ALBE_LAMBDA_SI;
271 ap.gamma[0]=ALBE_GAMMA_SI;
281 ap.lambda[1]=ALBE_LAMBDA_C;
283 ap.gamma[1]=ALBE_GAMMA_C;
288 ap.Smixed=ALBE_S_SIC;
289 ap.Rmixed=ALBE_R_SIC;
290 ap.Amixed=ALBE_A_SIC;
291 ap.Bmixed=ALBE_B_SIC;
292 ap.r0_mixed=ALBE_R0_SIC;
293 ap.lambda_m=ALBE_LAMBDA_SIC;
295 ap.gamma_m=ALBE_GAMMA_SIC;
296 ap.c_mixed=ALBE_C_SIC;
297 ap.d_mixed=ALBE_D_SIC;
298 ap.h_mixed=ALBE_H_SIC;
300 albe_mult_complete_params(&ap);
302 /* set (initial) dimensions of simulation volume */
304 set_dim(&md,LCNTX*ALBE_LC_SI,LCNTY*ALBE_LC_SI,LCNTZ*ALBE_LC_SI,TRUE);
305 //set_dim(&md,LCNTX*ALBE_LC_C,LCNTY*ALBE_LC_C,LCNTZ*ALBE_LC_C,TRUE);
306 //set_dim(&md,LCNTX*ALBE_LC_SIC,LCNTY*ALBE_LC_SIC,LCNTZ*ALBE_LC_SIC,TRUE);
308 set_dim(&md,LCNTX*LC_SI,LCNTY*LC_SI,LCNTZ*LC_SI,TRUE);
309 //set_dim(&md,LCNTX*LC_C,LCNTY*LC_C,LCNTZ*LC_C,TRUE);
310 //set_dim(&md,LCNTX*TM_LC_SIC,LCNTY*TM_LC_SIC,LCNTZ*TM_LC_SIC,TRUE);
313 /* set periodic boundary conditions in all directions */
314 set_pbc(&md,TRUE,TRUE,TRUE);
316 /* create the lattice / place atoms */
319 create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,M_SI,
320 //create_lattice(&md,DIAMOND,ALBE_LC_C,C,M_C,
322 create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
324 ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
325 // ATOM_ATTR_2BP|ATOM_ATTR_HB,
326 0,LCNTX,LCNTY,LCNTZ,NULL);
327 // 1,LCNTX,LCNTY,LCNTZ,NULL);
329 /* create zinkblende structure */
332 r.x=0.5*0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
333 create_lattice(&md,FCC,ALBE_LC_SIC,SI,M_SI,
334 ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
335 0,LCNTX,LCNTY,LCNTZ,&r);
336 r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
337 create_lattice(&md,FCC,ALBE_LC_SIC,C,M_C,
338 ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
339 1,LCNTX,LCNTY,LCNTZ,&r);
341 r.x=0.5*0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
342 create_lattice(&md,FCC,TM_LC_SIC,SI,M_SI,
343 ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
344 0,LCNTX,LCNTY,LCNTZ,&r);
345 r.x+=0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
346 create_lattice(&md,FCC,TM_LC_SIC,C,M_C,
347 ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
348 1,LCNTX,LCNTY,LCNTZ,&r);
352 /* check for right atom placing */
353 moldyn_bc_check(&md);
355 /* testing configuration */
356 //r.x=0.27*sqrt(3.0)*LC_SI/2.0; v.x=0;
357 //r.x=(TM_S_SI+TM_R_SI)/4.0; v.x=0;
360 //add_atom(&md,SI,M_SI,0,
361 // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
362 // ATOM_ATTR_2BP|ATOM_ATTR_HB,
364 //r.x=-r.x; v.x=-v.x;
367 //add_atom(&md,SI,M_SI,0,
368 // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
369 // ATOM_ATTR_2BP|ATOM_ATTR_HB,
371 //r.z=0.27*sqrt(3.0)*LC_SI/2.0; v.z=0;
372 //r.x=(TM_S_SI+TM_R_SI)/4.0; v.x=0;
375 //add_atom(&md,SI,M_SI,0,
376 // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
377 // ATOM_ATTR_2BP|ATOM_ATTR_HB,
379 //r.z=-r.z; v.z=-v.z;
382 //add_atom(&md,SI,M_SI,0,
383 // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
384 // ATOM_ATTR_2BP|ATOM_ATTR_HB,
387 /* set temperature & pressure */
388 set_temperature(&md,atof(argv[2])+273.0);
389 set_pressure(&md,BAR);
391 /* set amount of steps to skip before average calc */
392 set_avg_skip(&md,AVG_SKIP);
394 /* set p/t scaling */
395 //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
396 //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,
397 // T_SCALE_BERENDSEN,100.0);
398 //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
399 //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,0,0);
401 /* initial thermal fluctuations of particles (in equilibrium) */
402 thermal_init(&md,TRUE);
404 /* create the simulation schedule */
405 moldyn_add_schedule(&md,PRERUN,PRE_TAU);
407 /* schedule hook function */
408 memset(&hookparam,0,sizeof(t_hp));
411 moldyn_set_schedule_hook(&md,&sic_hook,&hookparam);
412 //moldyn_set_schedule_hook(&md,&hook_del_atom,&hookparam);
413 //moldyn_add_schedule(&md,POSTRUN,1.0);
415 /* activate logging */
416 moldyn_set_log_dir(&md,argv[1]);
417 moldyn_set_report(&md,"Frank Zirkelbach",R_TITLE);
418 moldyn_set_log(&md,LOG_TOTAL_ENERGY,LOG_E);
419 moldyn_set_log(&md,LOG_TEMPERATURE,LOG_T);
420 moldyn_set_log(&md,LOG_PRESSURE,LOG_P);
421 moldyn_set_log(&md,VISUAL_STEP,LOG_V);
422 moldyn_set_log(&md,SAVE_STEP,LOG_S);
423 moldyn_set_log(&md,CREATE_REPORT,0);
426 * let's do the actual md algorithm now
428 * integration of newtons equations
430 moldyn_integrate(&md);
436 * post processing the data
440 moldyn_shutdown(&md);