]> www.hackdaworld.org Git - physik/posic.git/blob - sic.c
's/for/while/' <- runtime schedule adding possible now
[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 #include "posic.h"
12
13 /* potential */
14 #include "potentials/harmonic_oscillator.h"
15 #include "potentials/lennard_jones.h"
16 #include "potentials/albe.h"
17
18 #ifdef TERSOFF_ORIG
19 #include "potentials/tersoff_orig.h"
20 #else
21 #include "potentials/tersoff.h"
22 #endif
23
24 #define INJECT          1
25 #define NR_ATOMS        4
26
27 int hook(void *moldyn,void *hook_params) {
28
29         t_moldyn *md;
30         t_3dvec r,v,dist;
31         double d;
32         unsigned char run;
33         int i,j;
34         t_atom *atom;
35
36         md=moldyn;
37
38         if(!(md->schedule.count%2)) {
39                 /* add carbon at random place, and enable t scaling */
40                 for(j=0;j<NR_ATOMS;j++) {
41                 run=1;
42                 while(run) {
43                         r.x=rand_get_double(&(md->random))*md->dim.x;
44                         r.y=rand_get_double(&(md->random))*md->dim.y;
45                         r.z=rand_get_double(&(md->random))*md->dim.z;
46                         for(i=0;i<md->count;i++) {
47                                 atom=&(md->atom[i]);
48                                 v3_sub(&dist,&(atom->r),&r);
49                                 d=v3_absolute_square(&dist);
50                                 if(d>TM_R_C)
51                                         run=0;
52                         }
53                 }
54                 v.x=0; v.y=0; v.z=0;
55                 add_atom(md,C,M_C,1,
56                          ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
57                          &r,&v);
58                 }
59                 printf("adding atoms & enable t scaling\n");
60                 set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0);
61         }
62         else {
63                 /* disable t scaling */
64                 printf("disabling t scaling\n");
65                 set_pt_scale(md,0,0,0,0);
66         }
67
68         return 0;
69 }
70
71 int main(int argc,char **argv) {
72
73         /* check argv */
74         if(argc!=3) {
75                 printf("[sic] usage: %s <logdir> <temperatur>\n",argv[0]);
76                 return -1;
77         }
78
79         /* main moldyn structure */
80         t_moldyn md;
81
82         /* potential parameters */
83         t_tersoff_mult_params tp;
84         t_albe_mult_params ap;
85
86         /* atom injection counter */
87         int inject;
88
89         /* testing location & velocity vector */
90         t_3dvec r,v;
91         memset(&r,0,sizeof(t_3dvec));
92         memset(&v,0,sizeof(t_3dvec));
93
94         /* initialize moldyn */
95         moldyn_init(&md,argc,argv);
96
97         /* choose integration algorithm */
98         set_int_alg(&md,MOLDYN_INTEGRATE_VERLET);
99
100         /* choose potential */
101 #ifdef ALBE
102         set_potential3b_j1(&md,albe_mult_3bp_j1);
103         set_potential3b_k1(&md,albe_mult_3bp_k1);
104         set_potential3b_j2(&md,albe_mult_3bp_j2);
105         set_potential3b_k2(&md,albe_mult_3bp_k2);
106 #else
107         set_potential1b(&md,tersoff_mult_1bp);
108         set_potential3b_j1(&md,tersoff_mult_3bp_j1);
109         set_potential3b_k1(&md,tersoff_mult_3bp_k1);
110         set_potential3b_j2(&md,tersoff_mult_3bp_j2);
111         set_potential3b_k2(&md,tersoff_mult_3bp_k2);
112 #endif
113
114 #ifdef ALBE
115         set_potential_params(&md,&ap);
116 #else
117         set_potential_params(&md,&tp);
118 #endif
119
120         /* cutoff radius */
121 #ifdef ALBE
122         set_cutoff(&md,ALBE_S_SI);
123 #else
124         set_cutoff(&md,TM_S_SI);
125 #endif
126
127         /*
128          * potential parameters
129          */
130
131         /*
132          * tersoff mult potential parameters for SiC
133          */
134         tp.S[0]=TM_S_SI;
135         tp.R[0]=TM_R_SI;
136         tp.A[0]=TM_A_SI;
137         tp.B[0]=TM_B_SI;
138         tp.lambda[0]=TM_LAMBDA_SI;
139         tp.mu[0]=TM_MU_SI;
140         tp.beta[0]=TM_BETA_SI;
141         tp.n[0]=TM_N_SI;
142         tp.c[0]=TM_C_SI;
143         tp.d[0]=TM_D_SI;
144         tp.h[0]=TM_H_SI;
145
146         tp.S[1]=TM_S_C;
147         tp.R[1]=TM_R_C;
148         tp.A[1]=TM_A_C;
149         tp.B[1]=TM_B_C;
150         tp.lambda[1]=TM_LAMBDA_C;
151         tp.mu[1]=TM_MU_C;
152         tp.beta[1]=TM_BETA_C;
153         tp.n[1]=TM_N_C;
154         tp.c[1]=TM_C_C;
155         tp.d[1]=TM_D_C;
156         tp.h[1]=TM_H_C;
157
158         tp.chi=TM_CHI_SIC;
159
160         tersoff_mult_complete_params(&tp);
161
162         /*
163          * albe mult potential parameters for SiC
164          */
165         ap.S[0]=ALBE_S_SI;
166         ap.R[0]=ALBE_R_SI;
167         ap.A[0]=ALBE_A_SI;
168         ap.B[0]=ALBE_B_SI;
169         ap.r0[0]=ALBE_R0_SI;
170         ap.lambda[0]=ALBE_LAMBDA_SI;
171         ap.mu[0]=ALBE_MU_SI;
172         ap.gamma[0]=ALBE_GAMMA_SI;
173         ap.c[0]=ALBE_C_SI;
174         ap.d[0]=ALBE_D_SI;
175         ap.h[0]=ALBE_H_SI;
176
177         ap.S[1]=ALBE_S_C;
178         ap.R[1]=ALBE_R_C;
179         ap.A[1]=ALBE_A_C;
180         ap.B[1]=ALBE_B_C;
181         ap.r0[1]=ALBE_R0_C;
182         ap.lambda[1]=ALBE_LAMBDA_C;
183         ap.mu[1]=ALBE_MU_C;
184         ap.gamma[1]=ALBE_GAMMA_C;
185         ap.c[1]=ALBE_C_C;
186         ap.d[1]=ALBE_D_C;
187         ap.h[1]=ALBE_H_C;
188
189         ap.Smixed=ALBE_S_SIC;
190         ap.Rmixed=ALBE_R_SIC;
191         ap.Amixed=ALBE_A_SIC;
192         ap.Bmixed=ALBE_B_SIC;
193         ap.r0_mixed=ALBE_R0_SIC;
194         ap.lambda_m=ALBE_LAMBDA_SIC;
195         ap.mu_m=ALBE_MU_SIC;
196         ap.gamma_m=ALBE_GAMMA_SIC;
197         ap.c_mixed=ALBE_C_SIC;
198         ap.d_mixed=ALBE_D_SIC;
199         ap.h_mixed=ALBE_H_SIC;
200
201         albe_mult_complete_params(&ap);
202
203         /* set (initial) dimensions of simulation volume */
204 #ifdef ALBE
205         //set_dim(&md,8*ALBE_LC_SI,8*ALBE_LC_SI,8*ALBE_LC_SI,TRUE);
206         //set_dim(&md,8*ALBE_LC_C,8*ALBE_LC_C,8*ALBE_LC_C,TRUE);
207         set_dim(&md,8*ALBE_LC_SIC,8*ALBE_LC_SIC,8*ALBE_LC_SIC,TRUE);
208 #else
209         //set_dim(&md,8*LC_SI,8*LC_SI,8*LC_SI,TRUE);
210         //set_dim(&md,8*LC_C,8*LC_C,8*LC_C,TRUE);
211         set_dim(&md,8*TM_LC_SIC,8*TM_LC_SIC,8*TM_LC_SIC,TRUE);
212 #endif
213
214         /* set periodic boundary conditions in all directions */
215         set_pbc(&md,TRUE,TRUE,TRUE);
216
217         /* create the lattice / place atoms */
218 #ifdef ALBE
219         //create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,M_SI,
220         //create_lattice(&md,DIAMOND,ALBE_LC_C,C,M_C,
221 #else
222         //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
223 #endif
224         //               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
225         //               ATOM_ATTR_2BP|ATOM_ATTR_HB,
226         //               0,8,8,8,NULL);
227         //               1,8,8,8,NULL);
228
229         /* create zinkblende structure */
230         /**/
231 #ifdef ALBE
232         r.x=0.5*0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
233         create_lattice(&md,FCC,ALBE_LC_SIC,SI,M_SI,
234                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
235                        0,8,8,8,&r);
236         r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
237         create_lattice(&md,FCC,ALBE_LC_SIC,C,M_C,
238                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
239                        1,8,8,8,&r);
240 #else
241         r.x=0.5*0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
242         create_lattice(&md,FCC,TM_LC_SIC,SI,M_SI,
243                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
244                        0,8,8,8,&r);
245         r.x+=0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
246         create_lattice(&md,FCC,TM_LC_SIC,C,M_C,
247                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
248                        1,8,8,8,&r);
249 #endif
250         /**/
251
252         /* check for right atom placing */
253         moldyn_bc_check(&md);
254
255         /* testing configuration */
256         //r.x=0.27*sqrt(3.0)*LC_SI/2.0; v.x=0;
257         //r.x=(TM_S_SI+TM_R_SI)/4.0;    v.x=0;
258         //r.y=0;                v.y=0;
259         //r.z=0;                v.z=0;
260         //add_atom(&md,SI,M_SI,0,
261         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
262         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
263         //           &r,&v);
264         //r.x=-r.x;     v.x=-v.x;
265         //r.y=0;                v.y=0;
266         //r.z=0;                v.z=0;
267         //add_atom(&md,SI,M_SI,0,
268         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
269         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
270         //           &r,&v);
271         //r.z=0.27*sqrt(3.0)*LC_SI/2.0; v.z=0;
272         //r.x=(TM_S_SI+TM_R_SI)/4.0;    v.x=0;
273         //r.y=0;                v.y=0;
274         //r.x=0;                v.x=0;
275         //add_atom(&md,SI,M_SI,0,
276         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
277         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
278         //           &r,&v);
279         //r.z=-r.z;     v.z=-v.z;
280         //r.y=0;                v.y=0;
281         //r.x=0;                v.x=0;
282         //add_atom(&md,SI,M_SI,0,
283         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
284         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
285         //           &r,&v);
286
287         /* set temperature & pressure */
288         set_temperature(&md,atof(argv[2])+273.0);
289         set_pressure(&md,BAR);
290
291         /* set amount of steps to skip before average calc */
292         set_avg_skip(&md,1000);
293
294         /* set p/t scaling */
295         //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
296         //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,
297         //                 T_SCALE_BERENDSEN,100.0);
298         //set_pt_scale(&md,0,0,T_SCALE_dIRECT,1.0);
299         //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,0,0);
300         
301         /* initial thermal fluctuations of particles (in equilibrium) */
302         thermal_init(&md,TRUE);
303
304         /* create the simulation schedule */
305         /* initial configuration */
306         moldyn_add_schedule(&md,10000,1.0);
307         //moldyn_add_schedule(&md,1000,1.0);
308         //moldyn_add_schedule(&md,1000,1.0);
309         //moldyn_add_schedule(&md,1000,1.0);
310         //moldyn_add_schedule(&md,1000,1.0);
311         //moldyn_add_schedule(&md,1000,1.0);
312         /* adding atoms */
313         //for(inject=0;inject<INJECT;inject++) {
314         //      /* injecting atom and run with enabled t scaling */
315         //      moldyn_add_schedule(&md,900,1.0);
316         //      /* continue running with disabled t scaling */
317         //      moldyn_add_schedule(&md,1100,1.0);
318         //}
319
320
321         /* schedule hook function */
322         moldyn_set_schedule_hook(&md,&hook,NULL);
323
324         /* activate logging */
325         moldyn_set_log_dir(&md,argv[1]);
326         moldyn_set_report(&md,"Frank Zirkelbach","Test 1");
327         moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
328         moldyn_set_log(&md,LOG_TEMPERATURE,1);
329         moldyn_set_log(&md,LOG_PRESSURE,1);
330         moldyn_set_log(&md,VISUAL_STEP,100);
331         moldyn_set_log(&md,SAVE_STEP,100);
332         moldyn_set_log(&md,CREATE_REPORT,0);
333
334         /*
335          * let's do the actual md algorithm now
336          *
337          * integration of newtons equations
338          */
339         moldyn_integrate(&md);
340 #ifdef dEBUG
341 return 0;
342 #endif
343
344         /*
345          * post processing the data
346          */
347
348         /* close */
349         moldyn_shutdown(&md);
350         
351         return 0;
352 }
353