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 #include <math.h>
  13
  14 #include "math/math.h"
  15 #include "init/init.h"
  16 #include "random/random.h"
  17
  18
  19 int create_lattice(unsigned char type,int element,double mass,double lc,
  20                    int a,int b,int c,t_atom **atom) {
  21
  22         int count;
  23         int ret;
  24         t_3dvec origin;
  25
  26         count=a*b*c;
  27
  28         if(type==FCC) count*=4;
  29         if(type==DIAMOND) count*=8;
  30
  31         *atom=malloc(count*sizeof(t_atom));
  32         if(*atom==NULL) {
  33                 perror("malloc (atoms)");
  34                 return -1;
  35         }
  36
  37         v3_zero(&origin);
  38
  39         switch(type) {
  40                 case FCC:
  41                         ret=fcc_init(a,b,c,lc,*atom,&origin);
  42                         break;
  43                 case DIAMOND:
  44                         ret=diamond_init(a,b,c,lc,*atom,&origin);
  45                         break;
  46                 default:
  47                         printf("unknown lattice type (%02x)\n",type);
  48                         return -1;
  49         }
  50
  51         /* debug */
  52         if(ret!=count) {
  53                 printf("ok, there is something wrong ...\n");
  54                 printf("calculated -> %d atoms\n",count);
  55                 printf("created -> %d atoms\n",ret);
  56                 return -1;
  57         }
  58
  59         while(count) {
  60                 (*atom)[count-1].element=element;
  61                 (*atom)[count-1].mass=mass;
  62                 count-=1;
  63         }
  64
  65         return ret;
  66 }
  67
  68 int destroy_lattice(t_atom *atom) {
  69
  70         if(atom) free(atom);
  71
  72         return 0;
  73 }
  74
  75 int thermal_init(t_atom *atom,t_random *random,int count,double t) {
  76
  77         /*
  78          * - gaussian distribution of velocities
  79          * - zero total momentum
  80          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
  81          */
  82
  83         int i;
  84         double v,sigma;
  85         t_3dvec p_total,delta;
  86
  87         /* gaussian distribution of velocities */
  88         v3_zero(&p_total);
  89         for(i=0;i<count;i++) {
  90                 sigma=sqrt(2.0*K_BOLTZMANN*t/atom[i].mass);
  91                 /* x direction */
  92                 v=sigma*rand_get_gauss(random);
  93                 atom[i].v.x=v;
  94                 p_total.x+=atom[i].mass*v;
  95                 /* y direction */
  96                 v=sigma*rand_get_gauss(random);
  97                 atom[i].v.y=v;
  98                 p_total.y+=atom[i].mass*v;
  99                 /* z direction */
 100                 v=sigma*rand_get_gauss(random);
 101                 atom[i].v.z=v;
 102                 p_total.z+=atom[i].mass*v;
 103         }
 104
 105         /* zero total momentum */
 106         v3_scale(&p_total,&p_total,1.0/count);
 107         for(i=0;i<count;i++) {
 108                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 109                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 110         }
 111
 112         /* velocity scaling */
 113         scale_velocity(atom,count,t);
 114
 115         return 0;
 116 }
 117
 118 int scale_velocity(t_atom *atom,int count,double t) {
 119
 120         int i;
 121         double e,c;
 122
 123         /*
 124          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 125          */
 126         e=0.0;
 127         for(i=0;i<count;i++)
 128                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 129         c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*t));
 130         for(i=0;i<count;i++)
 131                 v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
 132
 133         return 0;
 134 }
 135
 136 double get_e_kin(t_atom *atom,int count) {
 137
 138         int i;
 139         double e;
 140
 141         e=0.0;
 142
 143         for(i=0;i<count;i++) {
 144                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 145         }
 146
 147         return e;
 148 }
 149
 150 double get_e_pot(t_moldyn *moldyn) {
 151
 152         return(moldyn->potential(moldyn));
 153 }
 154
 155 double get_total_energy(t_moldyn *moldyn) {
 156
 157         double e;
 158
 159         e=get_e_kin(moldyn->atom,moldyn->count);
 160         e+=get_e_pot(moldyn);
 161
 162         return e;
 163 }
 164
 165 t_3dvec get_total_p(t_atom *atom, int count) {
 166
 167         t_3dvec p,p_total;
 168         int i;
 169
 170         v3_zero(&p_total);
 171         for(i=0;i<count;i++) {
 172                 v3_scale(&p,&(atom[i].v),atom[i].mass);
 173                 v3_add(&p_total,&p_total,&p);
 174         }
 175
 176         return p_total;
 177 }
 178
 179
 180 /*
 181  *
 182  * potentials & corresponding forces
 183  *
 184  */
 185
 186 /* lennard jones potential & force for one sort of atoms */
 187
 188 double potential_lennard_jones(t_moldyn *moldyn) {
 189
 190         t_lj_params *params;
 191         t_atom *atom;
 192         int i,j;
 193         int count;
 194         t_3dvec distance;
 195         double d,help;
 196         double u;
 197         double eps,sig6,sig12;
 198
 199         params=moldyn->pot_params;
 200         atom=moldyn->atom;
 201         count=moldyn->count;
 202         eps=params->epsilon;
 203         sig6=params->sigma6;
 204         sig12=params->sigma12;
 205
 206         u=0.0;
 207         for(i=0;i<count;i++) {
 208                 for(j=0;j<i;j++) {
 209                         v3_sub(&distance,&(atom[j].r),&(atom[i].r));
 210                         d=1.0/v3_absolute_square(&distance);    /* 1/r^2 */
 211                         help=d*d;                               /* 1/r^4 */
 212                         help*=d;                                /* 1/r^6 */
 213                         d=help*help;                            /* 1/r^12 */
 214                         u+=eps*(sig12*d-sig6*help);
 215                 }
 216         }
 217
 218         return u;
 219 }
 220
 221 int force_lennard_jones(t_moldyn *moldyn) {
 222
 223         t_lj_params *params;
 224         int i,j,count;
 225         t_atom *atom;
 226         t_3dvec distance;
 227         t_3dvec force;
 228         double d,h1,h2;
 229
 230         atom=moldyn->atom;
 231         count=moldyn->count;
 232         params=moldyn->pot_params;
 233
 234         for(i=0;i<count;i++) {
 235                 for(j=0;j<i;j++) {
 236                         v3_sub(&distance,&(atom[j].r),&(atom[i].r));
 237                         v3_per_bound(&distance,&(moldyn->dim));
 238                         d=v3_absolute_square(&distance);
 239                         if(d<=moldyn->cutoff_square) {
 240                                 h1=1.0/d;                       /* 1/r^2 */
 241                                 d=h1*h1;                        /* 1/r^4 */
 242                                 h2=d*d;                         /* 1/r^8 */
 243                                 h1*=d;                          /* 1/r^6 */
 244                                 h1*=h2;                         /* 1/r^14 */
 245                         }
 246                 }
 247         }
 248
 249         return 0;
 250 }
 251