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