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 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <sys/types.h>
  13 #include <sys/stat.h>
  14 #include <fcntl.h>
  15 #include <unistd.h>
  16 #include <math.h>
  17
  18 #include "moldyn.h"
  19
  20 #include "math/math.h"
  21 #include "init/init.h"
  22 #include "random/random.h"
  23 #include "visual/visual.h"
  24
  25 int moldyn_usage(char **argv) {
  26
  27         printf("\n%s usage:\n\n",argv[0]);
  28         printf("--- general options ---\n");
  29         printf("-E <steps> <file> (log total energy)\n");
  30         printf("-M <steps> <file> (log total momentum)\n");
  31         printf("-D <steps> <file> (dump total information)\n");
  32         printf("-S <steps> <filebase> (single save file)\n");
  33         printf("--- physics options ---\n");
  34         printf("-T <temperature> [K] (%f)\n",MOLDYN_TEMP);
  35         printf("-t <timestep tau> [s] (%f)\n",MOLDYN_TAU);
  36         printf("-R <runs> (%d)\n",MOLDYN_RUNS);
  37         printf("\n");
  38
  39         return 0;
  40 }
  41
  42 int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
  43
  44         int i;
  45
  46         memset(moldyn,0,sizeof(t_moldyn));
  47
  48         /* default values */
  49         moldyn->t=MOLDYN_TEMP;
  50         moldyn->tau=MOLDYN_TAU;
  51         moldyn->time_steps=MOLDYN_RUNS;
  52
  53         /* parse argv */
  54         for(i=1;i<argc;i++) {
  55                 if(argv[i][0]=='-') {
  56                         switch(argv[i][1]){
  57                                 case 'E':
  58                                         moldyn->ewrite=atoi(argv[++i]);
  59                                         strncpy(moldyn->efb,argv[++i],64);
  60                                         break;
  61                                 case 'M':
  62                                         moldyn->mwrite=atoi(argv[++i]);
  63                                         strncpy(moldyn->mfb,argv[++i],64);
  64                                         break;
  65                                 case 'D':
  66                                         moldyn->dwrite=atoi(argv[++i]);
  67                                         strncpy(moldyn->dfb,argv[++i],64);
  68                                         break;
  69                                 case 'S':
  70                                         moldyn->swrite=atoi(argv[++i]);
  71                                         strncpy(moldyn->sfb,argv[++i],64);
  72                                         break;
  73                                 case 'T':
  74                                         break;
  75                                         moldyn->t=atof(argv[++i]);
  76                                 case 't':
  77                                         moldyn->tau=atof(argv[++i]);
  78                                         break;
  79                                 case 'R':
  80                                         moldyn->time_steps=atoi(argv[++i]);
  81                                         break;
  82                                 default:
  83                                         printf("unknown option %s\n",argv[i]);
  84                                         moldyn_usage(argv);
  85                                         return -1;
  86                         }
  87                 } else {
  88                         moldyn_usage(argv);
  89                         return -1;
  90                 }
  91         }
  92
  93         return 0;
  94 }
  95
  96 int moldyn_log_init(t_moldyn *moldyn) {
  97
  98         moldyn->lvstat=0;
  99
 100         if(moldyn->ewrite) {
 101                 moldyn->efd=open(moldyn->efb,O_WRONLY|O_CREAT|O_TRUNC);
 102                 if(moldyn->efd<0) {
 103                         perror("[moldyn] efd open");
 104                         return moldyn->efd;
 105                 }
 106                 dprintf(moldyn->efd,"# moldyn total energy logfile\n");
 107                 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_E;
 108         }
 109
 110         if(moldyn->mwrite) {
 111                 moldyn->mfd=open(moldyn->mfb,O_WRONLY|O_CREAT|O_TRUNC);
 112                 if(moldyn->mfd<0) {
 113                         perror("[moldyn] mfd open");
 114                         return moldyn->mfd;
 115                 }
 116                 dprintf(moldyn->mfd,"# moldyn total momentum logfile\n");
 117                 moldyn->lvstat|=MOLDYN_LVSTAT_TOTAL_M;
 118         }
 119
 120         if(moldyn->swrite)
 121                 moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
 122
 123         if(moldyn->dwrite) {
 124                 moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC);
 125                 if(moldyn->dfd<0) {
 126                         perror("[moldyn] dfd open");
 127                         return moldyn->dfd;
 128                 }
 129                 write(moldyn->dfd,moldyn,sizeof(t_moldyn));
 130                 moldyn->lvstat|=MOLDYN_LVSTAT_DUMP;
 131         }
 132
 133         if(moldyn->dwrite)
 134                 moldyn->lvstat|=MOLDYN_LVSTAT_VISUAL;
 135
 136         return 0;
 137 }
 138
 139 int create_lattice(unsigned char type,int element,double mass,double lc,
 140                    int a,int b,int c,t_atom **atom) {
 141
 142         int count;
 143         int ret;
 144         t_3dvec origin;
 145
 146         count=a*b*c;
 147
 148         if(type==FCC) count*=4;
 149         if(type==DIAMOND) count*=8;
 150
 151         *atom=malloc(count*sizeof(t_atom));
 152         if(*atom==NULL) {
 153                 perror("malloc (atoms)");
 154                 return -1;
 155         }
 156
 157         v3_zero(&origin);
 158
 159         switch(type) {
 160                 case FCC:
 161                         ret=fcc_init(a,b,c,lc,*atom,&origin);
 162                         break;
 163                 case DIAMOND:
 164                         ret=diamond_init(a,b,c,lc,*atom,&origin);
 165                         break;
 166                 default:
 167                         printf("unknown lattice type (%02x)\n",type);
 168                         return -1;
 169         }
 170
 171         /* debug */
 172         if(ret!=count) {
 173                 printf("ok, there is something wrong ...\n");
 174                 printf("calculated -> %d atoms\n",count);
 175                 printf("created -> %d atoms\n",ret);
 176                 return -1;
 177         }
 178
 179         while(count) {
 180                 (*atom)[count-1].element=element;
 181                 (*atom)[count-1].mass=mass;
 182                 count-=1;
 183         }
 184
 185         return ret;
 186 }
 187
 188 int destroy_lattice(t_atom *atom) {
 189
 190         if(atom) free(atom);
 191
 192         return 0;
 193 }
 194
 195 int thermal_init(t_moldyn *moldyn,t_random *random,int count) {
 196
 197         /*
 198          * - gaussian distribution of velocities
 199          * - zero total momentum
 200          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 201          */
 202
 203         int i;
 204         double v,sigma;
 205         t_3dvec p_total,delta;
 206         t_atom *atom;
 207
 208         atom=moldyn->atom;
 209
 210         /* gaussian distribution of velocities */
 211         v3_zero(&p_total);
 212         for(i=0;i<count;i++) {
 213                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t/atom[i].mass);
 214                 /* x direction */
 215                 v=sigma*rand_get_gauss(random);
 216                 atom[i].v.x=v;
 217                 p_total.x+=atom[i].mass*v;
 218                 /* y direction */
 219                 v=sigma*rand_get_gauss(random);
 220                 atom[i].v.y=v;
 221                 p_total.y+=atom[i].mass*v;
 222                 /* z direction */
 223                 v=sigma*rand_get_gauss(random);
 224                 atom[i].v.z=v;
 225                 p_total.z+=atom[i].mass*v;
 226         }
 227
 228         /* zero total momentum */
 229         v3_scale(&p_total,&p_total,1.0/count);
 230         for(i=0;i<count;i++) {
 231                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 232                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 233         }
 234
 235         /* velocity scaling */
 236         scale_velocity(moldyn,count);
 237
 238         return 0;
 239 }
 240
 241 int scale_velocity(t_moldyn *moldyn,int count) {
 242
 243         int i;
 244         double e,c;
 245         t_atom *atom;
 246
 247         atom=moldyn->atom;
 248
 249         /*
 250          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 251          */
 252         e=0.0;
 253         for(i=0;i<count;i++)
 254                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 255         c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*moldyn->t));
 256         for(i=0;i<count;i++)
 257                 v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
 258
 259         return 0;
 260 }
 261
 262 double get_e_kin(t_atom *atom,int count) {
 263
 264         int i;
 265         double e;
 266
 267         e=0.0;
 268
 269         for(i=0;i<count;i++) {
 270                 e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 271         }
 272
 273         return e;
 274 }
 275
 276 double get_e_pot(t_moldyn *moldyn) {
 277
 278         return(moldyn->potential(moldyn));
 279 }
 280
 281 double get_total_energy(t_moldyn *moldyn) {
 282
 283         double e;
 284
 285         e=get_e_kin(moldyn->atom,moldyn->count);
 286         e+=get_e_pot(moldyn);
 287
 288         return e;
 289 }
 290
 291 t_3dvec get_total_p(t_atom *atom, int count) {
 292
 293         t_3dvec p,p_total;
 294         int i;
 295
 296         v3_zero(&p_total);
 297         for(i=0;i<count;i++) {
 298                 v3_scale(&p,&(atom[i].v),atom[i].mass);
 299                 v3_add(&p_total,&p_total,&p);
 300         }
 301
 302         return p_total;
 303 }
 304
 305 double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t) {
 306
 307         double tau;
 308
 309         tau=0.05*nn_dist/(sqrt(3.0*K_BOLTZMANN*t/moldyn->atom[0].mass));
 310         tau*=1.0E-9;
 311         if(tau<moldyn->tau)
 312                 printf("[moldyn] warning: time step  (%f > %.15f)\n",
 313                        moldyn->tau,tau);
 314
 315         return tau;
 316 }
 317
 318
 319 /*
 320  *
 321  * 'integration of newtons equation' - algorithms
 322  *
 323  */
 324
 325 /* start the integration */
 326
 327 int moldyn_integrate(t_moldyn *moldyn) {
 328
 329         int i;
 330         unsigned int e,m,s,d,v;
 331         unsigned char lvstat;
 332         t_3dvec p;
 333
 334         int fd;
 335         char fb[128];
 336
 337         e=moldyn->ewrite;
 338         m=moldyn->mwrite;
 339         s=moldyn->swrite;
 340         d=moldyn->dwrite;
 341         v=moldyn->vwrite;
 342
 343         if(!(lvstat&MOLDYN_LVSTAT_INITIALIZED)) {
 344                 printf("[moldyn] warning, lv system not initialized\n");
 345                 return -1;
 346         }
 347
 348         /* calculate initial forces */
 349         moldyn->force(moldyn);
 350
 351         for(i=0;i<moldyn->time_steps;i++) {
 352                 /* integration step */
 353                 moldyn->integrate(moldyn);
 354
 355                 /* check for log & visualiziation */
 356                 if(e) {
 357                         if(!(i%e))
 358                                 dprintf(moldyn->efd,
 359                                         "%.15f %.45f\n",i*moldyn->tau,
 360                                         get_total_energy(moldyn));
 361                 }
 362                 if(m) {
 363                         if(!(i%m)) {
 364                                 p=get_total_p(moldyn->atom,moldyn->count);
 365                                 dprintf(moldyn->mfd,
 366                                         "%.15f %.45f\n",i*moldyn->tau,
 367                                         v3_norm(&p));
 368                         }
 369                 }
 370                 if(s) {
 371                         if(!(i%s)) {
 372                                 snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
 373                                          moldyn->t,i*moldyn->tau);
 374                                 fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
 375                                 if(fd<0) perror("[moldyn] save fd open");
 376                                 else {
 377                                         write(fd,moldyn,sizeof(t_moldyn));
 378                                         write(fd,moldyn->atom,
 379                                               moldyn->count*sizeof(t_atom));
 380                                 }
 381                         }
 382                 }
 383                 if(d) {
 384                         if(!(i%d))
 385                                 write(moldyn->dfd,moldyn->atom,
 386                                       moldyn->count*sizeof(t_atom));
 387
 388                 }
 389                 if(v) {
 390                         if(!(i%v))
 391                                 visual_atoms(moldyn->visual,i*moldyn->tau,
 392                                              moldyn->atom,moldyn->count);
 393                 }
 394         }
 395
 396         return 0;
 397 }
 398
 399 /* velocity verlet */
 400
 401 int velocity_verlet(t_moldyn *moldyn) {
 402
 403         int i,count;
 404         double tau,tau_square;
 405         t_3dvec delta;
 406         t_atom *atom;
 407
 408         atom=moldyn->atom;
 409         count=moldyn->count;
 410         tau=moldyn->tau;
 411
 412         tau_square=tau*tau;
 413
 414         for(i=0;i<count;i++) {
 415                 /* new positions */
 416                 v3_scale(&delta,&(atom[i].v),tau);
 417                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 418                 v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
 419                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 420                 v3_per_bound(&(atom[i].r),&(moldyn->dim));
 421
 422                 /* velocities */
 423                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 424                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 425         }
 426
 427         /* forces depending on chosen potential */
 428         moldyn->force(moldyn);
 429
 430         for(i=0;i<count;i++) {
 431                 /* again velocities */
 432                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 433                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 434         }
 435
 436         return 0;
 437 }
 438
 439
 440 /*
 441  *
 442  * potentials & corresponding forces
 443  *
 444  */
 445
 446 /* harmonic oscillator potential and force */
 447
 448 double potential_harmonic_oscillator(t_moldyn *moldyn) {
 449
 450         t_ho_params *params;
 451         t_atom *atom;
 452         int i,j;
 453         int count;
 454         t_3dvec distance;
 455         double d,u;
 456         double sc,equi_dist;
 457
 458         params=moldyn->pot_params;
 459         atom=moldyn->atom;
 460         sc=params->spring_constant;
 461         equi_dist=params->equilibrium_distance;
 462         count=moldyn->count;
 463
 464         u=0.0;
 465         for(i=0;i<count;i++) {
 466                 for(j=0;j<i;j++) {
 467                         v3_sub(&distance,&(atom[i].r),&(atom[j].r));
 468                         d=v3_norm(&distance);
 469                         u+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
 470                 }
 471         }
 472
 473         return u;
 474 }
 475
 476 int force_harmonic_oscillator(t_moldyn *moldyn) {
 477
 478         t_ho_params *params;
 479         int i,j,count;
 480         t_atom *atom;
 481         t_3dvec distance;
 482         t_3dvec force;
 483         double d;
 484         double sc,equi_dist;
 485
 486         atom=moldyn->atom;
 487         count=moldyn->count;
 488         params=moldyn->pot_params;
 489         sc=params->spring_constant;
 490         equi_dist=params->equilibrium_distance;
 491
 492         for(i=0;i<count;i++) v3_zero(&(atom[i].f));
 493
 494         for(i=0;i<count;i++) {
 495                 for(j=0;j<i;j++) {
 496                         v3_sub(&distance,&(atom[i].r),&(atom[j].r));
 497                         v3_per_bound(&distance,&(moldyn->dim));
 498                         d=v3_norm(&distance);
 499                         if(d<=moldyn->cutoff) {
 500                                 v3_scale(&force,&distance,
 501                                          (-sc*(1.0-(equi_dist/d))));
 502                                 v3_add(&(atom[i].f),&(atom[i].f),&force);
 503                                 v3_sub(&(atom[j].f),&(atom[j].f),&force);
 504                         }
 505                 }
 506         }
 507
 508         return 0;
 509 }
 510
 511
 512 /* lennard jones potential & force for one sort of atoms */
 513
 514 double potential_lennard_jones(t_moldyn *moldyn) {
 515
 516         t_lj_params *params;
 517         t_atom *atom;
 518         int i,j;
 519         int count;
 520         t_3dvec distance;
 521         double d,help;
 522         double u;
 523         double eps,sig6,sig12;
 524
 525         params=moldyn->pot_params;
 526         atom=moldyn->atom;
 527         count=moldyn->count;
 528         eps=params->epsilon;
 529         sig6=params->sigma6;
 530         sig12=params->sigma12;
 531
 532         u=0.0;
 533         for(i=0;i<count;i++) {
 534                 for(j=0;j<i;j++) {
 535                         v3_sub(&distance,&(atom[j].r),&(atom[i].r));
 536                         d=1.0/v3_absolute_square(&distance);    /* 1/r^2 */
 537                         help=d*d;                               /* 1/r^4 */
 538                         help*=d;                                /* 1/r^6 */
 539                         d=help*help;                            /* 1/r^12 */
 540                         u+=eps*(sig12*d-sig6*help);
 541                 }
 542         }
 543
 544         return u;
 545 }
 546
 547 int force_lennard_jones(t_moldyn *moldyn) {
 548
 549         t_lj_params *params;
 550         int i,j,count;
 551         t_atom *atom;
 552         t_3dvec distance;
 553         t_3dvec force;
 554         double d,h1,h2;
 555         double eps,sig6,sig12;
 556
 557         atom=moldyn->atom;
 558         count=moldyn->count;
 559         params=moldyn->pot_params;
 560         eps=params->epsilon;
 561         sig6=params->sigma6;
 562         sig12=params->sigma12;
 563
 564         for(i=0;i<count;i++) v3_zero(&(atom[i].f));
 565
 566         for(i=0;i<count;i++) {
 567                 for(j=0;j<i;j++) {
 568                         v3_sub(&distance,&(atom[j].r),&(atom[i].r));
 569                         v3_per_bound(&distance,&(moldyn->dim));
 570                         d=v3_absolute_square(&distance);
 571                         if(d<=moldyn->cutoff_square) {
 572                                 h1=1.0/d;                       /* 1/r^2 */
 573                                 d=h1*h1;                        /* 1/r^4 */
 574                                 h2=d*d;                         /* 1/r^8 */
 575                                 h1*=d;                          /* 1/r^6 */
 576                                 h1*=h2;                         /* 1/r^14 */
 577                                 h1*=sig12;
 578                                 h2*=sig6;
 579                                 d=12.0*h1-6.0*h2;
 580                                 d*=eps;
 581                                 v3_scale(&force,&distance,d);
 582                                 v3_add(&(atom[j].f),&(atom[j].f),&force);
 583                                 v3_sub(&(atom[i].f),&(atom[i].f),&force);
 584                         }
 585                 }
 586         }
 587
 588         return 0;
 589 }
 590