/*
* posic.c - precipitation process of silicon carbide in silicon
*
- * author: Frank Zirkelbach <hackbard@hackdaworld.org>
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
*
*/
+
+#include <math.h>
#include "moldyn.h"
#include "math/math.h"
int count;
t_lj_params lj;
+ t_ho_params ho;
- char fb[32]="saves/fcc_test";
+ char fb[32]="saves/lj_test";
/* init */
b=LEN_Y;
c=LEN_Z;
+ vis.dim.x=a*LC_SI;
+ vis.dim.y=b*LC_SI;
+ vis.dim.z=c*LC_SI;
+
t=TEMPERATURE;
printf("placing silicon atoms ... ");
//printf("(%d) ok!\n",count);
count=2;
si=malloc(2*sizeof(t_atom));
- si[0].r.x=2.0;
+ si[0].r.x=0.16e-9;
si[0].r.y=0;
si[0].r.z=0;
- si[0].element=Si;
- si[0].mass=14.0;
- si[1].r.x=-2.0;
+ si[0].element=SI;
+ si[0].mass=M_SI;
+ si[1].r.x=-0.16e-9;
si[1].r.y=0;
si[1].r.z=0;
- si[1].element=Si;
- si[1].mass=14.0;
+ si[1].element=SI;
+ si[1].mass=M_SI;
printf("setting thermal fluctuations\n");
//thermal_init(si,&random,count,t);
/* check kinetic energy */
e=get_e_kin(si,count);
- printf("kinetic energy: %f\n",e);
- printf("3/2 N k T = %f\n",1.5*count*K_BOLTZMANN*t);
+ printf("kinetic energy: %.40f [J]\n",e);
+ printf("3/2 N k T = %.40f [J]\n",1.5*count*K_BOLTZMANN*t);
/* check total momentum */
p=get_total_p(si,count);
- printf("total momentum: %f\n",v3_norm(&p));
+ printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
/* check potential energy */
md.count=count;
md.atom=si;
md.potential=potential_lennard_jones;
+ //md.potential=potential_harmonic_oscillator;
md.force=force_lennard_jones;
+ //md.force=force_harmonic_oscillator;
//md.cutoff_square=((LC_SI/4.0)*(LC_SI/4.0));
- md.cutoff_square=36.0;
+ md.cutoff=(0.4e-9);
+ md.cutoff_square=(0.6e-9*0.6e-9);
md.pot_params=&lj;
+ //md.pot_params=&ho;
md.integrate=velocity_verlet;
md.time_steps=RUNS;
md.tau=TAU;
md.status=0;
md.visual=&vis;
+ md.write=WRITE_FILE;
- lj.sigma6=3.0/16.0*LC_SI*LC_SI;
+ lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
help=lj.sigma6*lj.sigma6;
lj.sigma6*=help;
lj.sigma12=lj.sigma6*lj.sigma6;
- lj.epsilon=10000;
+ lj.epsilon=LJ_EPSILON_SI;
+
+ ho.equilibrium_distance=0.3e-9;
+ ho.spring_constant=1.0e-9;
u=get_e_pot(&md);
- printf("potential energy: %f\n",u);
- printf("total energy (1): %f\n",e+u);
- printf("total energy (2): %f\n",get_total_energy(&md));
+ printf("potential energy: %.40f [J]\n",u);
+ printf("total energy (1): %.40f [J]\n",e+u);
+ printf("total energy (2): %.40f [J]\n",get_total_energy(&md));
md.dim.x=a*LC_SI;
md.dim.y=b*LC_SI;
md.dim.z=c*LC_SI;
+ printf("estimated accurate time step: %.30f [s]\n",
+ estimate_time_step(&md,3.0,t));
+
+
/*
* let's do the actual md algorithm now
*
* integration of newtons equations
*/
- /* visualize */
- //visual_atoms(&vis,0.0,si,count);
-
-
moldyn_integrate(&md);
- printf("total energy (after integration): %f\n",get_total_energy(&md));
+ printf("total energy (after integration): %.40f [J]\n",
+ get_total_energy(&md));
/* close */