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 moldyn_init(&md,argc,argv);
40 * the following overrides possibly set interaction methods by argv !!
44 lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
45 help=lj.sigma6*lj.sigma6;
47 lj.sigma12=lj.sigma6*lj.sigma6;
48 lj.epsilon4=4.0*LJ_EPSILON_SI;
49 ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
50 ho.spring_constant=1.0;
52 md.potential_force_function=lennard_jones;
53 //md.potential_force_function=harmonic_oscillator;
57 md.cutoff=R_CUTOFF*LC_SI;
60 * testing random numbers
63 #ifdef DEBUG_RANDOM_NUMBER
64 for(a=0;a<1000000;a++)
65 printf("%f %f\n",rand_get_gauss(&(md.random)),
66 rand_get_gauss(&(md.random)));
71 * geometry & particles
74 /* simulation cell volume in lattice constants */
82 /* (un)set to (not) get visualized 'bounding atoms' */
93 #ifndef SIMPLE_TESTING
94 md.count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&(md.atom));
95 printf("created silicon lattice (#atoms = %d)\n",md.count);
98 moldyn->atom=malloc(2*sizeof(t_atom));
99 moldyn->atom[0].r.x=0.13*sqrt(3.0)*LC_SI/2.0;
100 moldyn->atom[0].r.y=0;
101 moldyn->atom[0].r.z=0;
102 moldyn->atom[0].element=SI;
103 moldyn->atom[0].mass=M_SI;
104 moldyn->atom[1].r.x=-si[0].r.x;
105 moldyn->atom[1].r.y=0;
106 moldyn->atom[1].r.z=0;
107 moldyn->atom[1].element=SI;
108 moldyn->atom[1].mass=M_SI;
111 /* initial thermal fluctuations of particles */
113 #ifndef SIMPLE_TESTING
114 printf("setting thermal fluctuations (T=%f K)\n",md.t);
117 for(a=0;a<md.count;a++) v3_zero(&(si[0].v));
120 /* check kinetic energy */
121 e=get_e_kin(md.atom,md.count);
122 printf("kinetic energy: %.40f [J]\n",e);
123 printf("3/2 N k T = %.40f [J] (T=%f [K])\n",
124 1.5*md.count*K_BOLTZMANN*md.t,md.t);
126 /* check total momentum */
127 p=get_total_p(md.atom,md.count);
128 printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
130 /* check time step */
131 printf("estimated accurate time step: %.30f [s]\n",
132 estimate_time_step(&md,3.0,md.t));
134 /* initialize linked list / cell method */
135 printf("initializing linked cells\n");
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);