init/init.c

   1 /*
   2  * init.c - create initial conditions
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #include "../math/math.h"
   9 #include "../moldyn.h"
  10
  11 /* fcc lattice init */
  12 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
  13
  14         int count;
  15         int i,j;
  16         t_3dvec o,r,n;
  17         t_3dvec basis[3];
  18         double help[3];
  19         double x,y,z;
  20
  21         x=a*lc;
  22         y=b*lc;
  23         z=c*lc;
  24
  25         if(origin) v3_copy(&o,origin);
  26         else v3_zero(&o);
  27
  28         /* construct the basis */
  29         for(i=0;i<3;i++) {
  30                 for(j=0;j<3;j++) {
  31                         if(i!=j) help[j]=0.5*lc;
  32                         else help[j]=.0;
  33                 }
  34                 v3_set(&basis[i],help);
  35         }
  36
  37         v3_zero(&r);
  38         count=0;
  39
  40         /* fill up the room */
  41         r.x=o.x;
  42         while(r.x<x) {
  43                 r.y=o.y;
  44                 while(r.y<y) {
  45                         r.z=o.z;
  46                         while(r.z<z) {
  47                                 v3_copy(&(atom[count].r),&r);
  48                                 atom[count].element=1;
  49                                 count+=1;
  50                                 for(i=0;i<3;i++) {
  51                                         v3_add(&n,&r,&basis[i]);
  52                                         if((n.x<x+o.x)&&(n.y<y+o.y)&&(n.z<z+o.z)) {
  53                                                 v3_copy(&(atom[count].r),&n);
  54                                                 count+=1;
  55                                         }
  56                                 }
  57                                 r.z+=lc;
  58                         }
  59                         r.y+=lc;
  60                 }
  61                 r.x+=lc;
  62         }
  63
  64         /* coordinate transformation */
  65         help[0]=x/2.0;
  66         help[1]=y/2.0;
  67         help[2]=z/2.0;
  68         v3_set(&n,help);
  69         for(i=0;i<count;i++)
  70                 v3_sub(&(atom[i].r),&(atom[i].r),&n);
  71
  72         return count;
  73 }
  74
  75 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
  76
  77         int count;
  78         t_3dvec o;
  79
  80         count=fcc_init(a,b,c,lc,atom,origin);
  81
  82         o.x=0.25*lc;
  83         o.y=0.25*lc;
  84         o.z=0.25*lc;
  85
  86         if(origin) v3_add(&o,&o,origin);
  87
  88         count+=fcc_init(a,b,c,lc,&atom[count],&o);
  89
  90         return count;
  91 }