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) {
37 a=moldyn_init(&md,argc,argv);
41 * the following overrides possibly set interaction methods by argv !!
45 lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
46 help=lj.sigma6*lj.sigma6;
48 lj.sigma12=lj.sigma6*lj.sigma6;
49 lj.epsilon4=4.0*LJ_EPSILON_SI;
50 ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
53 md.potential_force_function=lennard_jones;
54 //md.potential_force_function=harmonic_oscillator;
58 md.cutoff=R_CUTOFF*LC_SI;
61 * testing random numbers
64 #ifdef DEBUG_RANDOM_NUMBER
65 for(a=0;a<1000000;a++)
66 printf("%f %f\n",rand_get_gauss(&(md.random)),
67 rand_get_gauss(&(md.random)));
72 * geometry & particles
75 /* simulation cell volume in lattice constants */
83 /* (un)set to (not) get visualized 'bounding atoms' */
94 #ifndef SIMPLE_TESTING
95 md.count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&(md.atom));
96 printf("created silicon lattice (#atoms = %d)\n",md.count);
99 md.atom=malloc(md.count*sizeof(t_atom));
100 md.atom[0].r.x=0.21*sqrt(3.0)*LC_SI/2.0;
103 md.atom[0].element=SI;
104 md.atom[0].mass=M_SI;
105 md.atom[1].r.x=-md.atom[0].r.x;
108 md.atom[1].element=SI;
109 md.atom[1].mass=M_SI;
111 md.atom[2].r.x=0.5*(a-1)*LC_SI;
112 md.atom[2].r.y=0.5*(b-1)*LC_SI;
114 md.atom[2].element=C;
117 //md.atom[3].r.x=0.5*(a-1)*LC_SI;
120 //md.atom[3].element=SI;
121 //md.atom[3].mass=M_SI;
124 /* initial thermal fluctuations of particles */
126 #ifndef SIMPLE_TESTING
127 printf("setting thermal fluctuations (T=%f K)\n",md.t);
130 for(a=0;a<md.count;a++) v3_zero(&(md.atom[0].v));
135 /* check kinetic energy */
136 e=get_e_kin(md.atom,md.count);
137 printf("kinetic energy: %.40f [J]\n",e);
138 printf("3/2 N k T = %.40f [J] (T=%f [K])\n",
139 1.5*md.count*K_BOLTZMANN*md.t,md.t);
141 /* check total momentum */
142 p=get_total_p(md.atom,md.count);
143 printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
145 /* check time step */
146 printf("estimated accurate time step: %.30f [s]\n",
147 estimate_time_step(&md,3.0,md.t));
150 * let's do the actual md algorithm now
152 * integration of newtons equations
155 moldyn_integrate(&md);
157 printf("total energy (after integration): %.40f [J]\n",
158 get_total_energy(&md));
162 link_cell_shutdown(&md);
164 moldyn_shutdown(&md);