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         t_atom *si;
  22
  23         t_visual vis;
  24
  25         t_random random;
  26
  27         int a,b,c;
  28         double t,e,u;
  29         double help;
  30         t_3dvec p;
  31         int count;
  32
  33         t_lj_params lj;
  34         t_ho_params ho;
  35
  36         char fb[32]="saves/lj_test";
  37
  38         /* init */
  39
  40         rand_init(&random,NULL,1);
  41         random.status|=RAND_STAT_VERBOSE;
  42
  43         /* testing random numbers */
  44         //for(a=0;a<1000000;a++)
  45         //      printf("%f %f\n",rand_get_gauss(&random),
  46         //                       rand_get_gauss(&random));
  47
  48         visual_init(&vis,fb);
  49
  50         a=LEN_X;
  51         b=LEN_Y;
  52         c=LEN_Z;
  53
  54         vis.dim.x=a*LC_SI;
  55         vis.dim.y=b*LC_SI;
  56         vis.dim.z=c*LC_SI;
  57
  58         t=TEMPERATURE;
  59
  60         printf("placing silicon atoms ... ");
  61         //count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si);
  62         //printf("(%d) ok!\n",count);
  63         count=2;
  64         si=malloc(2*sizeof(t_atom));
  65         si[0].r.x=0.16e-9;
  66         si[0].r.y=0;
  67         si[0].r.z=0;
  68         si[0].element=SI;
  69         si[0].mass=M_SI;
  70         si[1].r.x=-0.16e-9;
  71         si[1].r.y=0;
  72         si[1].r.z=0;
  73         si[1].element=SI;
  74         si[1].mass=M_SI;
  75
  76         printf("setting thermal fluctuations\n");
  77         //thermal_init(si,&random,count,t);
  78         v3_zero(&(si[0].v));
  79         v3_zero(&(si[1].v));
  80
  81         /* check kinetic energy */
  82
  83         e=get_e_kin(si,count);
  84         printf("kinetic energy: %.40f [J]\n",e);
  85         printf("3/2 N k T = %.40f [J]\n",1.5*count*K_BOLTZMANN*t);
  86
  87         /* check total momentum */
  88         p=get_total_p(si,count);
  89         printf("total momentum: %.30f [Ns]\n",v3_norm(&p));
  90
  91         /* check potential energy */
  92         md.count=count;
  93         md.atom=si;
  94         md.potential=potential_lennard_jones;
  95         //md.potential=potential_harmonic_oscillator;
  96         md.force=force_lennard_jones;
  97         //md.force=force_harmonic_oscillator;
  98         //md.cutoff_square=((LC_SI/4.0)*(LC_SI/4.0));
  99         md.cutoff=(0.4e-9);
 100         md.cutoff_square=(0.6e-9*0.6e-9);
 101         md.pot_params=&lj;
 102         //md.pot_params=&ho;
 103         md.integrate=velocity_verlet;
 104         md.time_steps=RUNS;
 105         md.tau=TAU;
 106         md.status=0;
 107         md.visual=&vis;
 108         md.write=WRITE_FILE;
 109
 110         lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI;
 111         help=lj.sigma6*lj.sigma6;
 112         lj.sigma6*=help;
 113         lj.sigma12=lj.sigma6*lj.sigma6;
 114         lj.epsilon=LJ_EPSILON_SI;
 115
 116         ho.equilibrium_distance=0.3e-9;
 117         ho.spring_constant=1.0e-9;
 118
 119         u=get_e_pot(&md);
 120
 121         printf("potential energy: %.40f [J]\n",u);
 122         printf("total energy (1): %.40f [J]\n",e+u);
 123         printf("total energy (2): %.40f [J]\n",get_total_energy(&md));
 124
 125         md.dim.x=a*LC_SI;
 126         md.dim.y=b*LC_SI;
 127         md.dim.z=c*LC_SI;
 128
 129         printf("estimated accurate time step: %.30f [s]\n",
 130                estimate_time_step(&md,3.0,t));
 131
 132
 133         /*
 134          * let's do the actual md algorithm now
 135          *
 136          * integration of newtons equations
 137          */
 138
 139         moldyn_integrate(&md);
 140
 141         printf("total energy (after integration): %.40f [J]\n",
 142                get_total_energy(&md));
 143
 144         /* close */
 145
 146         visual_tini(&vis);
 147
 148         rand_close(&random);
 149
 150
 151         //printf("starting velocity verlet: ");
 152         //fflush(stdout);
 153
 154         //for(runs=0;runs<RUNS;runs++) {
 155
 156         /*
 157          * velocity verlet
 158          *
 159          * r(t+h) = r(t) + h * dr/dt|t + h^2/2m * F(t)
 160          * dr/dt|(t+h) = dr/dt|t + h/2m * (F(t) + F(t+h))
 161          *
 162          */
 163         //for(i=0;i<amount_si;i++) {
 164 //              /* calculation of new positions r(t+h) */
 165 //              si[i].x+=si[i].vx*tau;
 166 //              si[i].y+=si[i].vy*tau;
 167 //              si[i].z+=si[i].vz*tau;
 168 //              si[i].x+=(tau2*si[i].fx/m2);
 169 //              if(si[i].x>LX) si[i].x-=LEN_X;
 170 //              else if(si[i].x<-LX) si[i].x+=LEN_X;
 171 //              si[i].y+=(tau2*si[i].fy/m2);
 172 //              if(si[i].y>LY) si[i].y-=LEN_Y;
 173 //              else if(si[i].y<-LY) si[i].y+=LEN_Y;
 174 //              si[i].z+=(tau2*si[i].fz/m2);
 175 //              if(si[i].z>LZ) si[i].z-=LEN_Z;
 176 //              else if(si[i].z<-LZ) si[i].z+=LEN_Z;
 177 //              /* calculation of velocities v(t+h/2) */
 178 //              si[i].vx+=(tau*si[i].fx/m2);
 179 //              si[i].vy+=(tau*si[i].fy/m2);
 180 //              si[i].vz+=(tau*si[i].fz/m2);
 181 //              /* reset of forces */
 182 //              si[i].fx=.0;
 183 //              si[i].fy=.0;
 184 //              si[i].fz=.0;
 185 //      }
 186 //      for(i=0;i<amount_si;i++) {
 187 //              /* calculation of forces at new positions r(t+h) */
 188 //              for(j=0;j<i;j++) {
 189 //                      deltax=si[i].x-si[j].x;
 190 //                      if(deltax>LX) deltax-=LEN_X;
 191 //                      else if(-deltax>LX) deltax+=LEN_X;
 192 //                      deltax2=deltax*deltax;
 193 //                      deltay=si[i].y-si[j].y;
 194 //                      if(deltay>LY) deltay-=LEN_Y;
 195 //                      else if(-deltay>LY) deltay+=LEN_Y;
 196 //                      deltay2=deltay*deltay;
 197 //                      deltaz=si[i].z-si[j].z;
 198 //                      if(deltaz>LZ) deltaz-=LEN_Z;
 199 //                      else if(-deltaz>LZ) deltaz+=LEN_Z;
 200 //                      deltaz2=deltaz*deltaz;
 201 //                      distance=deltax2+deltay2+deltaz2;
 202 //                      if(distance<=R2_CUTOFF) {
 203 //                              tmp=1.0/distance; // 1/r^2
 204 //                              lj1=tmp; // 1/r^2
 205 //                              tmp*=tmp; // 1/r^4
 206 //                              lj1*=tmp; // 1/r^6
 207 //                              tmp*=tmp; // 1/r^8
 208 //                              lj2=tmp; // 1/r^8
 209 //                              lj1*=tmp; // 1/r^14
 210 //                              lj1*=LJ_SIGMA_12;
 211 //                              lj2*=LJ_SIGMA_06;
 212 //                              lj=-2*lj1+lj2;
 213 //                              si[i].fx-=lj*deltax;
 214 //                              si[i].fy-=lj*deltay;
 215 //                              si[i].fz-=lj*deltaz;
 216 //                              si[j].fx+=lj*deltax;
 217 //                              si[j].fy+=lj*deltay;
 218 //                              si[j].fz+=lj*deltaz;
 219 //                      }
 220 //              }
 221 //      }
 222 //      for(i=0;i<amount_si;i++) {
 223 //              /* calculation of new velocities v(t+h) */
 224 //              si[i].vx+=(tau*si[i].fx/m2);
 225 //              si[i].vy+=(tau*si[i].fy/m2);
 226 //              si[i].vz+=(tau*si[i].fz/m2);
 227 //      }
 228 //
 229 //      if(!(runs%150)) {
 230 //
 231 //      /* rasmol script & xyz file */
 232 //      sprintf(xyz,"./saves/si-%.15f.xyz",time);
 233 //      sprintf(ppm,"./video/si-%.15f.ppm",time);
 234 //      fd1=open(xyz,O_WRONLY|O_CREAT|O_TRUNC);
 235 //      if(fd1<0) {
 236 //              perror("rasmol xyz file open");
 237 //              return -1;
 238 //      }
 239 //      dprintf(fd2,"load xyz %s\n",xyz);
 240 //      dprintf(fd2,"spacefill 200\n");
 241 //      dprintf(fd2,"rotate x 11\n");
 242 //      dprintf(fd2,"rotate y 13\n");
 243 //      dprintf(fd2,"set ambient 20\n");
 244 //      dprintf(fd2,"set specular on\n");
 245 //      dprintf(fd2,"zoom 400\n");
 246 //      dprintf(fd2,"write ppm %s\n",ppm);
 247 //      dprintf(fd2,"zap\n");
 248 //      dprintf(fd1,"%d\nsilicon\n",amount_si+9);
 249 //      for(i=0;i<amount_si;i++)
 250 //              dprintf(fd1,"Si %f %f %f %f\n",
 251 //                      si[i].x,si[i].y,si[i].z,time);
 252 //      dprintf(fd1,"H 0.0 0.0 0.0 %f\n",time);
 253 //      dprintf(fd1,"He %f %f %f %f\n",LX,LY,LZ,time);
 254 //      dprintf(fd1,"He %f %f %f %f\n",-LX,LY,LZ,time);
 255 //      dprintf(fd1,"He %f %f %f %f\n",LX,-LY,LZ,time);
 256 //      dprintf(fd1,"He %f %f %f %f\n",LX,LY,-LZ,time);
 257 //      dprintf(fd1,"He %f %f %f %f\n",-LX,-LY,LZ,time);
 258 //      dprintf(fd1,"He %f %f %f %f\n",-LX,LY,-LZ,time);
 259 //      dprintf(fd1,"He %f %f %f %f\n",LX,-LY,-LZ,time);
 260 //      dprintf(fd1,"He %f %f %f %f\n",-LX,-LY,-LZ,time);
 261 //      close(fd1);
 262 //
 263 //      }
 264 //
 265 //      /* increase time */
 266 //      time+=tau;
 267 //      printf(".");
 268 //      fflush(stdout);
 269 //
 270 //      }
 271 //
 272 //      printf(" done\n");
 273 //      close(fd2);
 274 //      free(si);
 275 //
 276
 277         return 0;
 278 }
 279