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         v3_add(&r,&r,&o);
  39         count=0;
  40
  41         /* fill up the room */
  42         while(r.x<x) {
  43                 r.y=.0;
  44                 while(r.y<y) {
  45                         r.z=.0;
  46                         while(r.z<z) {
  47                                 v3_copy(&(atom[count].r),&r);
  48                                 count+=1;
  49                                 for(i=0;i<3;i++) {
  50                                         v3_add(&n,&r,&basis[i]);
  51                                         if((n.x<x+o.x)&&(n.y<y+o.y)&&(n.z<z+o.z)) {
  52                                                 v3_copy(&(atom[count].r),&n);
  53                                                 count+=1;
  54                                         }
  55                                 }
  56                                 r.z+=lc;
  57                         }
  58                         r.y+=lc;
  59                 }
  60                 r.x+=lc;
  61         }
  62
  63         /* coordinate transformation */
  64         help[0]=x/2.0;
  65         help[1]=y/2.0;
  66         help[2]=z/2.0;
  67         v3_set(&n,help);
  68         for(i=0;i<count;i++)
  69                 v3_sub(&(atom[i].r),&(atom[i].r),&n);
  70
  71         return count;
  72 }
  73
  74 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
  75
  76         int count;
  77         t_3dvec o;
  78
  79         count=fcc_init(a,b,c,lc,atom,origin);
  80
  81         o.x=0.25*lc;
  82         o.y=0.25*lc;
  83         o.z=0.25*lc;
  84
  85         if(origin) v3_add(&o,&o,origin);
  86
  87         count+=fcc_init(a,b,c,lc,&atom[count],&o);
  88
  89         return count;
  90 }