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         t_atom *si;
  21         t_visual vis;
  22         t_random random;
  23
  24         int a,b,c;
  25         double e;
  26         double help;
  27         t_3dvec p;
  28
  29         t_lj_params lj;
  30         t_ho_params ho;
  31
  32         /* parse arguments */
  33         a=moldyn_parse_argv(&md,argc,argv);
  34         if(a<0) return -1;
  35
  36         /* init */
  37         moldyn_log_init(&md,&vis);
  38         rand_init(&random,NULL,1);
  39         random.status|=RAND_STAT_VERBOSE;
  40
  41         /* testing random numbers */
  42         //for(a=0;a<1000000;a++)
  43         //      printf("%f %f\n",rand_get_gauss(&random),
  44         //                       rand_get_gauss(&random));
  45
  46         a=LEN_X;
  47         b=LEN_Y;
  48         c=LEN_Z;
  49
  50         /* set for 'bounding atoms' */
  51         vis.dim.x=a*LC_SI;
  52         vis.dim.y=b*LC_SI;
  53         vis.dim.z=c*LC_SI;
  54
  55         /* init lattice
  56         printf("placing silicon atoms ... ");
  57         md.count=create_lattice(DIAMOND,SI,M_SI,LC_SI,a,b,c,&si);
  58         printf("(%d) ok!\n",md.count);
  59         testing purpose */
  60         md.count=2;
  61         si=malloc(2*sizeof(t_atom));
  62         si[0].r.x=0.13*sqrt(3.0)*LC_SI/2.0;
  63         si[0].r.y=0;
  64         si[0].r.z=0;
  65         si[0].element=SI;
  66         si[0].mass=M_SI;
  67         si[1].r.x=-si[0].r.x;
  68         si[1].r.y=0;
  69         si[1].r.z=0;
  70         si[1].element=SI;
  71         si[1].mass=M_SI;
  72         /* */
  73
  74         /* moldyn init (now si is a valid address) */
  75         md.atom=si;
  76         md.potential_force_function=lennard_jones;
  77         //md.potential_force_function=harmonic_oscillator;
  78         md.cutoff=R_CUTOFF*LC_SI;
  79         md.pot_params=&lj;
  80         //md.pot_params=&ho;
  81         md.status=0;
  82         md.visual=&vis;
  83         /* dimensions of the simulation cell */
  84         md.dim.x=a*LC_SI;
  85         md.dim.y=b*LC_SI;
  86         md.dim.z=c*LC_SI;
  87
  88         printf("setting thermal fluctuations (T=%f K)\n",md.t);
  89         // thermal_init(&md,&random);
  90         for(a=0;a<md.count;a++) v3_zero(&(si[0].v));
  91
  92         /* check kinetic energy */
  93
  94         e=get_e_kin(si,md.count);
  95         printf("kinetic energy: %.40f [J]\n",e);
  96         printf("3/2 N k T = %.40f [J] (T=%f [K])\n",
  97                1.5*md.count*K_BOLTZMANN*md.t,md.t);
  98
  99         /* check total momentum */
 100         p=get_total_p(si,md.count);
 101         printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
 102
 103         /* potential paramters */
 104         lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
 105         help=lj.sigma6*lj.sigma6;
 106         lj.sigma6*=help;
 107         lj.sigma12=lj.sigma6*lj.sigma6;
 108         lj.epsilon4=4.0*LJ_EPSILON_SI;
 109
 110         ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
 111         ho.spring_constant=1.0;
 112
 113         printf("estimated accurate time step: %.30f [s]\n",
 114                estimate_time_step(&md,3.0,md.t));
 115
 116         /*
 117          * let's do the actual md algorithm now
 118          *
 119          * integration of newtons equations
 120          */
 121
 122         moldyn_integrate(&md);
 123
 124         printf("total energy (after integration): %.40f [J]\n",
 125                get_total_energy(&md));
 126
 127         /* close */
 128
 129         link_cell_shutdown(&md);
 130
 131         rand_close(&random);
 132
 133         moldyn_shutdown(&md);
 134
 135         return 0;
 136 }
 137