2 * posic.c - precipitation process of silicon carbide in silicon
4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
11 #include "math/math.h"
12 #include "init/init.h"
13 #include "visual/visual.h"
17 int main(int argc,char **argv) {
33 a=moldyn_parse_argv(&md,argc,argv);
37 moldyn_log_init(&md,&vis);
38 rand_init(&random,NULL,1);
39 random.status|=RAND_STAT_VERBOSE;
41 /* testing random numbers */
42 //for(a=0;a<1000000;a++)
43 // printf("%f %f\n",rand_get_gauss(&random),
44 // rand_get_gauss(&random));
50 /* set for 'bounding atoms' */
56 printf("placing silicon atoms ... ");
57 md.count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si);
58 printf("(%d) ok!\n",md.count);
61 si=malloc(2*sizeof(t_atom));
62 si[0].r.x=0.13*sqrt(3.0)*LC_SI/2.0;
74 /* moldyn init (now si is a valid address) */
76 md.potential_force_function=lennard_jones;
77 //md.potential_force_function=harmonic_oscillator;
78 md.cutoff=R_CUTOFF*LC_SI;
83 /* dimensions of the simulation cell */
88 printf("setting thermal fluctuations (T=%f K)\n",md.t);
89 // thermal_init(&md,&random);
90 for(a=0;a<md.count;a++) v3_zero(&(si[0].v));
92 /* check kinetic energy */
94 e=get_e_kin(si,md.count);
95 printf("kinetic energy: %.40f [J]\n",e);
96 printf("3/2 N k T = %.40f [J] (T=%f [K])\n",
97 1.5*md.count*K_BOLTZMANN*md.t,md.t);
99 /* check total momentum */
100 p=get_total_p(si,md.count);
101 printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
103 /* potential paramters */
104 lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
105 help=lj.sigma6*lj.sigma6;
107 lj.sigma12=lj.sigma6*lj.sigma6;
108 lj.epsilon4=4.0*LJ_EPSILON_SI;
110 ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
111 ho.spring_constant=1.0;
113 printf("estimated accurate time step: %.30f [s]\n",
114 estimate_time_step(&md,3.0,md.t));
117 * let's do the actual md algorithm now
119 * integration of newtons equations
122 moldyn_integrate(&md);
124 printf("total energy (after integration): %.40f [J]\n",
125 get_total_energy(&md));
129 link_cell_shutdown(&md);
133 moldyn_shutdown(&md);