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         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.13*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*LC_SI;
  83         md.cutoff_square=md.cutoff*md.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         /* dimensions of the simulation cell */
  92         md.dim.x=a*LC_SI;
  93         md.dim.y=b*LC_SI;
  94         md.dim.z=c*LC_SI;
  95
  96         /* verlet list init */
  97         // later integrated in moldyn_init function!
  98         verlet_list_init(&md);
  99
 100         printf("setting thermal fluctuations (T=%f K)\n",md.t);
 101         thermal_init(&md,&random,count);
 102         //for(a=0;a<count;a++) v3_zero(&(si[0].v));
 103
 104         /* check kinetic energy */
 105
 106         e=get_e_kin(si,count);
 107         printf("kinetic energy: %.40f [J]\n",e);
 108         printf("3/2 N k T = %.40f [J] (T=%f [K])\n",
 109                1.5*count*K_BOLTZMANN*md.t,md.t);
 110
 111         /* check total momentum */
 112         p=get_total_p(si,count);
 113         printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
 114
 115         /* potential paramters */
 116         lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
 117         help=lj.sigma6*lj.sigma6;
 118         lj.sigma6*=help;
 119         lj.sigma12=lj.sigma6*lj.sigma6;
 120         lj.epsilon4=4.0*LJ_EPSILON_SI;
 121
 122         ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
 123         ho.spring_constant=1.0;
 124
 125         printf("estimated accurate time step: %.30f [s]\n",
 126                estimate_time_step(&md,3.0,md.t));
 127
 128
 129         /*
 130          * let's do the actual md algorithm now
 131          *
 132          * integration of newtons equations
 133          */
 134
 135         moldyn_integrate(&md);
 136
 137         printf("total energy (after integration): %.40f [J]\n",
 138                get_total_energy(&md));
 139
 140         /* close */
 141
 142         verlet_list_shutdown(&md);
 143
 144         rand_close(&random);
 145
 146         moldyn_shutdown(&md);
 147
 148         return 0;
 149 }
 150