posic.c

   1 /*
   2  * posic.c - precipitation process of silicon carbide in silicon
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #include <math.h>
   9
  10 #include "moldyn.h"
  11 #include "math/math.h"
  12 #include "init/init.h"
  13 #include "visual/visual.h"
  14
  15 #include "posic.h"
  16
  17 int main(int argc,char **argv) {
  18
  19         t_moldyn md;
  20
  21         t_lj_params lj;
  22         t_ho_params ho;
  23         t_tersoff_mult_params tp;
  24
  25         int a,b,c;
  26         double e;
  27         double help;
  28         t_3dvec p;
  29
  30         /*
  31          *  moldyn init
  32          *
  33          * - parsing argv
  34          * - log init
  35          * - random init
  36          *
  37          */
  38         a=moldyn_init(&md,argc,argv);
  39         if(a<0) return a;
  40
  41         /*
  42          * the following overrides possibly set interaction methods by argv !!
  43          */
  44
  45         /* params */
  46         lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
  47         help=lj.sigma6*lj.sigma6;
  48         lj.sigma6*=help;
  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;
  52         ho.spring_constant=1;
  53         /* assignement */
  54         md.potential_force_function=lennard_jones;
  55         //md.potential_force_function=harmonic_oscillator;
  56         md.pot_params=&lj;
  57         //md.pot_params=&ho;
  58         /* cutoff radius */
  59         md.cutoff=R_CUTOFF*LC_SI;
  60
  61         /*
  62          * testing random numbers
  63          */
  64
  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)));
  69         return 0;
  70 #endif
  71
  72         /*
  73          * geometry & particles
  74          */
  75
  76         /* simulation cell volume in lattice constants */
  77         a=LEN_X;
  78         b=LEN_Y;
  79         c=LEN_Z;
  80         md.dim.x=a*LC_SI;
  81         md.dim.y=b*LC_SI;
  82         md.dim.z=c*LC_SI;
  83
  84         /* (un)set to (not) get visualized 'bounding atoms' */
  85         md.vis.dim.x=a*LC_SI;
  86         md.vis.dim.y=b*LC_SI;
  87         md.vis.dim.z=c*LC_SI;
  88
  89         /*
  90          * particles
  91          */
  92
  93         /* lattice init */
  94
  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);
  98 #else
  99         md.count=2;
 100         md.atom=malloc(md.count*sizeof(t_atom));
 101         md.atom[0].r.x=0.23*sqrt(3.0)*LC_SI/2.0;
 102         md.atom[0].r.y=0;
 103         md.atom[0].r.z=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;
 107         md.atom[1].r.y=0;
 108         md.atom[1].r.z=0;
 109         md.atom[1].element=SI;
 110         md.atom[1].mass=M_SI;
 111
 112         //md.atom[2].r.x=0.5*(a-1)*LC_SI;
 113         //md.atom[2].r.y=0.5*(b-1)*LC_SI;
 114         //md.atom[2].r.z=0;
 115         //md.atom[2].element=C;
 116         //md.atom[2].mass=M_C;
 117
 118         //md.atom[3].r.x=0.5*(a-1)*LC_SI;
 119         //md.atom[3].r.y=0;
 120         //md.atom[3].r.z=0;
 121         //md.atom[3].element=SI;
 122         //md.atom[3].mass=M_SI;
 123 #endif
 124
 125         /* initial thermal fluctuations of particles */
 126
 127 #ifndef SIMPLE_TESTING
 128         printf("setting thermal fluctuations (T=%f K)\n",md.t);
 129         thermal_init(&md);
 130 #else
 131         for(a=0;a<md.count;a++) v3_zero(&(md.atom[0].v));
 132         md.atom[2].v.x=-320;
 133         md.atom[2].v.y=-320;
 134 #endif
 135
 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);
 141
 142         /* check total momentum */
 143         p=get_total_p(md.atom,md.count);
 144         printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
 145
 146         /* check time step */
 147         printf("estimated accurate time step: %.30f [s]\n",
 148                estimate_time_step(&md,3.0,md.t));
 149
 150         /*
 151          * let's do the actual md algorithm now
 152          *
 153          * integration of newtons equations
 154          */
 155
 156         moldyn_integrate(&md);
 157
 158         printf("total energy (after integration): %.40f [J]\n",
 159                get_total_energy(&md));
 160
 161         /* close */
 162
 163         link_cell_shutdown(&md);
 164
 165         moldyn_shutdown(&md);
 166
 167         return 0;
 168 }
 169