]> hackdaworld.org Git - physik/posic.git/blob - sic.c
77e3c79b4a9912808c23269d3800842baae73e90
[physik/posic.git] / 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