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