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