moldyn.c

   1 /*
   2  * moldyn.c - molecular dynamics library main file
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #include "moldyn.h"
   9
  10 #include <stdio.h>
  11 #include <stdlib.h>
  12
  13 #include "math/math.h"
  14 #include "init/init.h"
  15
  16
  17 int create_lattice(unsigned char type,int element,double mass,double lc,
  18                    int a,int b,int c,t_atom **atom) {
  19
  20         int count;
  21         int ret;
  22         t_3dvec origin;
  23
  24         count=a*b*c;
  25
  26         if(type==FCC) count*=4;
  27         if(type==DIAMOND) count*=8;
  28
  29         *atom=malloc(count*sizeof(t_atom));
  30         if(*atom==NULL) {
  31                 perror("malloc (atoms)");
  32                 return -1;
  33         }
  34
  35         v3_zero(&origin);
  36
  37         switch(type) {
  38                 case FCC:
  39                         ret=fcc_init(a,b,c,lc,*atom,&origin);
  40                         break;
  41                 case DIAMOND:
  42                         ret=diamond_init(a,b,c,lc,*atom,&origin);
  43                         break;
  44                 default:
  45                         printf("unknown lattice type (%02x)\n",type);
  46                         return -1;
  47         }
  48
  49         /* debug */
  50         if(ret!=count) {
  51                 printf("ok, there is something wrong ...\n");
  52                 printf("calculated -> %d atoms\n",count);
  53                 printf("created -> %d atoms\n",ret);
  54                 return -1;
  55         }
  56
  57         while(count) {
  58                 (*atom)[count-1].element=element;
  59                 (*atom)[count-1].mass=mass;
  60                 count-=1;
  61         }
  62
  63         return ret;
  64 }
  65
  66 int thermal_init(t_atom *atom,int count,double t) {
  67
  68         /*
  69          * - gaussian distribution of velocities
  70          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
  71          */
  72
  73         int i;
  74         double e,c,v;
  75
  76         e=.0;
  77
  78         for(i=0;i<count;i++) {
  79                 /* x direction */
  80                 v=gauss_rand();
  81                 atom[count].v.x=v;
  82                 e+=0.5*atom[count].mass*v*v;
  83                 /* y direction */
  84                 v=gauss_rand();
  85                 atom[count].v.y=v;
  86                 e+=0.5*atom[count].mass*v*v;
  87                 /* z direction */
  88                 v=gauss_rand();
  89                 atom[count].v.z=v;
  90                 e+=0.5*atom[count].mass*v*v;
  91         }
  92
  93         c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN));
  94
  95         for(i=0;i<count;i++)
  96                 v3_scale(&(atom[count].v),&(atom[count].v),(1.0/c));
  97
  98         return 0;
  99 }
 100