]> hackdaworld.org Git - physik/posic.git/blob - posic.c
added post proc stuff
[physik/posic.git] / posic.c
1 /*
2  * posic.c - precipitation process of silicon carbide in silicon
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 "math/math.h"
12 #include "init/init.h"
13 #include "visual/visual.h"
14
15 #include "posic.h"
16
17 int main(int argc,char **argv) {
18
19         t_moldyn md;
20
21         t_lj_params lj;
22         t_ho_params ho;
23         t_tersoff_mult_params tp;
24
25         int a,b,c;
26         double e;
27         double help;
28         t_3dvec p;
29
30         /*
31          *  moldyn init
32          *
33          * - parsing argv
34          * - log init
35          * - random init
36          *
37          */
38         a=moldyn_init(&md,argc,argv);
39         if(a<0) return a;
40
41         /*
42          * the following overrides possibly set interaction methods by argv !!
43          */
44
45         /* params */
46         lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
47         help=lj.sigma6*lj.sigma6;
48         lj.sigma6*=help;
49         lj.sigma12=lj.sigma6*lj.sigma6;
50         lj.epsilon4=4.0*LJ_EPSILON_SI;
51         ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
52         ho.spring_constant=1;
53         /* assignement */
54         md.potential_force_function=lennard_jones;
55         //md.potential_force_function=harmonic_oscillator;
56         md.pot_params=&lj;
57         //md.pot_params=&ho;
58         /* cutoff radius */
59         md.cutoff=R_CUTOFF*LC_SI;
60
61         /*
62          * testing random numbers
63          */
64
65 #ifdef DEBUG_RANDOM_NUMBER
66         for(a=0;a<1000000;a++)
67                 printf("%f %f\n",rand_get_gauss(&(md.random)),
68                                  rand_get_gauss(&(md.random)));
69         return 0;
70 #endif
71
72         /*
73          * geometry & particles
74          */
75
76         /* simulation cell volume in lattice constants */
77         a=LEN_X;
78         b=LEN_Y;
79         c=LEN_Z;
80         md.dim.x=a*LC_SI;
81         md.dim.y=b*LC_SI;
82         md.dim.z=c*LC_SI;
83
84         /* (un)set to (not) get visualized 'bounding atoms' */
85         md.vis.dim.x=a*LC_SI;
86         md.vis.dim.y=b*LC_SI;
87         md.vis.dim.z=c*LC_SI;
88
89         /*
90          * particles
91          */
92
93         /* lattice init */
94
95 #ifndef SIMPLE_TESTING
96         md.count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&(md.atom));
97         printf("created silicon lattice (#atoms = %d)\n",md.count);
98 #else
99         md.count=2;
100         md.atom=malloc(md.count*sizeof(t_atom));
101         md.atom[0].r.x=0.23*sqrt(3.0)*LC_SI/2.0;
102         md.atom[0].r.y=0;
103         md.atom[0].r.z=0;
104         md.atom[0].element=SI;
105         md.atom[0].mass=M_SI;
106         md.atom[1].r.x=-md.atom[0].r.x;
107         md.atom[1].r.y=0;
108         md.atom[1].r.z=0;
109         md.atom[1].element=SI;
110         md.atom[1].mass=M_SI;
111
112         //md.atom[2].r.x=0.5*(a-1)*LC_SI;
113         //md.atom[2].r.y=0.5*(b-1)*LC_SI;
114         //md.atom[2].r.z=0;
115         //md.atom[2].element=C;
116         //md.atom[2].mass=M_C;
117
118         //md.atom[3].r.x=0.5*(a-1)*LC_SI;
119         //md.atom[3].r.y=0;
120         //md.atom[3].r.z=0;
121         //md.atom[3].element=SI;
122         //md.atom[3].mass=M_SI;
123 #endif
124
125         /* initial thermal fluctuations of particles */
126
127 #ifndef SIMPLE_TESTING
128         printf("setting thermal fluctuations (T=%f K)\n",md.t);
129         thermal_init(&md);
130 #else
131         for(a=0;a<md.count;a++) v3_zero(&(md.atom[0].v));
132         md.atom[2].v.x=-320;
133         md.atom[2].v.y=-320;
134 #endif
135
136         /* check kinetic energy */
137         e=get_e_kin(md.atom,md.count);
138         printf("kinetic energy: %.40f [J]\n",e);
139         printf("3/2 N k T = %.40f [J] (T=%f [K])\n",
140                1.5*md.count*K_BOLTZMANN*md.t,md.t);
141
142         /* check total momentum */
143         p=get_total_p(md.atom,md.count);
144         printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
145
146         /* check time step */
147         printf("estimated accurate time step: %.30f [s]\n",
148                estimate_time_step(&md,3.0,md.t));
149
150         /*
151          * let's do the actual md algorithm now
152          *
153          * integration of newtons equations
154          */
155
156         moldyn_integrate(&md);
157
158         printf("total energy (after integration): %.40f [J]\n",
159                get_total_energy(&md));
160
161         /* close */
162
163         link_cell_shutdown(&md);
164
165         moldyn_shutdown(&md);
166         
167         return 0;
168 }
169