]> hackdaworld.org Git - physik/posic.git/blob - sic.c
added Posselt formation energy calc
[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                         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.1)*ALBE_LC_SI;
79                         moldyn->atom[4372].r.z=(-0.4)*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 #ifdef INS_DYNAMIC_LEN
90                         r.x=(rand_get_double(&(moldyn->random))-0.5)*\
91                             moldyn->dim.x;
92                         r.y=(rand_get_double(&(moldyn->random))-0.5)*\
93                             moldyn->dim.y;
94                         r.z=(rand_get_double(&(moldyn->random))-0.5)*\
95                             moldyn->dim.z;
96 #else
97                         r.x=(rand_get_double(&(moldyn->random))-0.5)*INS_LENX;
98                         r.y=(rand_get_double(&(moldyn->random))-0.5)*INS_LENY;
99                         r.z=(rand_get_double(&(moldyn->random))-0.5)*INS_LENZ;
100 #endif
101 #endif
102                         // offset
103                         r.x+=INS_OFFSET;
104                         r.y+=INS_OFFSET;
105                         r.z+=INS_OFFSET;
106                         /* assume valid coordinates */
107                         run=0;
108                         dmin=10000000000.0;             // for sure too high!
109                         for(i=0;i<moldyn->count;i++) {
110                                 atom=&(moldyn->atom[i]);
111                                 v3_sub(&dist,&(atom->r),&r);
112                                 check_per_bound(moldyn,&dist);
113                                 d=v3_absolute_square(&dist);
114                                 /* reject coordinates */
115                                 if(d<INS_R_C) {
116                                         //printf("atom %d - %f\n",i,d);
117                                         run=1;
118                                         break;
119                                 }
120                                 if(d<dmin)
121                                         dmin=d;
122                         }
123                 }
124                 add_atom(moldyn,INS_TYPE,INS_BRAND,
125                          ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|\
126                          INS_ATTR,
127                          &r,&v);
128                 printf(" %02d: atom %d | %f %f %f | %f\n",
129                        j,moldyn->count-1,r.x,r.y,r.z,dmin);
130         }
131
132         return 0;
133 }
134
135 int sic_hook(void *moldyn,void *hook_params) {
136
137         t_hp *hp;
138         t_moldyn *md;
139         int steps;
140         double tau;
141         double dt;
142         double dp;
143
144         hp=hook_params;
145         md=moldyn;
146
147         tau=1.0;
148         steps=0;
149
150         /* switch on t scaling */
151         if(md->schedule.count==0)
152                 set_pt_scale(md,P_SCALE_BERENDSEN,P_SCALE_TAU,
153                                 T_SCALE_BERENDSEN,T_SCALE_TAU);
154
155         /* my lousy state machine ! */
156
157         /* switch to insert state immediately */
158         if(hp->state==STATE_PRERUN)
159                 hp->state=STATE_INSERT;
160         
161         switch(hp->state) {
162                 case STATE_INSERT:
163                         goto insert;
164                         break;
165                 case STATE_POSTRUN:
166                         goto postrun;
167                         break;
168                 default:
169                         printf("[sic hook] unknown state\n");
170                         return -1;
171         }
172
173         /* act according to state */
174
175 insert:
176
177         /* assigne values */
178         steps=INS_RELAX;
179         tau=INS_TAU;
180
181         /* check temperature */
182         dt=md->t_avg-md->t_ref;
183         dp=md->p_avg-md->p_ref;
184         if(dt<0)
185                 dt=-dt;
186         if(dp<0)
187                 dp=-dp;
188         if((dt>INS_DELTA_TC)|(dp>INS_DELTA_PC))
189                 goto addsched;
190
191         /* immediately go on if no job is to be done */
192         if(hp->insert_count==INS_RUNS) {
193                 printf("    --- leaving insert state ---\n");
194                 hp->state=STATE_POSTRUN;
195                 goto postrun;
196         }
197
198         /* else -> insert atoms */
199         hp->insert_count+=1;
200         printf("   ### insert atoms (%d/%d) ###\n",
201                hp->insert_count*INS_ATOMS,INS_RUNS*INS_ATOMS);
202         insert_atoms(md);
203         goto addsched;
204
205 postrun:
206
207         /* assigne values */
208         steps=POST_RELAX;
209         tau=POST_TAU;
210
211         /* check temperature */
212         dt=md->t_avg-md->t_ref;
213         dp=md->p_avg-md->p_ref;
214         if(dt<0)
215                 dt=-dt;
216         if(dp<0)
217                 dp=-dp;
218         if((dt>POST_DELTA_TC)|(dp>POST_DELTA_PC))
219                 goto addsched;
220
221         /* immediately return if no job is to be done */
222         if(hp->postrun_count==POST_RUNS) {
223                 printf("    --- leaving post run state ---\n");
224                 return 0;
225         }
226
227         /* postrun action */
228         hp->postrun_count+=1;
229         printf(" ### postrun (%d/%d) ###\n",
230                hp->postrun_count,POST_RUNS);
231         set_temperature(md,md->t_ref-POST_DT);
232
233 addsched:
234
235         /* reset the average counters */
236         average_reset(md);
237
238         /* add schedule */
239         moldyn_add_schedule(md,steps,tau);
240
241         return 0;
242 }
243
244 int main(int argc,char **argv) {
245
246         /* main moldyn structure */
247         t_moldyn md;
248
249         /* hook parameter structure */
250         t_hp hookparam;
251
252         /* testing location & velocity vector */
253         t_3dvec r,v;
254         memset(&r,0,sizeof(t_3dvec));
255         memset(&v,0,sizeof(t_3dvec));
256
257         /* initialize moldyn */
258         moldyn_init(&md,argc,argv);
259
260         /* choose integration algorithm */
261         set_int_alg(&md,MOLDYN_INTEGRATE_VERLET);
262
263         /* choose potential */
264 #ifdef ALBE
265         if(set_potential(&md,MOLDYN_POTENTIAL_AM)<0)
266                 return -1;
267 #else
268         if(set_potential(&md,MOLDYN_POTENTIAL_TM)<0)
269                 return -1;
270 #endif
271
272         /* cutoff radius & bondlen */
273 #ifdef ALBE
274         set_cutoff(&md,ALBE_S_SI);
275         //set_cutoff(&md,ALBE_S_C);
276 #else
277         set_cutoff(&md,TM_S_SI);
278         //set_cutoff(&md,TM_S_C);
279 #endif
280
281         /*
282          * potential parameters
283          */
284
285 #ifndef ALBE
286         /*
287          * tersoff mult potential parameters for SiC
288          */
289         tersoff_mult_set_params(&md,SI,C);
290 #else
291         /*
292          * albe mult potential parameters for SiC
293          */
294         albe_mult_set_params(&md,SI,C);
295 #endif
296
297         /* set (initial) dimensions of simulation volume */
298 #ifdef ALBE
299  #ifdef INIT_SI
300         set_dim(&md,LCNTX*ALBE_LC_SI,LCNTY*ALBE_LC_SI,LCNTZ*ALBE_LC_SI,TRUE);
301  #endif
302  #ifdef INIT_C
303         set_dim(&md,LCNTX*ALBE_LC_C,LCNTY*ALBE_LC_C,LCNTZ*ALBE_LC_C,TRUE);
304  #endif
305  #ifdef INIT_3CSIC
306         set_dim(&md,LCNTX*ALBE_LC_SIC,LCNTY*ALBE_LC_SIC,LCNTZ*ALBE_LC_SIC,TRUE);
307  #endif
308 #else
309  #ifdef INIT_SI
310         set_dim(&md,LCNTX*LC_SI,LCNTY*LC_SI,LCNTZ*LC_SI,TRUE);
311  #endif
312  #ifdef INIT_C
313         set_dim(&md,LCNTX*LC_C,LCNTY*LC_C,LCNTZ*LC_C,TRUE);
314  #endif
315  #ifdef INIT_3CSIC
316         set_dim(&md,LCNTX*TM_LC_SIC,LCNTY*TM_LC_SIC,LCNTZ*TM_LC_SIC,TRUE);
317  #endif
318 #endif
319
320         /* set periodic boundary conditions in all directions */
321         set_pbc(&md,TRUE,TRUE,TRUE);
322
323         /* create the lattice / place atoms */
324
325         // diamond
326 #ifdef ALBE
327  #ifdef INIT_SI
328         create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,
329                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
330                        0,LCNTX,LCNTY,LCNTZ,NULL,0,NULL);
331  #endif
332  #ifdef INIT_C
333         create_lattice(&md,DIAMOND,ALBE_LC_C,C,
334                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
335                        1,LCNTX,LCNTY,LCNTZ,NULL,0,NULL);
336  #endif
337 #else
338  #ifdef INIT_SI
339         create_lattice(&md,DIAMOND,LC_SI,SI,
340                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
341                        0,LCNTX,LCNTY,LCNTZ,NULL,0,NULL);
342  #endif
343  #ifdef INIT_C
344         create_lattice(&md,DIAMOND,LC_C,SI,
345                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
346                        1,LCNTX,LCNTY,LCNTZ,NULL,0,NULL);
347  #endif
348 #endif
349
350         // zinkblende 
351 #ifdef INIT_3CSIC
352  #ifdef ALBE
353         r.x=0.5*0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
354         create_lattice(&md,FCC,ALBE_LC_SIC,SI,
355                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
356                        0,LCNTX,LCNTY,LCNTZ,&r,0,NULL);
357         r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
358         create_lattice(&md,FCC,ALBE_LC_SIC,C,
359                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB|ATOM_ATTR_VB,
360                        1,LCNTX,LCNTY,LCNTZ,&r,0,NULL);
361  #else
362         r.x=0.5*0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
363         create_lattice(&md,FCC,TM_LC_SIC,SI,
364                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
365                        0,LCNTX,LCNTY,LCNTZ,&r,0,NULL);
366         r.x+=0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
367         create_lattice(&md,FCC,TM_LC_SIC,C,
368                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
369                        1,LCNTX,LCNTY,LCNTZ,&r,0,NULL);
370  #endif
371 #endif
372
373         /* check for right atom placing */
374         moldyn_bc_check(&md);
375
376         /* testing configuration */
377         //r.x=0.27*sqrt(3.0)*LC_SI/2.0; v.x=0;
378         //r.x=(TM_S_SI+TM_R_SI)/4.0;    v.x=0;
379         //r.y=0;                v.y=0;
380         //r.z=0;                v.z=0;
381         //add_atom(&md,SI,M_SI,0,
382         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
383         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
384         //           &r,&v);
385         //r.x=-r.x;     v.x=-v.x;
386         //r.y=0;                v.y=0;
387         //r.z=0;                v.z=0;
388         //add_atom(&md,SI,M_SI,0,
389         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
390         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
391         //           &r,&v);
392         //r.z=0.27*sqrt(3.0)*LC_SI/2.0; v.z=0;
393         //r.x=(TM_S_SI+TM_R_SI)/4.0;    v.x=0;
394         //r.y=0;                v.y=0;
395         //r.x=0;                v.x=0;
396         //add_atom(&md,SI,M_SI,0,
397         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
398         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
399         //           &r,&v);
400         //r.z=-r.z;     v.z=-v.z;
401         //r.y=0;                v.y=0;
402         //r.x=0;                v.x=0;
403         //add_atom(&md,SI,M_SI,0,
404         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
405         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
406         //           &r,&v);
407
408         /* set temperature & pressure */
409         set_temperature(&md,atof(argv[2])+273.0);
410         set_pressure(&md,0.0);
411
412         /* set amount of steps to skip before average calc */
413         set_avg_skip(&md,AVG_SKIP);
414
415         /* set p/t scaling */
416         //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
417         //set_pt_scale(&md,P_SCALE_BERENDSEN,0.01/(100*GPA),
418         //                 T_SCALE_BERENDSEN,100.0);
419         //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
420         //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,0,0);
421         
422         /* initial thermal fluctuations of particles (in equilibrium) */
423         thermal_init(&md,TRUE);
424
425         /* create the simulation schedule */
426         moldyn_add_schedule(&md,PRERUN,PRE_TAU);
427
428         /* schedule hook function */
429         memset(&hookparam,0,sizeof(t_hp));
430         hookparam.argc=argc;
431         hookparam.argv=argv;
432         moldyn_set_schedule_hook(&md,&sic_hook,&hookparam);
433         //moldyn_set_schedule_hook(&md,&hook_del_atom,&hookparam);
434         //moldyn_add_schedule(&md,POSTRUN,1.0);
435
436         /* activate logging */
437         moldyn_set_log_dir(&md,argv[1]);
438         moldyn_set_report(&md,"Frank Zirkelbach",R_TITLE);
439         moldyn_set_log(&md,LOG_TOTAL_ENERGY,LOG_E);
440         moldyn_set_log(&md,LOG_TEMPERATURE,LOG_T);
441         moldyn_set_log(&md,LOG_PRESSURE,LOG_P);
442         moldyn_set_log(&md,LOG_VOLUME,LOG_V);
443         moldyn_set_log(&md,VISUAL_STEP,LOG_A);
444         moldyn_set_log(&md,SAVE_STEP,LOG_S);
445         moldyn_set_log(&md,CREATE_REPORT,0);
446
447         /* next neighbour distance for critical checking */
448         set_nn_dist(&md,0.25*ALBE_LC_SI*sqrt(3.0));
449
450         /*
451          * let's do the actual md algorithm now
452          *
453          * integration of newtons equations
454          */
455         moldyn_integrate(&md);
456 #ifdef dEBUG
457 return 0;
458 #endif
459
460         /*
461          * post processing the data
462          */
463
464         /* close */
465         moldyn_shutdown(&md);
466         
467         return 0;
468 }
469