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;
51 ho.spring_constant=1.0;
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(2*sizeof(t_atom));
100 md.atom[0].r.x=0.13*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;
112 /* initial thermal fluctuations of particles */
114 #ifndef SIMPLE_TESTING
115 printf("setting thermal fluctuations (T=%f K)\n",md.t);
118 for(a=0;a<md.count;a++) v3_zero(&(md.atom[0].v));
121 /* check kinetic energy */
122 e=get_e_kin(md.atom,md.count);
123 printf("kinetic energy: %.40f [J]\n",e);
124 printf("3/2 N k T = %.40f [J] (T=%f [K])\n",
125 1.5*md.count*K_BOLTZMANN*md.t,md.t);
127 /* check total momentum */
128 p=get_total_p(md.atom,md.count);
129 printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
131 /* check time step */
132 printf("estimated accurate time step: %.30f [s]\n",
133 estimate_time_step(&md,3.0,md.t));
135 /* initialize linked list / cell method */
139 * let's do the actual md algorithm now
141 * integration of newtons equations
144 moldyn_integrate(&md);
146 printf("total energy (after integration): %.40f [J]\n",
147 get_total_energy(&md));
151 link_cell_shutdown(&md);
153 moldyn_shutdown(&md);