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 #include "visual/visual.h"
  18
  19
  20 int create_lattice(unsigned char type,int element,double mass,double lc,
  21                    int a,int b,int c,t_atom **atom) {
  22
  23         int count;
  24         int ret;
  25         t_3dvec origin;
  26
  27         count=a*b*c;
  28
  29         if(type==FCC) count*=4;
  30         if(type==DIAMOND) count*=8;
  31
  32         *atom=malloc(count*sizeof(t_atom));
  33         if(*atom==NULL) {
  34                 perror("malloc (atoms)");
  35                 return -1;
  36         }
  37
  38         v3_zero(&origin);
  39
  40         switch(type) {
  41                 case FCC:
  42                         ret=fcc_init(a,b,c,lc,*atom,&origin);
  43                         break;
  44                 case DIAMOND:
  45                         ret=diamond_init(a,b,c,lc,*atom,&origin);
  46                         break;
  47                 default:
  48                         printf("unknown lattice type (%02x)\n",type);
  49                         return -1;
  50         }
  51
  52         /* debug */
  53         if(ret!=count) {
  54                 printf("ok, there is something wrong ...\n");
  55                 printf("calculated -> %d atoms\n",count);
  56                 printf("created -> %d atoms\n",ret);
  57                 return -1;
  58         }
  59
  60         while(count) {
  61                 (*atom)[count-1].element=element;
  62                 (*atom)[count-1].mass=mass;
  63                 count-=1;
  64         }
  65
  66         return ret;
  67 }
  68
  69 int destroy_lattice(t_atom *atom) {
  70
  71         if(atom) free(atom);
  72
  73         return 0;
  74 }
  75
  76 int thermal_init(t_atom *atom,t_random *random,int count,double t) {
  77
  78         /*
  79          * - gaussian distribution of velocities
  80          * - zero total momentum
  81          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
  82          */
  83
  84         int i;
  85         double v,sigma;
  86         t_3dvec p_total,delta;
  87
  88         /* gaussian distribution of velocities */
  89         v3_zero(&p_total);
  90         for(i=0;i<count;i++) {
  91                 sigma=sqrt(2.0*K_BOLTZMANN*t/atom[i].mass);
  92                 /* x direction */
  93                 v=sigma*rand_get_gauss(random);
  94                 atom[i].v.x=v;
  95                 p_total.x+=atom[i].mass*v;
  96                 /* y direction */
  97                 v=sigma*rand_get_gauss(random);
  98                 atom[i].v.y=v;
  99                 p_total.y+=atom[i].mass*v;
 100                 /* z direction */
 101                 v=sigma*rand_get_gauss(random);
 102                 atom[i].v.z=v;
 103                 p_total.z+=atom[i].mass*v;
 104         }
 105
 106         /* zero total momentum */
 107         v3_scale(&p_total,&p_total,1.0/count);
 108         for(i=0;i<count;i++) {
 109                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 110                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 111         }
 112
 113         /* velocity scaling */
 114         scale_velocity(atom,count,t);
 115
 116         return 0;
 117 }
 118
 119 int scale_velocity(t_atom *atom,int count,double t) {
 120
 121         int i;
 122         double e,c;
 123
 124         /*
 125          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 126          */
 127         e=0.0;
 128         for(i=0;i<count;i++)
 129                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 130         c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*t));
 131         for(i=0;i<count;i++)
 132                 v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
 133
 134         return 0;
 135 }
 136
 137 double get_e_kin(t_atom *atom,int count) {
 138
 139         int i;
 140         double e;
 141
 142         e=0.0;
 143
 144         for(i=0;i<count;i++) {
 145                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 146         }
 147
 148         return e;
 149 }
 150
 151 double get_e_pot(t_moldyn *moldyn) {
 152
 153         return(moldyn->potential(moldyn));
 154 }
 155
 156 double get_total_energy(t_moldyn *moldyn) {
 157
 158         double e;
 159
 160         e=get_e_kin(moldyn->atom,moldyn->count);
 161         e+=get_e_pot(moldyn);
 162
 163         return e;
 164 }
 165
 166 t_3dvec get_total_p(t_atom *atom, int count) {
 167
 168         t_3dvec p,p_total;
 169         int i;
 170
 171         v3_zero(&p_total);
 172         for(i=0;i<count;i++) {
 173                 v3_scale(&p,&(atom[i].v),atom[i].mass);
 174                 v3_add(&p_total,&p_total,&p);
 175         }
 176
 177         return p_total;
 178 }
 179
 180 double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
 181
 182         double tau;
 183
 184         tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
 185         tau*=1.0E-9;
 186         if(tau<moldyn->tau)
 187                 printf("[moldyn] warning: time step  (%f > %.15f)\n",
 188                        moldyn->tau,tau);
 189
 190         return tau;
 191 }
 192
 193
 194 /*
 195  *
 196  * 'integration of newtons equation' - algorithms
 197  *
 198  */
 199
 200 /* start the integration */
 201
 202 int moldyn_integrate(t_moldyn *moldyn) {
 203
 204         int i;
 205         int write;
 206
 207         write=moldyn->write;
 208
 209         /* calculate initial forces */
 210         moldyn->force(moldyn);
 211
 212         for(i=0;i<moldyn->time_steps;i++) {
 213                 /* integration step */
 214                 moldyn->integrate(moldyn);
 215
 216                 /* check for visualiziation */
 217                 if(!(i%write)) {
 218                         visual_atoms(moldyn->visual,i*moldyn->tau,
 219                                      moldyn->atom,moldyn->count);
 220                         printf("finished %d / %d\n",i,moldyn->time_steps);
 221                 }
 222         }
 223
 224         return 0;
 225 }
 226
 227 /* velocity verlet */
 228
 229 int velocity_verlet(t_moldyn *moldyn) {
 230
 231         int i,count;
 232         double tau,tau_square;
 233         t_3dvec delta;
 234         t_atom *atom;
 235
 236         atom=moldyn->atom;
 237         count=moldyn->count;
 238         tau=moldyn->tau;
 239
 240         tau_square=tau*tau;
 241
 242         for(i=0;i<count;i++) {
 243                 /* new positions */
 244                 v3_scale(&delta,&(atom[i].v),tau);
 245                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 246                 v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
 247                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 248                 v3_per_bound(&(atom[i].r),&(moldyn->dim));
 249
 250                 /* velocities */
 251                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 252                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 253         }
 254
 255         /* forces depending on chosen potential */
 256         moldyn->force(moldyn);
 257
 258         for(i=0;i<count;i++) {
 259                 /* again velocities */
 260                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 261                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 262         }
 263
 264         return 0;
 265 }
 266
 267
 268 /*
 269  *
 270  * potentials & corresponding forces
 271  *
 272  */
 273
 274 /* harmonic oscillator potential and force */
 275
 276 double potential_harmonic_oscillator(t_moldyn *moldyn) {
 277
 278         t_ho_params *params;
 279         t_atom *atom;
 280         int i,j;
 281         int count;
 282         t_3dvec distance;
 283         double d,u;
 284         double sc,equi_dist;
 285
 286         params=moldyn->pot_params;
 287         atom=moldyn->atom;
 288         sc=params->spring_constant;
 289         equi_dist=params->equilibrium_distance;
 290         count=moldyn->count;
 291
 292         u=0.0;
 293         for(i=0;i<count;i++) {
 294                 for(j=0;j<i;j++) {
 295                         v3_sub(&distance,&(atom[i].r),&(atom[j].r));
 296                         d=v3_norm(&distance);
 297                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 298                 }
 299         }
 300
 301         return u;
 302 }
 303
 304 int force_harmonic_oscillator(t_moldyn *moldyn) {
 305
 306         t_ho_params *params;
 307         int i,j,count;
 308         t_atom *atom;
 309         t_3dvec distance;
 310         t_3dvec force;
 311         double d;
 312         double sc,equi_dist;
 313
 314         atom=moldyn->atom;
 315         count=moldyn->count;
 316         params=moldyn->pot_params;
 317         sc=params->spring_constant;
 318         equi_dist=params->equilibrium_distance;
 319
 320         for(i=0;i<count;i++) v3_zero(&(atom[i].f));
 321
 322         for(i=0;i<count;i++) {
 323                 for(j=0;j<i;j++) {
 324                         v3_sub(&distance,&(atom[i].r),&(atom[j].r));
 325                         v3_per_bound(&distance,&(moldyn->dim));
 326                         d=v3_norm(&distance);
 327                         if(d<=moldyn->cutoff) {
 328                                 v3_scale(&force,&distance,
 329                                          (-sc*(1.0-(equi_dist/d))));
 330                                 v3_add(&(atom[i].f),&(atom[i].f),&force);
 331                                 v3_sub(&(atom[j].f),&(atom[j].f),&force);
 332                         }
 333                 }
 334         }
 335
 336         return 0;
 337 }
 338
 339
 340 /* lennard jones potential & force for one sort of atoms */
 341
 342 double potential_lennard_jones(t_moldyn *moldyn) {
 343
 344         t_lj_params *params;
 345         t_atom *atom;
 346         int i,j;
 347         int count;
 348         t_3dvec distance;
 349         double d,help;
 350         double u;
 351         double eps,sig6,sig12;
 352
 353         params=moldyn->pot_params;
 354         atom=moldyn->atom;
 355         count=moldyn->count;
 356         eps=params->epsilon;
 357         sig6=params->sigma6;
 358         sig12=params->sigma12;
 359
 360         u=0.0;
 361         for(i=0;i<count;i++) {
 362                 for(j=0;j<i;j++) {
 363                         v3_sub(&distance,&(atom[j].r),&(atom[i].r));
 364                         d=1.0/v3_absolute_square(&distance);    /* 1/r^2 */
 365                         help=d*d;                               /* 1/r^4 */
 366                         help*=d;                                /* 1/r^6 */
 367                         d=help*help;                            /* 1/r^12 */
 368                         u+=eps*(sig12*d-sig6*help);
 369                 }
 370         }
 371
 372         return u;
 373 }
 374
 375 int force_lennard_jones(t_moldyn *moldyn) {
 376
 377         t_lj_params *params;
 378         int i,j,count;
 379         t_atom *atom;
 380         t_3dvec distance;
 381         t_3dvec force;
 382         double d,h1,h2;
 383         double eps,sig6,sig12;
 384
 385         atom=moldyn->atom;
 386         count=moldyn->count;
 387         params=moldyn->pot_params;
 388         eps=params->epsilon;
 389         sig6=params->sigma6;
 390         sig12=params->sigma12;
 391
 392         for(i=0;i<count;i++) v3_zero(&(atom[i].f));
 393
 394         for(i=0;i<count;i++) {
 395                 for(j=0;j<i;j++) {
 396                         v3_sub(&distance,&(atom[j].r),&(atom[i].r));
 397                         v3_per_bound(&distance,&(moldyn->dim));
 398                         d=v3_absolute_square(&distance);
 399                         if(d<=moldyn->cutoff_square) {
 400                                 h1=1.0/d;                       /* 1/r^2 */
 401                                 d=h1*h1;                        /* 1/r^4 */
 402                                 h2=d*d;                         /* 1/r^8 */
 403                                 h1*=d;                          /* 1/r^6 */
 404                                 h1*=h2;                         /* 1/r^14 */
 405                                 h1*=sig12;
 406                                 h2*=sig6;
 407                                 d=12.0*h1-6.0*h2;
 408                                 d*=eps;
 409                                 v3_scale(&force,&distance,d);
 410                                 v3_add(&(atom[j].f),&(atom[j].f),&force);
 411                                 v3_sub(&(atom[i].f),&(atom[i].f),&force);
 412                         }
 413                 }
 414         }
 415
 416         return 0;
 417 }
 418