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