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         moldyn_init(&md,argc,argv);
  38
  39         /*
  40          * the following overrides possibly set interaction methods by argv !!
  41          */
  42
  43         /* params */
  44         lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
  45         help=lj.sigma6*lj.sigma6;
  46         lj.sigma6*=help;
  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;
  51         /* assignement */
  52         md.potential_force_function=lennard_jones;
  53         //md.potential_force_function=harmonic_oscillator;
  54         md.pot_params=&lj;
  55         //md.pot_params=&ho;
  56         /* cutoff radius */
  57         md.cutoff=R_CUTOFF*LC_SI;
  58
  59         /*
  60          * testing random numbers
  61          */
  62
  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)));
  67         return 0;
  68 #endif
  69
  70         /*
  71          * geometry & particles
  72          */
  73
  74         /* simulation cell volume in lattice constants */
  75         a=LEN_X;
  76         b=LEN_Y;
  77         c=LEN_Z;
  78         md.dim.x=a*LC_SI;
  79         md.dim.y=b*LC_SI;
  80         md.dim.z=c*LC_SI;
  81
  82         /* (un)set to (not) get visualized 'bounding atoms' */
  83         md.vis.dim.x=a*LC_SI;
  84         md.vis.dim.y=b*LC_SI;
  85         md.vis.dim.z=c*LC_SI;
  86
  87         /*
  88          * particles
  89          */
  90
  91         /* lattice init */
  92
  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);
  96 #else
  97         md.count=2;
  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;
 109 #endif
 110
 111         /* initial thermal fluctuations of particles */
 112
 113 #ifndef SIMPLE_TESTING
 114         printf("setting thermal fluctuations (T=%f K)\n",md.t);
 115         thermal_init(&md);
 116 #else
 117         for(a=0;a<md.count;a++) v3_zero(&(si[0].v));
 118 #endif
 119
 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);
 125
 126         /* check total momentum */
 127         p=get_total_p(md.atom,md.count);
 128         printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
 129
 130         /* check time step */
 131         printf("estimated accurate time step: %.30f [s]\n",
 132                estimate_time_step(&md,3.0,md.t));
 133
 134         /* initialize linked list / cell method */
 135         printf("initializing linked cells\n");
 136         link_cell_init(&md);
 137
 138         /*
 139          * let's do the actual md algorithm now
 140          *
 141          * integration of newtons equations
 142          */
 143
 144         moldyn_integrate(&md);
 145
 146         printf("total energy (after integration): %.40f [J]\n",
 147                get_total_energy(&md));
 148
 149         /* close */
 150
 151         link_cell_shutdown(&md);
 152
 153         moldyn_shutdown(&md);
 154
 155         return 0;
 156 }
 157