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 <sys/time.h>
  17 #include <time.h>
  18 #include <math.h>
  19
  20 #include "moldyn.h"
  21 #include "report/report.h"
  22
  23 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
  24
  25         printf("[moldyn] init\n");
  26
  27         memset(moldyn,0,sizeof(t_moldyn));
  28
  29         moldyn->argc=argc;
  30         moldyn->args=argv;
  31
  32         rand_init(&(moldyn->random),NULL,1);
  33         moldyn->random.status|=RAND_STAT_VERBOSE;
  34
  35         return 0;
  36 }
  37
  38 int moldyn_shutdown(t_moldyn *moldyn) {
  39
  40         printf("[moldyn] shutdown\n");
  41
  42         moldyn_log_shutdown(moldyn);
  43         link_cell_shutdown(moldyn);
  44         rand_close(&(moldyn->random));
  45         free(moldyn->atom);
  46
  47         return 0;
  48 }
  49
  50 int set_int_alg(t_moldyn *moldyn,u8 algo) {
  51
  52         printf("[moldyn] integration algorithm: ");
  53
  54         switch(algo) {
  55                 case MOLDYN_INTEGRATE_VERLET:
  56                         moldyn->integrate=velocity_verlet;
  57                         printf("velocity verlet\n");
  58                         break;
  59                 default:
  60                         printf("unknown integration algorithm: %02x\n",algo);
  61                         printf("unknown\n");
  62                         return -1;
  63         }
  64
  65         return 0;
  66 }
  67
  68 int set_cutoff(t_moldyn *moldyn,double cutoff) {
  69
  70         moldyn->cutoff=cutoff;
  71
  72         printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff);
  73
  74         return 0;
  75 }
  76
  77 int set_temperature(t_moldyn *moldyn,double t_ref) {
  78
  79         moldyn->t_ref=t_ref;
  80
  81         printf("[moldyn] temperature [K]: %f\n",moldyn->t_ref);
  82
  83         return 0;
  84 }
  85
  86 int set_pressure(t_moldyn *moldyn,double p_ref) {
  87
  88         moldyn->p_ref=p_ref;
  89
  90         printf("[moldyn] pressure [bar]: %f\n",moldyn->p_ref/BAR);
  91
  92         return 0;
  93 }
  94
  95 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
  96
  97         moldyn->pt_scale=(ptype|ttype);
  98         moldyn->t_tc=ttc;
  99         moldyn->p_tc=ptc;
 100
 101         printf("[moldyn] p/t scaling:\n");
 102
 103         printf("  p: %s",ptype?"yes":"no ");
 104         if(ptype)
 105                 printf(" | type: %02x | factor: %f",ptype,ptc);
 106         printf("\n");
 107
 108         printf("  t: %s",ttype?"yes":"no ");
 109         if(ttype)
 110                 printf(" | type: %02x | factor: %f",ttype,ttc);
 111         printf("\n");
 112
 113         return 0;
 114 }
 115
 116 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
 117
 118         moldyn->dim.x=x;
 119         moldyn->dim.y=y;
 120         moldyn->dim.z=z;
 121
 122         moldyn->volume=x*y*z;
 123
 124         if(visualize) {
 125                 moldyn->vis.dim.x=x;
 126                 moldyn->vis.dim.y=y;
 127                 moldyn->vis.dim.z=z;
 128         }
 129
 130         moldyn->dv=0.000001*moldyn->volume;
 131
 132         printf("[moldyn] dimensions in A and A^3 respectively:\n");
 133         printf("  x: %f\n",moldyn->dim.x);
 134         printf("  y: %f\n",moldyn->dim.y);
 135         printf("  z: %f\n",moldyn->dim.z);
 136         printf("  volume: %f\n",moldyn->volume);
 137         printf("  visualize simulation box: %s\n",visualize?"yes":"no");
 138         printf("  delta volume (pressure calc): %f\n",moldyn->dv);
 139
 140         return 0;
 141 }
 142
 143 int set_nn_dist(t_moldyn *moldyn,double dist) {
 144
 145         moldyn->nnd=dist;
 146
 147         return 0;
 148 }
 149
 150 int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
 151
 152         printf("[moldyn] periodic boundary conditions:\n");
 153
 154         if(x)
 155                 moldyn->status|=MOLDYN_STAT_PBX;
 156
 157         if(y)
 158                 moldyn->status|=MOLDYN_STAT_PBY;
 159
 160         if(z)
 161                 moldyn->status|=MOLDYN_STAT_PBZ;
 162
 163         printf("  x: %s\n",x?"yes":"no");
 164         printf("  y: %s\n",y?"yes":"no");
 165         printf("  z: %s\n",z?"yes":"no");
 166
 167         return 0;
 168 }
 169
 170 int set_potential1b(t_moldyn *moldyn,pf_func1b func) {
 171
 172         moldyn->func1b=func;
 173
 174         return 0;
 175 }
 176
 177 int set_potential2b(t_moldyn *moldyn,pf_func2b func) {
 178
 179         moldyn->func2b=func;
 180
 181         return 0;
 182 }
 183
 184 int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) {
 185
 186         moldyn->func3b_j1=func;
 187
 188         return 0;
 189 }
 190
 191 int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) {
 192
 193         moldyn->func3b_j2=func;
 194
 195         return 0;
 196 }
 197
 198 int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) {
 199
 200         moldyn->func3b_j3=func;
 201
 202         return 0;
 203 }
 204
 205 int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) {
 206
 207         moldyn->func3b_k1=func;
 208
 209         return 0;
 210 }
 211
 212 int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) {
 213
 214         moldyn->func3b_k2=func;
 215
 216         return 0;
 217 }
 218
 219 int set_potential_params(t_moldyn *moldyn,void *params) {
 220
 221         moldyn->pot_params=params;
 222
 223         return 0;
 224 }
 225
 226 int set_avg_skip(t_moldyn *moldyn,int skip) {
 227
 228         printf("[moldyn] skip %d steps before starting average calc\n",skip);
 229         moldyn->avg_skip=skip;
 230
 231         return 0;
 232 }
 233
 234 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
 235
 236         strncpy(moldyn->vlsdir,dir,127);
 237
 238         return 0;
 239 }
 240
 241 int moldyn_set_report(t_moldyn *moldyn,char *author,char *title) {
 242
 243         strncpy(moldyn->rauthor,author,63);
 244         strncpy(moldyn->rtitle,title,63);
 245
 246         return 0;
 247 }
 248
 249 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 250
 251         char filename[128];
 252         int ret;
 253
 254         printf("[moldyn] set log: ");
 255
 256         switch(type) {
 257                 case LOG_TOTAL_ENERGY:
 258                         moldyn->ewrite=timer;
 259                         snprintf(filename,127,"%s/energy",moldyn->vlsdir);
 260                         moldyn->efd=open(filename,
 261                                          O_WRONLY|O_CREAT|O_EXCL,
 262                                          S_IRUSR|S_IWUSR);
 263                         if(moldyn->efd<0) {
 264                                 perror("[moldyn] energy log fd open");
 265                                 return moldyn->efd;
 266                         }
 267                         dprintf(moldyn->efd,"# total energy log file\n");
 268                         printf("total energy (%d)\n",timer);
 269                         break;
 270                 case LOG_TOTAL_MOMENTUM:
 271                         moldyn->mwrite=timer;
 272                         snprintf(filename,127,"%s/momentum",moldyn->vlsdir);
 273                         moldyn->mfd=open(filename,
 274                                          O_WRONLY|O_CREAT|O_EXCL,
 275                                          S_IRUSR|S_IWUSR);
 276                         if(moldyn->mfd<0) {
 277                                 perror("[moldyn] momentum log fd open");
 278                                 return moldyn->mfd;
 279                         }
 280                         dprintf(moldyn->efd,"# total momentum log file\n");
 281                         printf("total momentum (%d)\n",timer);
 282                         break;
 283                 case LOG_PRESSURE:
 284                         moldyn->pwrite=timer;
 285                         snprintf(filename,127,"%s/pressure",moldyn->vlsdir);
 286                         moldyn->pfd=open(filename,
 287                                          O_WRONLY|O_CREAT|O_EXCL,
 288                                          S_IRUSR|S_IWUSR);
 289                         if(moldyn->pfd<0) {
 290                                 perror("[moldyn] pressure log file\n");
 291                                 return moldyn->pfd;
 292                         }
 293                         dprintf(moldyn->pfd,"# pressure log file\n");
 294                         printf("pressure (%d)\n",timer);
 295                         break;
 296                 case LOG_TEMPERATURE:
 297                         moldyn->twrite=timer;
 298                         snprintf(filename,127,"%s/temperature",moldyn->vlsdir);
 299                         moldyn->tfd=open(filename,
 300                                          O_WRONLY|O_CREAT|O_EXCL,
 301                                          S_IRUSR|S_IWUSR);
 302                         if(moldyn->tfd<0) {
 303                                 perror("[moldyn] temperature log file\n");
 304                                 return moldyn->tfd;
 305                         }
 306                         dprintf(moldyn->tfd,"# temperature log file\n");
 307                         printf("temperature (%d)\n",timer);
 308                         break;
 309                 case SAVE_STEP:
 310                         moldyn->swrite=timer;
 311                         printf("save file (%d)\n",timer);
 312                         break;
 313                 case VISUAL_STEP:
 314                         moldyn->vwrite=timer;
 315                         ret=visual_init(&(moldyn->vis),moldyn->vlsdir);
 316                         if(ret<0) {
 317                                 printf("[moldyn] visual init failure\n");
 318                                 return ret;
 319                         }
 320                         printf("visual file (%d)\n",timer);
 321                         break;
 322                 case CREATE_REPORT:
 323                         snprintf(filename,127,"%s/report.tex",moldyn->vlsdir);
 324                         moldyn->rfd=open(filename,
 325                                          O_WRONLY|O_CREAT|O_EXCL,
 326                                          S_IRUSR|S_IWUSR);
 327                         if(moldyn->rfd<0) {
 328                                 perror("[moldyn] report fd open");
 329                                 return moldyn->rfd;
 330                         }
 331                         printf("report -> ");
 332                         if(moldyn->efd) {
 333                                 snprintf(filename,127,"%s/e_plot.scr",
 334                                          moldyn->vlsdir);
 335                                 moldyn->epfd=open(filename,
 336                                                  O_WRONLY|O_CREAT|O_EXCL,
 337                                                  S_IRUSR|S_IWUSR);
 338                                 if(moldyn->epfd<0) {
 339                                         perror("[moldyn] energy plot fd open");
 340                                         return moldyn->epfd;
 341                                 }
 342                                 dprintf(moldyn->epfd,e_plot_script);
 343                                 close(moldyn->epfd);
 344                                 printf("energy ");
 345                         }
 346                         if(moldyn->pfd) {
 347                                 snprintf(filename,127,"%s/pressure_plot.scr",
 348                                          moldyn->vlsdir);
 349                                 moldyn->ppfd=open(filename,
 350                                                   O_WRONLY|O_CREAT|O_EXCL,
 351                                                   S_IRUSR|S_IWUSR);
 352                                 if(moldyn->ppfd<0) {
 353                                         perror("[moldyn] p plot fd open");
 354                                         return moldyn->ppfd;
 355                                 }
 356                                 dprintf(moldyn->ppfd,pressure_plot_script);
 357                                 close(moldyn->ppfd);
 358                                 printf("pressure ");
 359                         }
 360                         if(moldyn->tfd) {
 361                                 snprintf(filename,127,"%s/temperature_plot.scr",
 362                                          moldyn->vlsdir);
 363                                 moldyn->tpfd=open(filename,
 364                                                   O_WRONLY|O_CREAT|O_EXCL,
 365                                                   S_IRUSR|S_IWUSR);
 366                                 if(moldyn->tpfd<0) {
 367                                         perror("[moldyn] t plot fd open");
 368                                         return moldyn->tpfd;
 369                                 }
 370                                 dprintf(moldyn->tpfd,temperature_plot_script);
 371                                 close(moldyn->tpfd);
 372                                 printf("temperature ");
 373                         }
 374                         dprintf(moldyn->rfd,report_start,
 375                                 moldyn->rauthor,moldyn->rtitle);
 376                         printf("\n");
 377                         break;
 378                 default:
 379                         printf("unknown log type: %02x\n",type);
 380                         return -1;
 381         }
 382
 383         return 0;
 384 }
 385
 386 int moldyn_log_shutdown(t_moldyn *moldyn) {
 387
 388         char sc[256];
 389
 390         printf("[moldyn] log shutdown\n");
 391         if(moldyn->efd) {
 392                 close(moldyn->efd);
 393                 if(moldyn->rfd) {
 394                         dprintf(moldyn->rfd,report_energy);
 395                         snprintf(sc,255,"cd %s && gnuplot e_plot.scr",
 396                                  moldyn->vlsdir);
 397                         system(sc);
 398                 }
 399         }
 400         if(moldyn->mfd) close(moldyn->mfd);
 401         if(moldyn->pfd) {
 402                 close(moldyn->pfd);
 403                 if(moldyn->rfd)
 404                         dprintf(moldyn->rfd,report_pressure);
 405                         snprintf(sc,255,"cd %s && gnuplot pressure_plot.scr",
 406                                  moldyn->vlsdir);
 407                         system(sc);
 408         }
 409         if(moldyn->tfd) {
 410                 close(moldyn->tfd);
 411                 if(moldyn->rfd)
 412                         dprintf(moldyn->rfd,report_temperature);
 413                         snprintf(sc,255,"cd %s && gnuplot temperature_plot.scr",
 414                                  moldyn->vlsdir);
 415                         system(sc);
 416         }
 417         if(moldyn->rfd) {
 418                 dprintf(moldyn->rfd,report_end);
 419                 close(moldyn->rfd);
 420                 snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
 421                          moldyn->vlsdir);
 422                 system(sc);
 423                 snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
 424                          moldyn->vlsdir);
 425                 system(sc);
 426                 snprintf(sc,255,"cd %s && dvipdf report >/dev/null 2>&1",
 427                          moldyn->vlsdir);
 428                 system(sc);
 429         }
 430         if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
 431
 432         return 0;
 433 }
 434
 435 /*
 436  * creating lattice functions
 437  */
 438
 439 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 440                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) {
 441
 442         int new,count;
 443         int ret;
 444         t_3dvec orig;
 445         void *ptr;
 446         t_atom *atom;
 447
 448         new=a*b*c;
 449         count=moldyn->count;
 450
 451         /* how many atoms do we expect */
 452         if(type==CUBIC) new*=1;
 453         if(type==FCC) new*=4;
 454         if(type==DIAMOND) new*=8;
 455
 456         /* allocate space for atoms */
 457         ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom));
 458         if(!ptr) {
 459                 perror("[moldyn] realloc (create lattice)");
 460                 return -1;
 461         }
 462         moldyn->atom=ptr;
 463         atom=&(moldyn->atom[count]);
 464
 465         /* no atoms on the boundaries (only reason: it looks better!) */
 466         if(!origin) {
 467                 orig.x=0.5*lc;
 468                 orig.y=0.5*lc;
 469                 orig.z=0.5*lc;
 470         }
 471         else {
 472                 orig.x=origin->x;
 473                 orig.y=origin->y;
 474                 orig.z=origin->z;
 475         }
 476
 477         switch(type) {
 478                 case CUBIC:
 479                         set_nn_dist(moldyn,lc);
 480                         ret=cubic_init(a,b,c,lc,atom,&orig);
 481                         break;
 482                 case FCC:
 483                         if(!origin)
 484                                 v3_scale(&orig,&orig,0.5);
 485                         set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
 486                         ret=fcc_init(a,b,c,lc,atom,&orig);
 487                         break;
 488                 case DIAMOND:
 489                         if(!origin)
 490                                 v3_scale(&orig,&orig,0.25);
 491                         set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
 492                         ret=diamond_init(a,b,c,lc,atom,&orig);
 493                         break;
 494                 default:
 495                         printf("unknown lattice type (%02x)\n",type);
 496                         return -1;
 497         }
 498
 499         /* debug */
 500         if(ret!=new) {
 501                 printf("[moldyn] creating lattice failed\n");
 502                 printf("  amount of atoms\n");
 503                 printf("  - expected: %d\n",new);
 504                 printf("  - created: %d\n",ret);
 505                 return -1;
 506         }
 507
 508         moldyn->count+=new;
 509         printf("[moldyn] created lattice with %d atoms\n",new);
 510
 511         for(ret=0;ret<new;ret++) {
 512                 atom[ret].element=element;
 513                 atom[ret].mass=mass;
 514                 atom[ret].attr=attr;
 515                 atom[ret].brand=brand;
 516                 atom[ret].tag=count+ret;
 517                 check_per_bound(moldyn,&(atom[ret].r));
 518                 atom[ret].r_0=atom[ret].r;
 519         }
 520
 521         /* update total system mass */
 522         total_mass_calc(moldyn);
 523
 524         return ret;
 525 }
 526
 527 int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 528              t_3dvec *r,t_3dvec *v) {
 529
 530         t_atom *atom;
 531         void *ptr;
 532         int count;
 533
 534         atom=moldyn->atom;
 535         count=(moldyn->count)++;
 536
 537         ptr=realloc(atom,(count+1)*sizeof(t_atom));
 538         if(!ptr) {
 539                 perror("[moldyn] realloc (add atom)");
 540                 return -1;
 541         }
 542         moldyn->atom=ptr;
 543
 544         atom=moldyn->atom;
 545         atom[count].r=*r;
 546         atom[count].v=*v;
 547         atom[count].element=element;
 548         atom[count].mass=mass;
 549         atom[count].brand=brand;
 550         atom[count].tag=count;
 551         atom[count].attr=attr;
 552         check_per_bound(moldyn,&(atom[count].r));
 553         atom[count].r_0=atom[count].r;
 554
 555         /* update total system mass */
 556         total_mass_calc(moldyn);
 557
 558         return 0;
 559 }
 560
 561 /* cubic init */
 562 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 563
 564         int count;
 565         t_3dvec r;
 566         int i,j,k;
 567         t_3dvec o;
 568
 569         count=0;
 570         if(origin)
 571                 v3_copy(&o,origin);
 572         else
 573                 v3_zero(&o);
 574
 575         r.x=o.x;
 576         for(i=0;i<a;i++) {
 577                 r.y=o.y;
 578                 for(j=0;j<b;j++) {
 579                         r.z=o.z;
 580                         for(k=0;k<c;k++) {
 581                                 v3_copy(&(atom[count].r),&r);
 582                                 count+=1;
 583                                 r.z+=lc;
 584                         }
 585                         r.y+=lc;
 586                 }
 587                 r.x+=lc;
 588         }
 589
 590         for(i=0;i<count;i++) {
 591                 atom[i].r.x-=(a*lc)/2.0;
 592                 atom[i].r.y-=(b*lc)/2.0;
 593                 atom[i].r.z-=(c*lc)/2.0;
 594         }
 595
 596         return count;
 597 }
 598
 599 /* fcc lattice init */
 600 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 601
 602         int count;
 603         int i,j,k,l;
 604         t_3dvec o,r,n;
 605         t_3dvec basis[3];
 606
 607         count=0;
 608         if(origin)
 609                 v3_copy(&o,origin);
 610         else
 611                 v3_zero(&o);
 612
 613         /* construct the basis */
 614         memset(basis,0,3*sizeof(t_3dvec));
 615         basis[0].x=0.5*lc;
 616         basis[0].y=0.5*lc;
 617         basis[1].x=0.5*lc;
 618         basis[1].z=0.5*lc;
 619         basis[2].y=0.5*lc;
 620         basis[2].z=0.5*lc;
 621
 622         /* fill up the room */
 623         r.x=o.x;
 624         for(i=0;i<a;i++) {
 625                 r.y=o.y;
 626                 for(j=0;j<b;j++) {
 627                         r.z=o.z;
 628                         for(k=0;k<c;k++) {
 629                                 /* first atom */
 630                                 v3_copy(&(atom[count].r),&r);
 631                                 count+=1;
 632                                 r.z+=lc;
 633                                 /* the three face centered atoms */
 634                                 for(l=0;l<3;l++) {
 635                                         v3_add(&n,&r,&basis[l]);
 636                                         v3_copy(&(atom[count].r),&n);
 637                                         count+=1;
 638                                 }
 639                         }
 640                         r.y+=lc;
 641                 }
 642                 r.x+=lc;
 643         }
 644
 645         /* coordinate transformation */
 646         for(i=0;i<count;i++) {
 647                 atom[i].r.x-=(a*lc)/2.0;
 648                 atom[i].r.y-=(b*lc)/2.0;
 649                 atom[i].r.z-=(c*lc)/2.0;
 650         }
 651
 652         return count;
 653 }
 654
 655 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 656
 657         int count;
 658         t_3dvec o;
 659
 660         count=fcc_init(a,b,c,lc,atom,origin);
 661
 662         o.x=0.25*lc;
 663         o.y=0.25*lc;
 664         o.z=0.25*lc;
 665
 666         if(origin) v3_add(&o,&o,origin);
 667
 668         count+=fcc_init(a,b,c,lc,&atom[count],&o);
 669
 670         return count;
 671 }
 672
 673 int destroy_atoms(t_moldyn *moldyn) {
 674
 675         if(moldyn->atom) free(moldyn->atom);
 676
 677         return 0;
 678 }
 679
 680 int thermal_init(t_moldyn *moldyn,u8 equi_init) {
 681
 682         /*
 683          * - gaussian distribution of velocities
 684          * - zero total momentum
 685          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 686          */
 687
 688         int i;
 689         double v,sigma;
 690         t_3dvec p_total,delta;
 691         t_atom *atom;
 692         t_random *random;
 693
 694         atom=moldyn->atom;
 695         random=&(moldyn->random);
 696
 697         printf("[moldyn] thermal init (equi init: %s)\n",equi_init?"yes":"no");
 698
 699         /* gaussian distribution of velocities */
 700         v3_zero(&p_total);
 701         for(i=0;i<moldyn->count;i++) {
 702                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
 703                 /* x direction */
 704                 v=sigma*rand_get_gauss(random);
 705                 atom[i].v.x=v;
 706                 p_total.x+=atom[i].mass*v;
 707                 /* y direction */
 708                 v=sigma*rand_get_gauss(random);
 709                 atom[i].v.y=v;
 710                 p_total.y+=atom[i].mass*v;
 711                 /* z direction */
 712                 v=sigma*rand_get_gauss(random);
 713                 atom[i].v.z=v;
 714                 p_total.z+=atom[i].mass*v;
 715         }
 716
 717         /* zero total momentum */
 718         v3_scale(&p_total,&p_total,1.0/moldyn->count);
 719         for(i=0;i<moldyn->count;i++) {
 720                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 721                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 722         }
 723
 724         /* velocity scaling */
 725         scale_velocity(moldyn,equi_init);
 726
 727         return 0;
 728 }
 729
 730 double total_mass_calc(t_moldyn *moldyn) {
 731
 732         int i;
 733
 734         moldyn->mass=0.0;
 735
 736         for(i=0;i<moldyn->count;i++)
 737                 moldyn->mass+=moldyn->atom[i].mass;
 738
 739         return moldyn->mass;
 740 }
 741
 742 double temperature_calc(t_moldyn *moldyn) {
 743
 744         /* assume up to date kinetic energy, which is 3/2 N k_B T */
 745
 746         moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
 747
 748         return moldyn->t;
 749 }
 750
 751 double get_temperature(t_moldyn *moldyn) {
 752
 753         return moldyn->t;
 754 }
 755
 756 int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 757
 758         int i;
 759         double e,scale;
 760         t_atom *atom;
 761         int count;
 762
 763         atom=moldyn->atom;
 764
 765         /*
 766          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 767          */
 768
 769         /* get kinetic energy / temperature & count involved atoms */
 770         e=0.0;
 771         count=0;
 772         for(i=0;i<moldyn->count;i++) {
 773                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) {
 774                         e+=atom[i].mass*v3_absolute_square(&(atom[i].v));
 775                         count+=1;
 776                 }
 777         }
 778         e*=0.5;
 779         if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN);
 780         else return 0;  /* no atoms involved in scaling! */
 781
 782         /* (temporary) hack for e,t = 0 */
 783         if(e==0.0) {
 784         moldyn->t=0.0;
 785                 if(moldyn->t_ref!=0.0) {
 786                         thermal_init(moldyn,equi_init);
 787                         return 0;
 788                 }
 789                 else
 790                         return 0; /* no scaling needed */
 791         }
 792
 793
 794         /* get scaling factor */
 795         scale=moldyn->t_ref/moldyn->t;
 796         if(equi_init&TRUE)
 797                 scale*=2.0;
 798         else
 799                 if(moldyn->pt_scale&T_SCALE_BERENDSEN)
 800                         scale=1.0+(scale-1.0)/moldyn->t_tc;
 801         scale=sqrt(scale);
 802
 803         /* velocity scaling */
 804         for(i=0;i<moldyn->count;i++) {
 805                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
 806                         v3_scale(&(atom[i].v),&(atom[i].v),scale);
 807         }
 808
 809         return 0;
 810 }
 811
 812 double ideal_gas_law_pressure(t_moldyn *moldyn) {
 813
 814         double p;
 815
 816         p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume;
 817
 818         return p;
 819 }
 820
 821 double virial_sum(t_moldyn *moldyn) {
 822
 823         int i;
 824         double v;
 825         t_virial *virial;
 826
 827         /* virial (sum over atom virials) */
 828         v=0.0;
 829         for(i=0;i<moldyn->count;i++) {
 830                 virial=&(moldyn->atom[i].virial);
 831                 v+=(virial->xx+virial->yy+virial->zz);
 832         }
 833         moldyn->virial=v;
 834
 835         /* global virial (absolute coordinates) */
 836         virial=&(moldyn->gvir);
 837         moldyn->gv=virial->xx+virial->yy+virial->zz;
 838
 839         return moldyn->virial;
 840 }
 841
 842 double pressure_calc(t_moldyn *moldyn) {
 843
 844         /*
 845          * PV = NkT + <W>
 846          * with W = 1/3 sum_i f_i r_i (- skipped!)
 847          * virial = sum_i f_i r_i
 848          *
 849          * => P = (2 Ekin + virial) / (3V)
 850          */
 851
 852         /* assume up to date virial & up to date kinetic energy */
 853
 854         /* pressure (atom virials) */
 855         moldyn->p=2.0*moldyn->ekin+moldyn->virial;
 856         moldyn->p/=(3.0*moldyn->volume);
 857
 858         /* pressure (absolute coordinates) */
 859         moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
 860         moldyn->gp/=(3.0*moldyn->volume);
 861
 862         return moldyn->p;
 863 }
 864
 865 int average_and_fluctuation_calc(t_moldyn *moldyn) {
 866
 867         if(moldyn->total_steps<moldyn->avg_skip)
 868                 return 0;
 869
 870         int denom=moldyn->total_steps+1-moldyn->avg_skip;
 871
 872         /* assume up to date energies, temperature, pressure etc */
 873
 874         /* kinetic energy */
 875         moldyn->k_sum+=moldyn->ekin;
 876         moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
 877         moldyn->k_avg=moldyn->k_sum/denom;
 878         moldyn->k2_avg=moldyn->k2_sum/denom;
 879         moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg);
 880
 881         /* potential energy */
 882         moldyn->v_sum+=moldyn->energy;
 883         moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
 884         moldyn->v_avg=moldyn->v_sum/denom;
 885         moldyn->v2_avg=moldyn->v2_sum/denom;
 886         moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg);
 887
 888         /* temperature */
 889         moldyn->t_sum+=moldyn->t;
 890         moldyn->t_avg=moldyn->t_sum/denom;
 891
 892         /* virial */
 893         moldyn->virial_sum+=moldyn->virial;
 894         moldyn->virial_avg=moldyn->virial_sum/denom;
 895         moldyn->gv_sum+=moldyn->gv;
 896         moldyn->gv_avg=moldyn->gv_sum/denom;
 897
 898         /* pressure */
 899         moldyn->p_sum+=moldyn->p;
 900         moldyn->p_avg=moldyn->p_sum/denom;
 901         moldyn->gp_sum+=moldyn->gp;
 902         moldyn->gp_avg=moldyn->gp_sum/denom;
 903
 904         return 0;
 905 }
 906
 907 int get_heat_capacity(t_moldyn *moldyn) {
 908
 909         double temp2,ighc;
 910
 911         /* averages needed for heat capacity calc */
 912         if(moldyn->total_steps<moldyn->avg_skip)
 913                 return 0;
 914
 915         /* (temperature average)^2 */
 916         temp2=moldyn->t_avg*moldyn->t_avg;
 917         printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
 918                moldyn->t_avg);
 919
 920         /* ideal gas contribution */
 921         ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
 922         printf("  ideal gas contribution: %f\n",
 923                ighc/moldyn->mass*KILOGRAM/JOULE);
 924
 925         /* specific heat for nvt ensemble */
 926         moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc;
 927         moldyn->c_v_nvt/=moldyn->mass;
 928
 929         /* specific heat for nve ensemble */
 930         moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2)));
 931         moldyn->c_v_nve/=moldyn->mass;
 932
 933         printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
 934         printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
 935 printf("  --> <dV2> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->count*K_B2*moldyn->t_avg*moldyn->t_avg*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM)));
 936
 937         return 0;
 938 }
 939
 940 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 941
 942         t_3dvec dim,*tp;
 943         double u_up,u_down,dv;
 944         double scale,p;
 945         t_atom *store;
 946
 947         /*
 948          * dU = - p dV
 949          *
 950          * => p = - dU/dV
 951          *
 952          */
 953
 954         scale=0.00001;
 955         dv=8*scale*scale*scale*moldyn->volume;
 956
 957         store=malloc(moldyn->count*sizeof(t_atom));
 958         if(store==NULL) {
 959                 printf("[moldyn] allocating store mem failed\n");
 960                 return -1;
 961         }
 962
 963         /* save unscaled potential energy + atom/dim configuration */
 964         memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
 965         dim=moldyn->dim;
 966
 967         /* scale up dimension and atom positions */
 968         scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
 969         scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
 970         link_cell_shutdown(moldyn);
 971         link_cell_init(moldyn,QUIET);
 972         potential_force_calc(moldyn);
 973         u_up=moldyn->energy;
 974
 975         /* restore atomic configuration + dim */
 976         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 977         moldyn->dim=dim;
 978
 979         /* scale down dimension and atom positions */
 980         scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
 981         scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
 982         link_cell_shutdown(moldyn);
 983         link_cell_init(moldyn,QUIET);
 984         potential_force_calc(moldyn);
 985         u_down=moldyn->energy;
 986
 987         /* calculate pressure */
 988         p=-(u_up-u_down)/dv;
 989 printf("-------> %.10f %.10f %f\n",u_up/EV/moldyn->count,u_down/EV/moldyn->count,p/BAR);
 990
 991         /* restore atomic configuration + dim */
 992         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 993         moldyn->dim=dim;
 994
 995         /* restore energy */
 996         potential_force_calc(moldyn);
 997
 998         link_cell_shutdown(moldyn);
 999         link_cell_init(moldyn,QUIET);
1000
1001         return p;
1002 }
1003
1004 double get_pressure(t_moldyn *moldyn) {
1005
1006         return moldyn->p;
1007
1008 }
1009
1010 int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1011
1012         t_3dvec *dim;
1013
1014         dim=&(moldyn->dim);
1015
1016         if(dir==SCALE_UP)
1017                 scale=1.0+scale;
1018
1019         if(dir==SCALE_DOWN)
1020                 scale=1.0-scale;
1021
1022         if(x) dim->x*=scale;
1023         if(y) dim->y*=scale;
1024         if(z) dim->z*=scale;
1025
1026         return 0;
1027 }
1028
1029 int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1030
1031         int i;
1032         t_3dvec *r;
1033
1034         if(dir==SCALE_UP)
1035                 scale=1.0+scale;
1036
1037         if(dir==SCALE_DOWN)
1038                 scale=1.0-scale;
1039
1040         for(i=0;i<moldyn->count;i++) {
1041                 r=&(moldyn->atom[i].r);
1042                 if(x) r->x*=scale;
1043                 if(y) r->y*=scale;
1044                 if(z) r->z*=scale;
1045         }
1046
1047         return 0;
1048 }
1049
1050 int scale_volume(t_moldyn *moldyn) {
1051
1052         t_3dvec *dim,*vdim;
1053         double scale;
1054         t_linkcell *lc;
1055
1056         vdim=&(moldyn->vis.dim);
1057         dim=&(moldyn->dim);
1058         lc=&(moldyn->lc);
1059
1060         /* scaling factor */
1061         if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
1062                 scale=1.0-(moldyn->p_ref-moldyn->p)/moldyn->p_tc;
1063                 scale=pow(scale,ONE_THIRD);
1064         }
1065         else {
1066                 scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
1067         }
1068 moldyn->debug=scale;
1069
1070         /* scale the atoms and dimensions */
1071         scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1072         scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1073
1074         /* visualize dimensions */
1075         if(vdim->x!=0) {
1076                 vdim->x=dim->x;
1077                 vdim->y=dim->y;
1078                 vdim->z=dim->z;
1079         }
1080
1081         /* recalculate scaled volume */
1082         moldyn->volume=dim->x*dim->y*dim->z;
1083
1084         /* adjust/reinit linkcell */
1085         if(((int)(dim->x/moldyn->cutoff)!=lc->nx)||
1086            ((int)(dim->y/moldyn->cutoff)!=lc->ny)||
1087            ((int)(dim->z/moldyn->cutoff)!=lc->nx)) {
1088                 link_cell_shutdown(moldyn);
1089                 link_cell_init(moldyn,QUIET);
1090         } else {
1091                 lc->x*=scale;
1092                 lc->y*=scale;
1093                 lc->z*=scale;
1094         }
1095
1096         return 0;
1097
1098 }
1099
1100 double e_kin_calc(t_moldyn *moldyn) {
1101
1102         int i;
1103         t_atom *atom;
1104
1105         atom=moldyn->atom;
1106         moldyn->ekin=0.0;
1107
1108         for(i=0;i<moldyn->count;i++) {
1109                 atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
1110                 moldyn->ekin+=atom[i].ekin;
1111         }
1112
1113         return moldyn->ekin;
1114 }
1115
1116 double get_total_energy(t_moldyn *moldyn) {
1117
1118         return(moldyn->ekin+moldyn->energy);
1119 }
1120
1121 t_3dvec get_total_p(t_moldyn *moldyn) {
1122
1123         t_3dvec p,p_total;
1124         int i;
1125         t_atom *atom;
1126
1127         atom=moldyn->atom;
1128
1129         v3_zero(&p_total);
1130         for(i=0;i<moldyn->count;i++) {
1131                 v3_scale(&p,&(atom[i].v),atom[i].mass);
1132                 v3_add(&p_total,&p_total,&p);
1133         }
1134
1135         return p_total;
1136 }
1137
1138 double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
1139
1140         double tau;
1141
1142         /* nn_dist is the nearest neighbour distance */
1143
1144         tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
1145
1146         return tau;
1147 }
1148
1149 /*
1150  * numerical tricks
1151  */
1152
1153 /* linked list / cell method */
1154
1155 int link_cell_init(t_moldyn *moldyn,u8 vol) {
1156
1157         t_linkcell *lc;
1158         int i;
1159
1160         lc=&(moldyn->lc);
1161
1162         /* partitioning the md cell */
1163         lc->nx=moldyn->dim.x/moldyn->cutoff;
1164         lc->x=moldyn->dim.x/lc->nx;
1165         lc->ny=moldyn->dim.y/moldyn->cutoff;
1166         lc->y=moldyn->dim.y/lc->ny;
1167         lc->nz=moldyn->dim.z/moldyn->cutoff;
1168         lc->z=moldyn->dim.z/lc->nz;
1169
1170         lc->cells=lc->nx*lc->ny*lc->nz;
1171         lc->subcell=malloc(lc->cells*sizeof(t_list));
1172
1173         if(lc->cells<27)
1174                 printf("[moldyn] FATAL: less then 27 subcells!\n");
1175
1176         if(vol) {
1177                 printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
1178                 printf("  x: %d x %f A\n",lc->nx,lc->x);
1179                 printf("  y: %d x %f A\n",lc->ny,lc->y);
1180                 printf("  z: %d x %f A\n",lc->nz,lc->z);
1181         }
1182
1183         for(i=0;i<lc->cells;i++)
1184                 list_init_f(&(lc->subcell[i]));
1185
1186         link_cell_update(moldyn);
1187
1188         return 0;
1189 }
1190
1191 int link_cell_update(t_moldyn *moldyn) {
1192
1193         int count,i,j,k;
1194         int nx,ny;
1195         t_atom *atom;
1196         t_linkcell *lc;
1197         double x,y,z;
1198
1199         atom=moldyn->atom;
1200         lc=&(moldyn->lc);
1201
1202         nx=lc->nx;
1203         ny=lc->ny;
1204
1205         x=moldyn->dim.x/2;
1206         y=moldyn->dim.y/2;
1207         z=moldyn->dim.z/2;
1208
1209         for(i=0;i<lc->cells;i++)
1210                 list_destroy_f(&(lc->subcell[i]));
1211
1212         for(count=0;count<moldyn->count;count++) {
1213                 i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
1214                 j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
1215                 k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
1216                 list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
1217                                      &(atom[count]));
1218         }
1219
1220         return 0;
1221 }
1222
1223 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
1224
1225         t_linkcell *lc;
1226         int a;
1227         int count1,count2;
1228         int ci,cj,ck;
1229         int nx,ny,nz;
1230         int x,y,z;
1231         u8 bx,by,bz;
1232
1233         lc=&(moldyn->lc);
1234         nx=lc->nx;
1235         ny=lc->ny;
1236         nz=lc->nz;
1237         count1=1;
1238         count2=27;
1239         a=nx*ny;
1240
1241         cell[0]=lc->subcell[i+j*nx+k*a];
1242         for(ci=-1;ci<=1;ci++) {
1243                 bx=0;
1244                 x=i+ci;
1245                 if((x<0)||(x>=nx)) {
1246                         x=(x+nx)%nx;
1247                         bx=1;
1248                 }
1249                 for(cj=-1;cj<=1;cj++) {
1250                         by=0;
1251                         y=j+cj;
1252                         if((y<0)||(y>=ny)) {
1253                                 y=(y+ny)%ny;
1254                                 by=1;
1255                         }
1256                         for(ck=-1;ck<=1;ck++) {
1257                                 bz=0;
1258                                 z=k+ck;
1259                                 if((z<0)||(z>=nz)) {
1260                                         z=(z+nz)%nz;
1261                                         bz=1;
1262                                 }
1263                                 if(!(ci|cj|ck)) continue;
1264                                 if(bx|by|bz) {
1265                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
1266                                 }
1267                                 else {
1268                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
1269                                 }
1270                         }
1271                 }
1272         }
1273
1274         lc->dnlc=count1;
1275
1276         return count1;
1277 }
1278
1279 int link_cell_shutdown(t_moldyn *moldyn) {
1280
1281         int i;
1282         t_linkcell *lc;
1283
1284         lc=&(moldyn->lc);
1285
1286         for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
1287                 list_destroy_f(&(moldyn->lc.subcell[i]));
1288
1289         free(lc->subcell);
1290
1291         return 0;
1292 }
1293
1294 int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
1295
1296         int count;
1297         void *ptr;
1298         t_moldyn_schedule *schedule;
1299
1300         schedule=&(moldyn->schedule);
1301         count=++(schedule->total_sched);
1302
1303         ptr=realloc(schedule->runs,count*sizeof(int));
1304         if(!ptr) {
1305                 perror("[moldyn] realloc (runs)");
1306                 return -1;
1307         }
1308         schedule->runs=ptr;
1309         schedule->runs[count-1]=runs;
1310
1311         ptr=realloc(schedule->tau,count*sizeof(double));
1312         if(!ptr) {
1313                 perror("[moldyn] realloc (tau)");
1314                 return -1;
1315         }
1316         schedule->tau=ptr;
1317         schedule->tau[count-1]=tau;
1318
1319         printf("[moldyn] schedule added:\n");
1320         printf("  number: %d | runs: %d | tau: %f\n",count-1,runs,tau);
1321
1322
1323         return 0;
1324 }
1325
1326 int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
1327
1328         moldyn->schedule.hook=hook;
1329         moldyn->schedule.hook_params=hook_params;
1330
1331         return 0;
1332 }
1333
1334 /*
1335  *
1336  * 'integration of newtons equation' - algorithms
1337  *
1338  */
1339
1340 /* start the integration */
1341
1342 int moldyn_integrate(t_moldyn *moldyn) {
1343
1344         int i;
1345         unsigned int e,m,s,v,p,t;
1346         t_3dvec momentum;
1347         t_moldyn_schedule *sched;
1348         t_atom *atom;
1349         int fd;
1350         char dir[128];
1351         double ds;
1352         double energy_scale;
1353         struct timeval t1,t2;
1354         //double tp;
1355
1356         sched=&(moldyn->schedule);
1357         atom=moldyn->atom;
1358
1359         /* initialize linked cell method */
1360         link_cell_init(moldyn,VERBOSE);
1361
1362         /* logging & visualization */
1363         e=moldyn->ewrite;
1364         m=moldyn->mwrite;
1365         s=moldyn->swrite;
1366         v=moldyn->vwrite;
1367         p=moldyn->pwrite;
1368         t=moldyn->twrite;
1369
1370         /* sqaure of some variables */
1371         moldyn->tau_square=moldyn->tau*moldyn->tau;
1372         moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
1373
1374         /* get current time */
1375         gettimeofday(&t1,NULL);
1376
1377         /* calculate initial forces */
1378         potential_force_calc(moldyn);
1379 #ifdef DEBUG
1380 return 0;
1381 #endif
1382
1383         /* some stupid checks before we actually start calculating bullshit */
1384         if(moldyn->cutoff>0.5*moldyn->dim.x)
1385                 printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
1386         if(moldyn->cutoff>0.5*moldyn->dim.y)
1387                 printf("[moldyn] warning: cutoff > 0.5 x dim.y\n");
1388         if(moldyn->cutoff>0.5*moldyn->dim.z)
1389                 printf("[moldyn] warning: cutoff > 0.5 x dim.z\n");
1390         ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
1391         if(ds>0.05*moldyn->nnd)
1392                 printf("[moldyn] warning: forces too high / tau too small!\n");
1393
1394         /* zero absolute time */
1395         moldyn->time=0.0;
1396         moldyn->total_steps=0;
1397
1398         /* debugging, ignore */
1399         moldyn->debug=0;
1400
1401         /* tell the world */
1402         printf("[moldyn] integration start, go get a coffee ...\n");
1403
1404         /* executing the schedule */
1405         sched->count=0;
1406         while(sched->count<sched->total_sched) {
1407
1408                 /* setting amount of runs and finite time step size */
1409                 moldyn->tau=sched->tau[sched->count];
1410                 moldyn->tau_square=moldyn->tau*moldyn->tau;
1411                 moldyn->time_steps=sched->runs[sched->count];
1412
1413                 /* energy scaling factor (might change!) */
1414                 energy_scale=moldyn->count*EV;
1415
1416         /* integration according to schedule */
1417
1418         for(i=0;i<moldyn->time_steps;i++) {
1419
1420                 /* integration step */
1421                 moldyn->integrate(moldyn);
1422
1423                 /* calculate kinetic energy, temperature and pressure */
1424                 e_kin_calc(moldyn);
1425                 temperature_calc(moldyn);
1426                 virial_sum(moldyn);
1427                 pressure_calc(moldyn);
1428                 average_and_fluctuation_calc(moldyn);
1429
1430                 /* p/t scaling */
1431                 if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
1432                         scale_velocity(moldyn,FALSE);
1433                 if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
1434                         scale_volume(moldyn);
1435
1436                 /* check for log & visualization */
1437                 if(e) {
1438                         if(!(moldyn->total_steps%e))
1439                                 dprintf(moldyn->efd,
1440                                         "%f %f %f %f\n",
1441                                         moldyn->time,moldyn->ekin/energy_scale,
1442                                         moldyn->energy/energy_scale,
1443                                         get_total_energy(moldyn)/energy_scale);
1444                 }
1445                 if(m) {
1446                         if(!(moldyn->total_steps%m)) {
1447                                 momentum=get_total_p(moldyn);
1448                                 dprintf(moldyn->mfd,
1449                                         "%f %f %f %f %f\n",moldyn->time,
1450                                         momentum.x,momentum.y,momentum.z,
1451                                         v3_norm(&momentum));
1452                         }
1453                 }
1454                 if(p) {
1455                         if(!(moldyn->total_steps%p)) {
1456                                 dprintf(moldyn->pfd,
1457                                         "%f %f %f %f %f\n",moldyn->time,
1458                                          moldyn->p/BAR,moldyn->p_avg/BAR,
1459                                          moldyn->gp/BAR,moldyn->gp_avg/BAR);
1460                         }
1461                 }
1462                 if(t) {
1463                         if(!(moldyn->total_steps%t)) {
1464                                 dprintf(moldyn->tfd,
1465                                         "%f %f %f\n",
1466                                         moldyn->time,moldyn->t,moldyn->t_avg);
1467                         }
1468                 }
1469                 if(s) {
1470                         if(!(moldyn->total_steps%s)) {
1471                                 snprintf(dir,128,"%s/s-%07.f.save",
1472                                          moldyn->vlsdir,moldyn->time);
1473                                 fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
1474                                         S_IRUSR|S_IWUSR);
1475                                 if(fd<0) perror("[moldyn] save fd open");
1476                                 else {
1477                                         write(fd,moldyn,sizeof(t_moldyn));
1478                                         write(fd,moldyn->atom,
1479                                               moldyn->count*sizeof(t_atom));
1480                                 }
1481                                 close(fd);
1482                         }
1483                 }
1484                 if(v) {
1485                         if(!(moldyn->total_steps%v)) {
1486                                 visual_atoms(&(moldyn->vis),moldyn->time,
1487                                              moldyn->atom,moldyn->count);
1488                         }
1489                 }
1490
1491                 /* display progress */
1492                 if(!(moldyn->total_steps%10)) {
1493                         /* get current time */
1494                         gettimeofday(&t2,NULL);
1495
1496         printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
1497                sched->count,i,
1498                moldyn->t,moldyn->t_avg,
1499                moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
1500                moldyn->volume,
1501                (int)(t2.tv_sec-t1.tv_sec));
1502         fflush(stdout);
1503
1504                         /* copy over time */
1505                         t1=t2;
1506                 }
1507
1508                 /* increase absolute time */
1509                 moldyn->time+=moldyn->tau;
1510                 moldyn->total_steps+=1;
1511
1512         }
1513
1514                 /* check for hooks */
1515                 if(sched->hook) {
1516                         printf("\n ## schedule hook %d/%d start ##\n",
1517                                sched->count+1,sched->total_sched-1);
1518                         sched->hook(moldyn,sched->hook_params);
1519                         printf(" ## schedule hook end ##\n");
1520                 }
1521
1522                 /* increase the schedule counter */
1523                 sched->count+=1;
1524
1525         }
1526
1527         return 0;
1528 }
1529
1530 /* velocity verlet */
1531
1532 int velocity_verlet(t_moldyn *moldyn) {
1533
1534         int i,count;
1535         double tau,tau_square,h;
1536         t_3dvec delta;
1537         t_atom *atom;
1538
1539         atom=moldyn->atom;
1540         count=moldyn->count;
1541         tau=moldyn->tau;
1542         tau_square=moldyn->tau_square;
1543
1544         for(i=0;i<count;i++) {
1545                 /* new positions */
1546                 h=0.5/atom[i].mass;
1547                 v3_scale(&delta,&(atom[i].v),tau);
1548                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1549                 v3_scale(&delta,&(atom[i].f),h*tau_square);
1550                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1551                 check_per_bound(moldyn,&(atom[i].r));
1552
1553                 /* velocities [actually v(t+tau/2)] */
1554                 v3_scale(&delta,&(atom[i].f),h*tau);
1555                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1556         }
1557
1558         /* neighbour list update */
1559         link_cell_update(moldyn);
1560
1561         /* forces depending on chosen potential */
1562         potential_force_calc(moldyn);
1563
1564         for(i=0;i<count;i++) {
1565                 /* again velocities [actually v(t+tau)] */
1566                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
1567                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1568         }
1569
1570         return 0;
1571 }
1572
1573
1574 /*
1575  *
1576  * potentials & corresponding forces & virial routine
1577  *
1578  */
1579
1580 /* generic potential and force calculation */
1581
1582 int potential_force_calc(t_moldyn *moldyn) {
1583
1584         int i,j,k,count;
1585         t_atom *itom,*jtom,*ktom;
1586         t_virial *virial;
1587         t_linkcell *lc;
1588         t_list neighbour_i[27];
1589         t_list neighbour_i2[27];
1590         t_list *this,*that;
1591         u8 bc_ij,bc_ik;
1592         int dnlc;
1593
1594         count=moldyn->count;
1595         itom=moldyn->atom;
1596         lc=&(moldyn->lc);
1597
1598         /* reset energy */
1599         moldyn->energy=0.0;
1600
1601         /* reset global virial */
1602         memset(&(moldyn->gvir),0,sizeof(t_virial));
1603
1604         /* reset force, site energy and virial of every atom */
1605         for(i=0;i<count;i++) {
1606
1607                 /* reset force */
1608                 v3_zero(&(itom[i].f));
1609
1610                 /* reset virial */
1611                 virial=(&(itom[i].virial));
1612                 virial->xx=0.0;
1613                 virial->yy=0.0;
1614                 virial->zz=0.0;
1615                 virial->xy=0.0;
1616                 virial->xz=0.0;
1617                 virial->yz=0.0;
1618
1619                 /* reset site energy */
1620                 itom[i].e=0.0;
1621
1622         }
1623
1624         /* get energy, force and virial of every atom */
1625
1626         /* first (and only) loop over atoms i */
1627         for(i=0;i<count;i++) {
1628
1629                 /* single particle potential/force */
1630                 if(itom[i].attr&ATOM_ATTR_1BP)
1631                         if(moldyn->func1b)
1632                                 moldyn->func1b(moldyn,&(itom[i]));
1633
1634                 if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
1635                         continue;
1636
1637                 /* 2 body pair potential/force */
1638
1639                 link_cell_neighbour_index(moldyn,
1640                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
1641                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
1642                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
1643                                           neighbour_i);
1644
1645                 dnlc=lc->dnlc;
1646
1647                 /* first loop over atoms j */
1648                 if(moldyn->func2b) {
1649                         for(j=0;j<27;j++) {
1650
1651                                 this=&(neighbour_i[j]);
1652                                 list_reset_f(this);
1653
1654                                 if(this->start==NULL)
1655                                         continue;
1656
1657                                 bc_ij=(j<dnlc)?0:1;
1658
1659                                 do {
1660                                         jtom=this->current->data;
1661
1662                                         if(jtom==&(itom[i]))
1663                                                 continue;
1664
1665                                         if((jtom->attr&ATOM_ATTR_2BP)&
1666                                            (itom[i].attr&ATOM_ATTR_2BP)) {
1667                                                 moldyn->func2b(moldyn,
1668                                                                &(itom[i]),
1669                                                                jtom,
1670                                                                bc_ij);
1671                                         }
1672                                 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1673
1674                         }
1675                 }
1676
1677                 /* 3 body potential/force */
1678
1679                 if(!(itom[i].attr&ATOM_ATTR_3BP))
1680                         continue;
1681
1682                 /* copy the neighbour lists */
1683                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
1684
1685                 /* second loop over atoms j */
1686                 for(j=0;j<27;j++) {
1687
1688                         this=&(neighbour_i[j]);
1689                         list_reset_f(this);
1690
1691                         if(this->start==NULL)
1692                                 continue;
1693
1694                         bc_ij=(j<dnlc)?0:1;
1695
1696                         do {
1697                                 jtom=this->current->data;
1698
1699                                 if(jtom==&(itom[i]))
1700                                         continue;
1701
1702                                 if(!(jtom->attr&ATOM_ATTR_3BP))
1703                                         continue;
1704
1705                                 /* reset 3bp run */
1706                                 moldyn->run3bp=1;
1707
1708                                 if(moldyn->func3b_j1)
1709                                         moldyn->func3b_j1(moldyn,
1710                                                           &(itom[i]),
1711                                                           jtom,
1712                                                           bc_ij);
1713
1714                                 /* in first j loop, 3bp run can be skipped */
1715                                 if(!(moldyn->run3bp))
1716                                         continue;
1717
1718                                 /* first loop over atoms k */
1719                                 if(moldyn->func3b_k1) {
1720
1721                                 for(k=0;k<27;k++) {
1722
1723                                         that=&(neighbour_i2[k]);
1724                                         list_reset_f(that);
1725
1726                                         if(that->start==NULL)
1727                                                 continue;
1728
1729                                         bc_ik=(k<dnlc)?0:1;
1730
1731                                         do {
1732
1733                                                 ktom=that->current->data;
1734
1735                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
1736                                                         continue;
1737
1738                                                 if(ktom==jtom)
1739                                                         continue;
1740
1741                                                 if(ktom==&(itom[i]))
1742                                                         continue;
1743
1744                                                 moldyn->func3b_k1(moldyn,
1745                                                                   &(itom[i]),
1746                                                                   jtom,
1747                                                                   ktom,
1748                                                                   bc_ik|bc_ij);
1749
1750                                         } while(list_next_f(that)!=\
1751                                                 L_NO_NEXT_ELEMENT);
1752
1753                                 }
1754
1755                                 }
1756
1757                                 if(moldyn->func3b_j2)
1758                                         moldyn->func3b_j2(moldyn,
1759                                                           &(itom[i]),
1760                                                           jtom,
1761                                                           bc_ij);
1762
1763                                 /* second loop over atoms k */
1764                                 if(moldyn->func3b_k2) {
1765
1766                                 for(k=0;k<27;k++) {
1767
1768                                         that=&(neighbour_i2[k]);
1769                                         list_reset_f(that);
1770
1771                                         if(that->start==NULL)
1772                                                 continue;
1773
1774                                         bc_ik=(k<dnlc)?0:1;
1775
1776                                         do {
1777
1778                                                 ktom=that->current->data;
1779
1780                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
1781                                                         continue;
1782
1783                                                 if(ktom==jtom)
1784                                                         continue;
1785
1786                                                 if(ktom==&(itom[i]))
1787                                                         continue;
1788
1789                                                 moldyn->func3b_k2(moldyn,
1790                                                                   &(itom[i]),
1791                                                                   jtom,
1792                                                                   ktom,
1793                                                                   bc_ik|bc_ij);
1794
1795                                         } while(list_next_f(that)!=\
1796                                                 L_NO_NEXT_ELEMENT);
1797
1798                                 }
1799
1800                                 }
1801
1802                                 /* 2bp post function */
1803                                 if(moldyn->func3b_j3) {
1804                                         moldyn->func3b_j3(moldyn,
1805                                                           &(itom[i]),
1806                                                           jtom,bc_ij);
1807                                 }
1808
1809                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1810
1811                 }
1812
1813 #ifdef DEBUG
1814         //printf("\n\n");
1815 #endif
1816 #ifdef VDEBUG
1817         printf("\n\n");
1818 #endif
1819
1820         }
1821
1822 #ifdef DEBUG
1823         printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
1824 #endif
1825
1826         /* calculate global virial */
1827         for(i=0;i<count;i++) {
1828                 moldyn->gvir.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x;
1829                 moldyn->gvir.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y;
1830                 moldyn->gvir.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z;
1831                 moldyn->gvir.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x;
1832                 moldyn->gvir.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x;
1833                 moldyn->gvir.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y;
1834         }
1835
1836         return 0;
1837 }
1838
1839 /*
1840  * virial calculation
1841  */
1842
1843 //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
1844 int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
1845
1846         a->virial.xx+=f->x*d->x;
1847         a->virial.yy+=f->y*d->y;
1848         a->virial.zz+=f->z*d->z;
1849         a->virial.xy+=f->x*d->y;
1850         a->virial.xz+=f->x*d->z;
1851         a->virial.yz+=f->y*d->z;
1852
1853         return 0;
1854 }
1855
1856 /*
1857  * periodic boundary checking
1858  */
1859
1860 //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
1861 int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
1862
1863         double x,y,z;
1864         t_3dvec *dim;
1865
1866         dim=&(moldyn->dim);
1867
1868         x=dim->x/2;
1869         y=dim->y/2;
1870         z=dim->z/2;
1871
1872         if(moldyn->status&MOLDYN_STAT_PBX) {
1873                 if(a->x>=x) a->x-=dim->x;
1874                 else if(-a->x>x) a->x+=dim->x;
1875         }
1876         if(moldyn->status&MOLDYN_STAT_PBY) {
1877                 if(a->y>=y) a->y-=dim->y;
1878                 else if(-a->y>y) a->y+=dim->y;
1879         }
1880         if(moldyn->status&MOLDYN_STAT_PBZ) {
1881                 if(a->z>=z) a->z-=dim->z;
1882                 else if(-a->z>z) a->z+=dim->z;
1883         }
1884
1885         return 0;
1886 }
1887
1888 /*
1889  * debugging / critical check functions
1890  */
1891
1892 int moldyn_bc_check(t_moldyn *moldyn) {
1893
1894         t_atom *atom;
1895         t_3dvec *dim;
1896         int i;
1897         double x;
1898         u8 byte;
1899         int j,k;
1900
1901         atom=moldyn->atom;
1902         dim=&(moldyn->dim);
1903         x=dim->x/2;
1904
1905         for(i=0;i<moldyn->count;i++) {
1906                 if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {
1907                         printf("FATAL: atom %d: x: %.20f (%.20f)\n",
1908                                i,atom[i].r.x,dim->x/2);
1909                         printf("diagnostic:\n");
1910                         printf("-----------\natom.r.x:\n");
1911                         for(j=0;j<8;j++) {
1912                                 memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
1913                                 for(k=0;k<8;k++)
1914                                         printf("%d%c",
1915                                         ((byte)&(1<<k))?1:0,
1916                                         (k==7)?'\n':'|');
1917                         }
1918                         printf("---------------\nx=dim.x/2:\n");
1919                         for(j=0;j<8;j++) {
1920                                 memcpy(&byte,(u8 *)(&x)+j,1);
1921                                 for(k=0;k<8;k++)
1922                                         printf("%d%c",
1923                                         ((byte)&(1<<k))?1:0,
1924                                         (k==7)?'\n':'|');
1925                         }
1926                         if(atom[i].r.x==x) printf("the same!\n");
1927                         else printf("different!\n");
1928                 }
1929                 if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
1930                         printf("FATAL: atom %d: y: %.20f (%.20f)\n",
1931                                i,atom[i].r.y,dim->y/2);
1932                 if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
1933                         printf("FATAL: atom %d: z: %.20f (%.20f)\n",
1934                                i,atom[i].r.z,dim->z/2);
1935         }
1936
1937         return 0;
1938 }
1939
1940 /*
1941  * restore function
1942  */
1943
1944 int moldyn_load(t_moldyn *moldyn) {
1945
1946         // later ...
1947
1948         return 0;
1949 }
1950
1951 /*
1952  * post processing functions
1953  */
1954
1955 int get_line(int fd,char *line,int max) {
1956
1957         int count,ret;
1958
1959         count=0;
1960
1961         while(1) {
1962                 if(count==max) return count;
1963                 ret=read(fd,line+count,1);
1964                 if(ret<=0) return ret;
1965                 if(line[count]=='\n') {
1966                         line[count]='\0';
1967                         return count+1;
1968                 }
1969                 count+=1;
1970         }
1971 }
1972
1973 int analyze_bonds(t_moldyn *moldyn) {
1974
1975
1976
1977
1978         return 0;
1979 }
1980