]> hackdaworld.org Git - physik/posic.git/blob - posic.c
Makefile ok?
[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         t_atom *si;
21         t_visual vis;
22         t_random random;
23
24         int a,b,c;
25         double e;
26         double help;
27         t_3dvec p;
28         int count;
29
30         t_lj_params lj;
31         t_ho_params ho;
32
33         /* parse arguments */
34         a=moldyn_parse_argv(&md,argc,argv);
35         if(a<0) return -1;
36
37         /* init */
38         moldyn_log_init(&md,&vis);
39         rand_init(&random,NULL,1);
40         random.status|=RAND_STAT_VERBOSE;
41
42         /* testing random numbers */
43         //for(a=0;a<1000000;a++)
44         //      printf("%f %f\n",rand_get_gauss(&random),
45         //                       rand_get_gauss(&random));
46
47         a=LEN_X;
48         b=LEN_Y;
49         c=LEN_Z;
50
51         /* set for 'bounding atoms' */
52         vis.dim.x=a*LC_SI;
53         vis.dim.y=b*LC_SI;
54         vis.dim.z=c*LC_SI;
55
56         /* init lattice */
57         printf("placing silicon atoms ... ");
58         count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si);
59         printf("(%d) ok!\n",count);
60         /* testing purpose
61         count=2;
62         si=malloc(2*sizeof(t_atom));
63         si[0].r.x=0.13*sqrt(3.0)*LC_SI/2.0;
64         si[0].r.y=0;
65         si[0].r.z=0;
66         si[0].element=SI;
67         si[0].mass=M_SI;
68         si[1].r.x=-si[0].r.x;
69         si[1].r.y=0;
70         si[1].r.z=0;
71         si[1].element=SI;
72         si[1].mass=M_SI;
73         */
74
75         /* moldyn init (now si is a valid address) */
76         md.count=count;
77         md.atom=si;
78         md.potential=potential_lennard_jones;
79         md.force=force_lennard_jones;
80         //md.potential=potential_harmonic_oscillator;
81         //md.force=force_harmonic_oscillator;
82         md.cutoff=R_CUTOFF*LC_SI;
83         md.cutoff_square=md.cutoff*md.cutoff;
84         md.pot_params=&lj;
85         //md.pot_params=&ho;
86         md.integrate=velocity_verlet;
87         //md.time_steps=RUNS;
88         //md.tau=TAU;
89         md.status=0;
90         md.visual=&vis;
91         /* dimensions of the simulation cell */
92         md.dim.x=a*LC_SI;
93         md.dim.y=b*LC_SI;
94         md.dim.z=c*LC_SI;
95
96         /* verlet list init */
97         // later integrated in moldyn_init function!
98         verlet_list_init(&md);
99
100         printf("setting thermal fluctuations (T=%f K)\n",md.t);
101         thermal_init(&md,&random,count);
102         //for(a=0;a<count;a++) v3_zero(&(si[0].v));
103
104         /* check kinetic energy */
105
106         e=get_e_kin(si,count);
107         printf("kinetic energy: %.40f [J]\n",e);
108         printf("3/2 N k T = %.40f [J] (T=%f [K])\n",
109                1.5*count*K_BOLTZMANN*md.t,md.t);
110
111         /* check total momentum */
112         p=get_total_p(si,count);
113         printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
114
115         /* potential paramters */
116         lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
117         help=lj.sigma6*lj.sigma6;
118         lj.sigma6*=help;
119         lj.sigma12=lj.sigma6*lj.sigma6;
120         lj.epsilon4=4.0*LJ_EPSILON_SI;
121
122         ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
123         ho.spring_constant=1.0;
124
125         printf("estimated accurate time step: %.30f [s]\n",
126                estimate_time_step(&md,3.0,md.t));
127
128
129         /*
130          * let's do the actual md algorithm now
131          *
132          * integration of newtons equations
133          */
134
135         moldyn_integrate(&md);
136
137         printf("total energy (after integration): %.40f [J]\n",
138                get_total_energy(&md));
139
140         /* close */
141
142         verlet_list_shutdown(&md);
143
144         rand_close(&random);
145
146         moldyn_shutdown(&md);
147         
148         return 0;
149 }
150