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) {
23 t_tersoff_mult_params tp;
38 a=moldyn_init(&md,argc,argv);
42 * the following overrides possibly set interaction methods by argv !!
46 lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
47 help=lj.sigma6*lj.sigma6;
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;
54 md.potential_force_function=lennard_jones;
55 //md.potential_force_function=harmonic_oscillator;
59 md.cutoff=R_CUTOFF*LC_SI;
62 * testing random numbers
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)));
73 * geometry & particles
76 /* simulation cell volume in lattice constants */
84 /* (un)set to (not) get visualized 'bounding atoms' */
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);
100 md.atom=malloc(md.count*sizeof(t_atom));
101 md.atom[0].r.x=0.23*sqrt(3.0)*LC_SI/2.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;
109 md.atom[1].element=SI;
110 md.atom[1].mass=M_SI;
112 //md.atom[2].r.x=0.5*(a-1)*LC_SI;
113 //md.atom[2].r.y=0.5*(b-1)*LC_SI;
115 //md.atom[2].element=C;
116 //md.atom[2].mass=M_C;
118 //md.atom[3].r.x=0.5*(a-1)*LC_SI;
121 //md.atom[3].element=SI;
122 //md.atom[3].mass=M_SI;
125 /* initial thermal fluctuations of particles */
127 #ifndef SIMPLE_TESTING
128 printf("setting thermal fluctuations (T=%f K)\n",md.t);
131 for(a=0;a<md.count;a++) v3_zero(&(md.atom[0].v));
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);
142 /* check total momentum */
143 p=get_total_p(md.atom,md.count);
144 printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
146 /* check time step */
147 printf("estimated accurate time step: %.30f [s]\n",
148 estimate_time_step(&md,3.0,md.t));
151 * let's do the actual md algorithm now
153 * integration of newtons equations
156 moldyn_integrate(&md);
158 printf("total energy (after integration): %.40f [J]\n",
159 get_total_energy(&md));
163 link_cell_shutdown(&md);
165 moldyn_shutdown(&md);