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 /* potential includes */
  24 #include "potentials/harmonic_oscillator.h"
  25 #include "potentials/lennard_jones.h"
  26 #include "potentials/albe.h"
  27 #ifdef TERSOFF_ORIG
  28 #include "potentials/tersoff_orig.h"
  29 #else
  30 #include "potentials/tersoff.h"
  31 #endif
  32
  33
  34 /*
  35  * global variables, pse and atom colors (only needed here)
  36  */
  37
  38 static char *pse_name[]={
  39         "*",
  40         "H",
  41         "He",
  42         "Li",
  43         "Be",
  44         "B",
  45         "C",
  46         "N",
  47         "O",
  48         "F",
  49         "Ne",
  50         "Na",
  51         "Mg",
  52         "Al",
  53         "Si",
  54         "P",
  55         "S",
  56         "Cl",
  57         "Ar",
  58 };
  59
  60 static char *pse_col[]={
  61         "*",
  62         "White",
  63         "He",
  64         "Li",
  65         "Be",
  66         "B",
  67         "Gray",
  68         "N",
  69         "Blue",
  70         "F",
  71         "Ne",
  72         "Na",
  73         "Mg",
  74         "Al",
  75         "Yellow",
  76         "P",
  77         "S",
  78         "Cl",
  79         "Ar",
  80 };
  81
  82 /*
  83 static double pse_mass[]={
  84         0,
  85         0,
  86         0,
  87         0,
  88         0,
  89         0,
  90         M_C,
  91         0,
  92         0,
  93         0,
  94         0,
  95         0,
  96         0,
  97         0,
  98         M_SI,
  99         0,
 100         0,
 101         0,
 102         0,
 103 };
 104
 105 static double pse_lc[]={
 106         0,
 107         0,
 108         0,
 109         0,
 110         0,
 111         0,
 112         LC_C,
 113         0,
 114         0,
 115         0,
 116         0,
 117         0,
 118         0,
 119         0,
 120         LC_SI,
 121         0,
 122         0,
 123         0,
 124         0,
 125 };
 126 */
 127
 128 /*
 129  * the moldyn functions
 130  */
 131
 132 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 133
 134         printf("[moldyn] init\n");
 135
 136         memset(moldyn,0,sizeof(t_moldyn));
 137
 138         moldyn->argc=argc;
 139         moldyn->args=argv;
 140
 141         rand_init(&(moldyn->random),NULL,1);
 142         moldyn->random.status|=RAND_STAT_VERBOSE;
 143
 144         return 0;
 145 }
 146
 147 int moldyn_shutdown(t_moldyn *moldyn) {
 148
 149         printf("[moldyn] shutdown\n");
 150
 151         moldyn_log_shutdown(moldyn);
 152         link_cell_shutdown(moldyn);
 153         rand_close(&(moldyn->random));
 154         free(moldyn->atom);
 155
 156         return 0;
 157 }
 158
 159 int set_int_alg(t_moldyn *moldyn,u8 algo) {
 160
 161         printf("[moldyn] integration algorithm: ");
 162
 163         switch(algo) {
 164                 case MOLDYN_INTEGRATE_VERLET:
 165                         moldyn->integrate=velocity_verlet;
 166                         printf("velocity verlet\n");
 167                         break;
 168                 default:
 169                         printf("unknown integration algorithm: %02x\n",algo);
 170                         printf("unknown\n");
 171                         return -1;
 172         }
 173
 174         return 0;
 175 }
 176
 177 int set_cutoff(t_moldyn *moldyn,double cutoff) {
 178
 179         moldyn->cutoff=cutoff;
 180         moldyn->cutoff_square=cutoff*cutoff;
 181
 182         printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff);
 183
 184         return 0;
 185 }
 186
 187 int set_temperature(t_moldyn *moldyn,double t_ref) {
 188
 189         moldyn->t_ref=t_ref;
 190
 191         printf("[moldyn] temperature [K]: %f\n",moldyn->t_ref);
 192
 193         return 0;
 194 }
 195
 196 int set_pressure(t_moldyn *moldyn,double p_ref) {
 197
 198         moldyn->p_ref=p_ref;
 199
 200         printf("[moldyn] pressure [bar]: %f\n",moldyn->p_ref/BAR);
 201
 202         return 0;
 203 }
 204
 205 int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc) {
 206
 207         moldyn->pt_scale&=(~(P_SCALE_MASK));
 208         moldyn->pt_scale|=ptype;
 209         moldyn->p_tc=ptc;
 210
 211         printf("[moldyn] p/t scaling:\n");
 212
 213         printf("  p: %s",ptype?"yes":"no ");
 214         if(ptype)
 215                 printf(" | type: %02x | factor: %f",ptype,ptc);
 216         printf("\n");
 217
 218         return 0;
 219 }
 220
 221 int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc) {
 222
 223         moldyn->pt_scale&=(~(T_SCALE_MASK));
 224         moldyn->pt_scale|=ttype;
 225         moldyn->t_tc=ttc;
 226
 227         printf("[moldyn] p/t scaling:\n");
 228
 229         printf("  t: %s",ttype?"yes":"no ");
 230         if(ttype)
 231                 printf(" | type: %02x | factor: %f",ttype,ttc);
 232         printf("\n");
 233
 234         return 0;
 235 }
 236
 237 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
 238
 239         moldyn->pt_scale=(ptype|ttype);
 240         moldyn->t_tc=ttc;
 241         moldyn->p_tc=ptc;
 242
 243         printf("[moldyn] p/t scaling:\n");
 244
 245         printf("  p: %s",ptype?"yes":"no ");
 246         if(ptype)
 247                 printf(" | type: %02x | factor: %f",ptype,ptc);
 248         printf("\n");
 249
 250         printf("  t: %s",ttype?"yes":"no ");
 251         if(ttype)
 252                 printf(" | type: %02x | factor: %f",ttype,ttc);
 253         printf("\n");
 254
 255         return 0;
 256 }
 257
 258 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
 259
 260         moldyn->dim.x=x;
 261         moldyn->dim.y=y;
 262         moldyn->dim.z=z;
 263
 264         moldyn->volume=x*y*z;
 265
 266         if(visualize) {
 267                 moldyn->vis.dim.x=x;
 268                 moldyn->vis.dim.y=y;
 269                 moldyn->vis.dim.z=z;
 270         }
 271
 272         printf("[moldyn] dimensions in A and A^3 respectively:\n");
 273         printf("  x: %f\n",moldyn->dim.x);
 274         printf("  y: %f\n",moldyn->dim.y);
 275         printf("  z: %f\n",moldyn->dim.z);
 276         printf("  volume: %f\n",moldyn->volume);
 277         printf("  visualize simulation box: %s\n",visualize?"yes":"no");
 278
 279         return 0;
 280 }
 281
 282 int set_nn_dist(t_moldyn *moldyn,double dist) {
 283
 284         moldyn->nnd=dist;
 285
 286         return 0;
 287 }
 288
 289 int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
 290
 291         printf("[moldyn] periodic boundary conditions:\n");
 292
 293         if(x)
 294                 moldyn->status|=MOLDYN_STAT_PBX;
 295
 296         if(y)
 297                 moldyn->status|=MOLDYN_STAT_PBY;
 298
 299         if(z)
 300                 moldyn->status|=MOLDYN_STAT_PBZ;
 301
 302         printf("  x: %s\n",x?"yes":"no");
 303         printf("  y: %s\n",y?"yes":"no");
 304         printf("  z: %s\n",z?"yes":"no");
 305
 306         return 0;
 307 }
 308
 309 int set_potential(t_moldyn *moldyn,u8 type) {
 310
 311         switch(type) {
 312                 case MOLDYN_POTENTIAL_TM:
 313                         moldyn->func1b=tersoff_mult_1bp;
 314                         moldyn->func3b_j1=tersoff_mult_3bp_j1;
 315                         moldyn->func3b_k1=tersoff_mult_3bp_k1;
 316                         moldyn->func3b_j2=tersoff_mult_3bp_j2;
 317                         moldyn->func3b_k2=tersoff_mult_3bp_k2;
 318                         // missing: check 2b bond func
 319                         break;
 320                 case MOLDYN_POTENTIAL_AM:
 321                         moldyn->func3b_j1=albe_mult_3bp_j1;
 322                         moldyn->func3b_k1=albe_mult_3bp_k1;
 323                         moldyn->func3b_j2=albe_mult_3bp_j2;
 324                         moldyn->func3b_k2=albe_mult_3bp_k2;
 325                         moldyn->check_2b_bond=albe_mult_check_2b_bond;
 326                         break;
 327                 case MOLDYN_POTENTIAL_HO:
 328                         moldyn->func2b=harmonic_oscillator;
 329                         moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
 330                         break;
 331                 case MOLDYN_POTENTIAL_LJ:
 332                         moldyn->func2b=lennard_jones;
 333                         moldyn->check_2b_bond=lennard_jones_check_2b_bond;
 334                         break;
 335                 default:
 336                         printf("[moldyn] set potential: unknown type %02x\n",
 337                                type);
 338                         return -1;
 339         }
 340
 341         return 0;
 342 }
 343
 344 int set_avg_skip(t_moldyn *moldyn,int skip) {
 345
 346         printf("[moldyn] skip %d steps before starting average calc\n",skip);
 347         moldyn->avg_skip=skip;
 348
 349         return 0;
 350 }
 351
 352 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
 353
 354         strncpy(moldyn->vlsdir,dir,127);
 355
 356         return 0;
 357 }
 358
 359 int moldyn_set_report(t_moldyn *moldyn,char *author,char *title) {
 360
 361         strncpy(moldyn->rauthor,author,63);
 362         strncpy(moldyn->rtitle,title,63);
 363
 364         return 0;
 365 }
 366
 367 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 368
 369         char filename[128];
 370         int ret;
 371
 372         printf("[moldyn] set log: ");
 373
 374         switch(type) {
 375                 case LOG_TOTAL_ENERGY:
 376                         moldyn->ewrite=timer;
 377                         snprintf(filename,127,"%s/energy",moldyn->vlsdir);
 378                         moldyn->efd=open(filename,
 379                                          O_WRONLY|O_CREAT|O_EXCL,
 380                                          S_IRUSR|S_IWUSR);
 381                         if(moldyn->efd<0) {
 382                                 perror("[moldyn] energy log fd open");
 383                                 return moldyn->efd;
 384                         }
 385                         dprintf(moldyn->efd,"# total energy log file\n");
 386                         printf("total energy (%d)\n",timer);
 387                         break;
 388                 case LOG_TOTAL_MOMENTUM:
 389                         moldyn->mwrite=timer;
 390                         snprintf(filename,127,"%s/momentum",moldyn->vlsdir);
 391                         moldyn->mfd=open(filename,
 392                                          O_WRONLY|O_CREAT|O_EXCL,
 393                                          S_IRUSR|S_IWUSR);
 394                         if(moldyn->mfd<0) {
 395                                 perror("[moldyn] momentum log fd open");
 396                                 return moldyn->mfd;
 397                         }
 398                         dprintf(moldyn->efd,"# total momentum log file\n");
 399                         printf("total momentum (%d)\n",timer);
 400                         break;
 401                 case LOG_PRESSURE:
 402                         moldyn->pwrite=timer;
 403                         snprintf(filename,127,"%s/pressure",moldyn->vlsdir);
 404                         moldyn->pfd=open(filename,
 405                                          O_WRONLY|O_CREAT|O_EXCL,
 406                                          S_IRUSR|S_IWUSR);
 407                         if(moldyn->pfd<0) {
 408                                 perror("[moldyn] pressure log file\n");
 409                                 return moldyn->pfd;
 410                         }
 411                         dprintf(moldyn->pfd,"# pressure log file\n");
 412                         printf("pressure (%d)\n",timer);
 413                         break;
 414                 case LOG_TEMPERATURE:
 415                         moldyn->twrite=timer;
 416                         snprintf(filename,127,"%s/temperature",moldyn->vlsdir);
 417                         moldyn->tfd=open(filename,
 418                                          O_WRONLY|O_CREAT|O_EXCL,
 419                                          S_IRUSR|S_IWUSR);
 420                         if(moldyn->tfd<0) {
 421                                 perror("[moldyn] temperature log file\n");
 422                                 return moldyn->tfd;
 423                         }
 424                         dprintf(moldyn->tfd,"# temperature log file\n");
 425                         printf("temperature (%d)\n",timer);
 426                         break;
 427                 case LOG_VOLUME:
 428                         moldyn->vwrite=timer;
 429                         snprintf(filename,127,"%s/volume",moldyn->vlsdir);
 430                         moldyn->vfd=open(filename,
 431                                          O_WRONLY|O_CREAT|O_EXCL,
 432                                          S_IRUSR|S_IWUSR);
 433                         if(moldyn->vfd<0) {
 434                                 perror("[moldyn] volume log file\n");
 435                                 return moldyn->vfd;
 436                         }
 437                         dprintf(moldyn->vfd,"# volume log file\n");
 438                         printf("volume (%d)\n",timer);
 439                         break;
 440                 case SAVE_STEP:
 441                         moldyn->swrite=timer;
 442                         printf("save file (%d)\n",timer);
 443                         break;
 444                 case VISUAL_STEP:
 445                         moldyn->awrite=timer;
 446                         ret=visual_init(moldyn,moldyn->vlsdir);
 447                         if(ret<0) {
 448                                 printf("[moldyn] visual init failure\n");
 449                                 return ret;
 450                         }
 451                         printf("visual file (%d)\n",timer);
 452                         break;
 453                 case CREATE_REPORT:
 454                         snprintf(filename,127,"%s/report.tex",moldyn->vlsdir);
 455                         moldyn->rfd=open(filename,
 456                                          O_WRONLY|O_CREAT|O_EXCL,
 457                                          S_IRUSR|S_IWUSR);
 458                         if(moldyn->rfd<0) {
 459                                 perror("[moldyn] report fd open");
 460                                 return moldyn->rfd;
 461                         }
 462                         printf("report -> ");
 463                         if(moldyn->efd) {
 464                                 snprintf(filename,127,"%s/e_plot.scr",
 465                                          moldyn->vlsdir);
 466                                 moldyn->epfd=open(filename,
 467                                                  O_WRONLY|O_CREAT|O_EXCL,
 468                                                  S_IRUSR|S_IWUSR);
 469                                 if(moldyn->epfd<0) {
 470                                         perror("[moldyn] energy plot fd open");
 471                                         return moldyn->epfd;
 472                                 }
 473                                 dprintf(moldyn->epfd,e_plot_script);
 474                                 close(moldyn->epfd);
 475                                 printf("energy ");
 476                         }
 477                         if(moldyn->pfd) {
 478                                 snprintf(filename,127,"%s/pressure_plot.scr",
 479                                          moldyn->vlsdir);
 480                                 moldyn->ppfd=open(filename,
 481                                                   O_WRONLY|O_CREAT|O_EXCL,
 482                                                   S_IRUSR|S_IWUSR);
 483                                 if(moldyn->ppfd<0) {
 484                                         perror("[moldyn] p plot fd open");
 485                                         return moldyn->ppfd;
 486                                 }
 487                                 dprintf(moldyn->ppfd,pressure_plot_script);
 488                                 close(moldyn->ppfd);
 489                                 printf("pressure ");
 490                         }
 491                         if(moldyn->tfd) {
 492                                 snprintf(filename,127,"%s/temperature_plot.scr",
 493                                          moldyn->vlsdir);
 494                                 moldyn->tpfd=open(filename,
 495                                                   O_WRONLY|O_CREAT|O_EXCL,
 496                                                   S_IRUSR|S_IWUSR);
 497                                 if(moldyn->tpfd<0) {
 498                                         perror("[moldyn] t plot fd open");
 499                                         return moldyn->tpfd;
 500                                 }
 501                                 dprintf(moldyn->tpfd,temperature_plot_script);
 502                                 close(moldyn->tpfd);
 503                                 printf("temperature ");
 504                         }
 505                         dprintf(moldyn->rfd,report_start,
 506                                 moldyn->rauthor,moldyn->rtitle);
 507                         printf("\n");
 508                         break;
 509                 default:
 510                         printf("unknown log type: %02x\n",type);
 511                         return -1;
 512         }
 513
 514         return 0;
 515 }
 516
 517 int moldyn_log_shutdown(t_moldyn *moldyn) {
 518
 519         char sc[256];
 520
 521         printf("[moldyn] log shutdown\n");
 522         if(moldyn->efd) {
 523                 close(moldyn->efd);
 524                 if(moldyn->rfd) {
 525                         dprintf(moldyn->rfd,report_energy);
 526                         snprintf(sc,255,"cd %s && gnuplot e_plot.scr",
 527                                  moldyn->vlsdir);
 528                         system(sc);
 529                 }
 530         }
 531         if(moldyn->mfd) close(moldyn->mfd);
 532         if(moldyn->pfd) {
 533                 close(moldyn->pfd);
 534                 if(moldyn->rfd)
 535                         dprintf(moldyn->rfd,report_pressure);
 536                         snprintf(sc,255,"cd %s && gnuplot pressure_plot.scr",
 537                                  moldyn->vlsdir);
 538                         system(sc);
 539         }
 540         if(moldyn->tfd) {
 541                 close(moldyn->tfd);
 542                 if(moldyn->rfd)
 543                         dprintf(moldyn->rfd,report_temperature);
 544                         snprintf(sc,255,"cd %s && gnuplot temperature_plot.scr",
 545                                  moldyn->vlsdir);
 546                         system(sc);
 547         }
 548         if(moldyn->rfd) {
 549                 dprintf(moldyn->rfd,report_end);
 550                 close(moldyn->rfd);
 551                 snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
 552                          moldyn->vlsdir);
 553                 system(sc);
 554                 snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
 555                          moldyn->vlsdir);
 556                 system(sc);
 557                 snprintf(sc,255,"cd %s && dvipdf report >/dev/null 2>&1",
 558                          moldyn->vlsdir);
 559                 system(sc);
 560         }
 561
 562         return 0;
 563 }
 564
 565 /*
 566  * creating lattice functions
 567  */
 568
 569 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 570                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) {
 571
 572         int new,count;
 573         int ret;
 574         t_3dvec orig;
 575         void *ptr;
 576         t_atom *atom;
 577         char name[16];
 578
 579         new=a*b*c;
 580         count=moldyn->count;
 581
 582         /* how many atoms do we expect */
 583         if(type==CUBIC) new*=1;
 584         if(type==FCC) new*=4;
 585         if(type==DIAMOND) new*=8;
 586
 587         /* allocate space for atoms */
 588         ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom));
 589         if(!ptr) {
 590                 perror("[moldyn] realloc (create lattice)");
 591                 return -1;
 592         }
 593         moldyn->atom=ptr;
 594         atom=&(moldyn->atom[count]);
 595
 596         /* no atoms on the boundaries (only reason: it looks better!) */
 597         if(!origin) {
 598                 orig.x=0.5*lc;
 599                 orig.y=0.5*lc;
 600                 orig.z=0.5*lc;
 601         }
 602         else {
 603                 orig.x=origin->x;
 604                 orig.y=origin->y;
 605                 orig.z=origin->z;
 606         }
 607
 608         switch(type) {
 609                 case CUBIC:
 610                         set_nn_dist(moldyn,lc);
 611                         ret=cubic_init(a,b,c,lc,atom,&orig);
 612                         strcpy(name,"cubic");
 613                         break;
 614                 case FCC:
 615                         if(!origin)
 616                                 v3_scale(&orig,&orig,0.5);
 617                         set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
 618                         ret=fcc_init(a,b,c,lc,atom,&orig);
 619                         strcpy(name,"fcc");
 620                         break;
 621                 case DIAMOND:
 622                         if(!origin)
 623                                 v3_scale(&orig,&orig,0.25);
 624                         set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
 625                         ret=diamond_init(a,b,c,lc,atom,&orig);
 626                         strcpy(name,"diamond");
 627                         break;
 628                 default:
 629                         printf("unknown lattice type (%02x)\n",type);
 630                         return -1;
 631         }
 632
 633         /* debug */
 634         if(ret!=new) {
 635                 printf("[moldyn] creating lattice failed\n");
 636                 printf("  amount of atoms\n");
 637                 printf("  - expected: %d\n",new);
 638                 printf("  - created: %d\n",ret);
 639                 return -1;
 640         }
 641
 642         moldyn->count+=new;
 643         printf("[moldyn] created %s lattice with %d atoms\n",name,new);
 644
 645         for(ret=0;ret<new;ret++) {
 646                 atom[ret].element=element;
 647                 atom[ret].mass=mass;
 648                 atom[ret].attr=attr;
 649                 atom[ret].brand=brand;
 650                 atom[ret].tag=count+ret;
 651                 check_per_bound(moldyn,&(atom[ret].r));
 652                 atom[ret].r_0=atom[ret].r;
 653         }
 654
 655         /* update total system mass */
 656         total_mass_calc(moldyn);
 657
 658         return ret;
 659 }
 660
 661 int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 662              t_3dvec *r,t_3dvec *v) {
 663
 664         t_atom *atom;
 665         void *ptr;
 666         int count;
 667
 668         atom=moldyn->atom;
 669         count=(moldyn->count)++;        // asshole style!
 670
 671         ptr=realloc(atom,(count+1)*sizeof(t_atom));
 672         if(!ptr) {
 673                 perror("[moldyn] realloc (add atom)");
 674                 return -1;
 675         }
 676         moldyn->atom=ptr;
 677
 678         atom=moldyn->atom;
 679
 680         /* initialize new atom */
 681         memset(&(atom[count]),0,sizeof(t_atom));
 682         atom[count].r=*r;
 683         atom[count].v=*v;
 684         atom[count].element=element;
 685         atom[count].mass=mass;
 686         atom[count].brand=brand;
 687         atom[count].tag=count;
 688         atom[count].attr=attr;
 689         check_per_bound(moldyn,&(atom[count].r));
 690         atom[count].r_0=atom[count].r;
 691
 692         /* update total system mass */
 693         total_mass_calc(moldyn);
 694
 695         return 0;
 696 }
 697
 698 int del_atom(t_moldyn *moldyn,int tag) {
 699
 700         t_atom *new,*old;
 701         int cnt;
 702
 703         old=moldyn->atom;
 704
 705         new=(t_atom *)malloc((moldyn->count-1)*sizeof(t_atom));
 706         if(!new) {
 707                 perror("[moldyn]malloc (del atom)");
 708                 return -1;
 709         }
 710
 711         for(cnt=0;cnt<tag;cnt++)
 712                 new[cnt]=old[cnt];
 713
 714         for(cnt=tag+1;cnt<moldyn->count;cnt++) {
 715                 new[cnt-1]=old[cnt];
 716                 new[cnt-1].tag=cnt-1;
 717         }
 718
 719         moldyn->count-=1;
 720         moldyn->atom=new;
 721
 722         free(old);
 723
 724         return 0;
 725 }
 726
 727 /* cubic init */
 728 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 729
 730         int count;
 731         t_3dvec r;
 732         int i,j,k;
 733         t_3dvec o;
 734
 735         count=0;
 736         if(origin)
 737                 v3_copy(&o,origin);
 738         else
 739                 v3_zero(&o);
 740
 741         r.x=o.x;
 742         for(i=0;i<a;i++) {
 743                 r.y=o.y;
 744                 for(j=0;j<b;j++) {
 745                         r.z=o.z;
 746                         for(k=0;k<c;k++) {
 747                                 v3_copy(&(atom[count].r),&r);
 748                                 count+=1;
 749                                 r.z+=lc;
 750                         }
 751                         r.y+=lc;
 752                 }
 753                 r.x+=lc;
 754         }
 755
 756         for(i=0;i<count;i++) {
 757                 atom[i].r.x-=(a*lc)/2.0;
 758                 atom[i].r.y-=(b*lc)/2.0;
 759                 atom[i].r.z-=(c*lc)/2.0;
 760         }
 761
 762         return count;
 763 }
 764
 765 /* fcc lattice init */
 766 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 767
 768         int count;
 769         int i,j,k,l;
 770         t_3dvec o,r,n;
 771         t_3dvec basis[3];
 772
 773         count=0;
 774         if(origin)
 775                 v3_copy(&o,origin);
 776         else
 777                 v3_zero(&o);
 778
 779         /* construct the basis */
 780         memset(basis,0,3*sizeof(t_3dvec));
 781         basis[0].x=0.5*lc;
 782         basis[0].y=0.5*lc;
 783         basis[1].x=0.5*lc;
 784         basis[1].z=0.5*lc;
 785         basis[2].y=0.5*lc;
 786         basis[2].z=0.5*lc;
 787
 788         /* fill up the room */
 789         r.x=o.x;
 790         for(i=0;i<a;i++) {
 791                 r.y=o.y;
 792                 for(j=0;j<b;j++) {
 793                         r.z=o.z;
 794                         for(k=0;k<c;k++) {
 795                                 /* first atom */
 796                                 v3_copy(&(atom[count].r),&r);
 797                                 count+=1;
 798                                 r.z+=lc;
 799                                 /* the three face centered atoms */
 800                                 for(l=0;l<3;l++) {
 801                                         v3_add(&n,&r,&basis[l]);
 802                                         v3_copy(&(atom[count].r),&n);
 803                                         count+=1;
 804                                 }
 805                         }
 806                         r.y+=lc;
 807                 }
 808                 r.x+=lc;
 809         }
 810
 811         /* coordinate transformation */
 812         for(i=0;i<count;i++) {
 813                 atom[i].r.x-=(a*lc)/2.0;
 814                 atom[i].r.y-=(b*lc)/2.0;
 815                 atom[i].r.z-=(c*lc)/2.0;
 816         }
 817
 818         return count;
 819 }
 820
 821 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 822
 823         int count;
 824         t_3dvec o;
 825
 826         count=fcc_init(a,b,c,lc,atom,origin);
 827
 828         o.x=0.25*lc;
 829         o.y=0.25*lc;
 830         o.z=0.25*lc;
 831
 832         if(origin) v3_add(&o,&o,origin);
 833
 834         count+=fcc_init(a,b,c,lc,&atom[count],&o);
 835
 836         return count;
 837 }
 838
 839 int destroy_atoms(t_moldyn *moldyn) {
 840
 841         if(moldyn->atom) free(moldyn->atom);
 842
 843         return 0;
 844 }
 845
 846 int thermal_init(t_moldyn *moldyn,u8 equi_init) {
 847
 848         /*
 849          * - gaussian distribution of velocities
 850          * - zero total momentum
 851          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 852          */
 853
 854         int i;
 855         double v,sigma;
 856         t_3dvec p_total,delta;
 857         t_atom *atom;
 858         t_random *random;
 859
 860         atom=moldyn->atom;
 861         random=&(moldyn->random);
 862
 863         printf("[moldyn] thermal init (equi init: %s)\n",equi_init?"yes":"no");
 864
 865         /* gaussian distribution of velocities */
 866         v3_zero(&p_total);
 867         for(i=0;i<moldyn->count;i++) {
 868                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
 869                 /* x direction */
 870                 v=sigma*rand_get_gauss(random);
 871                 atom[i].v.x=v;
 872                 p_total.x+=atom[i].mass*v;
 873                 /* y direction */
 874                 v=sigma*rand_get_gauss(random);
 875                 atom[i].v.y=v;
 876                 p_total.y+=atom[i].mass*v;
 877                 /* z direction */
 878                 v=sigma*rand_get_gauss(random);
 879                 atom[i].v.z=v;
 880                 p_total.z+=atom[i].mass*v;
 881         }
 882
 883         /* zero total momentum */
 884         v3_scale(&p_total,&p_total,1.0/moldyn->count);
 885         for(i=0;i<moldyn->count;i++) {
 886                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 887                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 888         }
 889
 890         /* velocity scaling */
 891         scale_velocity(moldyn,equi_init);
 892
 893         return 0;
 894 }
 895
 896 double total_mass_calc(t_moldyn *moldyn) {
 897
 898         int i;
 899
 900         moldyn->mass=0.0;
 901
 902         for(i=0;i<moldyn->count;i++)
 903                 moldyn->mass+=moldyn->atom[i].mass;
 904
 905         return moldyn->mass;
 906 }
 907
 908 double temperature_calc(t_moldyn *moldyn) {
 909
 910         /* assume up to date kinetic energy, which is 3/2 N k_B T */
 911
 912         moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
 913
 914         return moldyn->t;
 915 }
 916
 917 double get_temperature(t_moldyn *moldyn) {
 918
 919         return moldyn->t;
 920 }
 921
 922 int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 923
 924         int i;
 925         double e,scale;
 926         t_atom *atom;
 927         int count;
 928
 929         atom=moldyn->atom;
 930
 931         /*
 932          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 933          */
 934
 935         /* get kinetic energy / temperature & count involved atoms */
 936         e=0.0;
 937         count=0;
 938         for(i=0;i<moldyn->count;i++) {
 939                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) {
 940                         e+=atom[i].mass*v3_absolute_square(&(atom[i].v));
 941                         count+=1;
 942                 }
 943         }
 944         e*=0.5;
 945         if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN);
 946         else return 0;  /* no atoms involved in scaling! */
 947
 948         /* (temporary) hack for e,t = 0 */
 949         if(e==0.0) {
 950         moldyn->t=0.0;
 951                 if(moldyn->t_ref!=0.0) {
 952                         thermal_init(moldyn,equi_init);
 953                         return 0;
 954                 }
 955                 else
 956                         return 0; /* no scaling needed */
 957         }
 958
 959
 960         /* get scaling factor */
 961         scale=moldyn->t_ref/moldyn->t;
 962         if(equi_init&TRUE)
 963                 scale*=2.0;
 964         else
 965                 if(moldyn->pt_scale&T_SCALE_BERENDSEN)
 966                         scale=1.0+(scale-1.0)/moldyn->t_tc;
 967         scale=sqrt(scale);
 968
 969         /* velocity scaling */
 970         for(i=0;i<moldyn->count;i++) {
 971                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
 972                         v3_scale(&(atom[i].v),&(atom[i].v),scale);
 973         }
 974
 975         return 0;
 976 }
 977
 978 double ideal_gas_law_pressure(t_moldyn *moldyn) {
 979
 980         double p;
 981
 982         p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume;
 983
 984         return p;
 985 }
 986
 987 double virial_sum(t_moldyn *moldyn) {
 988
 989         int i;
 990         t_virial *virial;
 991
 992         /* virial (sum over atom virials) */
 993         moldyn->virial=0.0;
 994         moldyn->vir.xx=0.0;
 995         moldyn->vir.yy=0.0;
 996         moldyn->vir.zz=0.0;
 997         moldyn->vir.xy=0.0;
 998         moldyn->vir.xz=0.0;
 999         moldyn->vir.yz=0.0;
1000         for(i=0;i<moldyn->count;i++) {
1001                 virial=&(moldyn->atom[i].virial);
1002                 moldyn->virial+=(virial->xx+virial->yy+virial->zz);
1003                 moldyn->vir.xx+=virial->xx;
1004                 moldyn->vir.yy+=virial->yy;
1005                 moldyn->vir.zz+=virial->zz;
1006                 moldyn->vir.xy+=virial->xy;
1007                 moldyn->vir.xz+=virial->xz;
1008                 moldyn->vir.yz+=virial->yz;
1009         }
1010
1011         /* global virial (absolute coordinates) */
1012         virial=&(moldyn->gvir);
1013         moldyn->gv=virial->xx+virial->yy+virial->zz;
1014
1015         return moldyn->virial;
1016 }
1017
1018 double pressure_calc(t_moldyn *moldyn) {
1019
1020         /*
1021          * PV = NkT + <W>
1022          * with W = 1/3 sum_i f_i r_i (- skipped!)
1023          * virial = sum_i f_i r_i
1024          *
1025          * => P = (2 Ekin + virial) / (3V)
1026          */
1027
1028         /* assume up to date virial & up to date kinetic energy */
1029
1030         /* pressure (atom virials) */
1031         moldyn->p=2.0*moldyn->ekin+moldyn->virial;
1032         moldyn->p/=(3.0*moldyn->volume);
1033
1034         /* pressure (absolute coordinates) */
1035         moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
1036         moldyn->gp/=(3.0*moldyn->volume);
1037
1038         return moldyn->p;
1039 }
1040
1041 int average_reset(t_moldyn *moldyn) {
1042
1043         printf("[moldyn] average reset\n");
1044
1045         /* update skip value */
1046         moldyn->avg_skip=moldyn->total_steps;
1047
1048         /* kinetic energy */
1049         moldyn->k_sum=0.0;
1050         moldyn->k2_sum=0.0;
1051
1052         /* potential energy */
1053         moldyn->v_sum=0.0;
1054         moldyn->v2_sum=0.0;
1055
1056         /* temperature */
1057         moldyn->t_sum=0.0;
1058
1059         /* virial */
1060         moldyn->virial_sum=0.0;
1061         moldyn->gv_sum=0.0;
1062
1063         /* pressure */
1064         moldyn->p_sum=0.0;
1065         moldyn->gp_sum=0.0;
1066         moldyn->tp_sum=0.0;
1067
1068         return 0;
1069 }
1070
1071 int average_and_fluctuation_calc(t_moldyn *moldyn) {
1072
1073         int denom;
1074
1075         if(moldyn->total_steps<moldyn->avg_skip)
1076                 return 0;
1077
1078         denom=moldyn->total_steps+1-moldyn->avg_skip;
1079
1080         /* assume up to date energies, temperature, pressure etc */
1081
1082         /* kinetic energy */
1083         moldyn->k_sum+=moldyn->ekin;
1084         moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
1085         moldyn->k_avg=moldyn->k_sum/denom;
1086         moldyn->k2_avg=moldyn->k2_sum/denom;
1087         moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg);
1088
1089         /* potential energy */
1090         moldyn->v_sum+=moldyn->energy;
1091         moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
1092         moldyn->v_avg=moldyn->v_sum/denom;
1093         moldyn->v2_avg=moldyn->v2_sum/denom;
1094         moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg);
1095
1096         /* temperature */
1097         moldyn->t_sum+=moldyn->t;
1098         moldyn->t_avg=moldyn->t_sum/denom;
1099
1100         /* virial */
1101         moldyn->virial_sum+=moldyn->virial;
1102         moldyn->virial_avg=moldyn->virial_sum/denom;
1103         moldyn->gv_sum+=moldyn->gv;
1104         moldyn->gv_avg=moldyn->gv_sum/denom;
1105
1106         /* pressure */
1107         moldyn->p_sum+=moldyn->p;
1108         moldyn->p_avg=moldyn->p_sum/denom;
1109         moldyn->gp_sum+=moldyn->gp;
1110         moldyn->gp_avg=moldyn->gp_sum/denom;
1111         moldyn->tp_sum+=moldyn->tp;
1112         moldyn->tp_avg=moldyn->tp_sum/denom;
1113
1114         return 0;
1115 }
1116
1117 int get_heat_capacity(t_moldyn *moldyn) {
1118
1119         double temp2,ighc;
1120
1121         /* averages needed for heat capacity calc */
1122         if(moldyn->total_steps<moldyn->avg_skip)
1123                 return 0;
1124
1125         /* (temperature average)^2 */
1126         temp2=moldyn->t_avg*moldyn->t_avg;
1127         printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
1128                moldyn->t_avg);
1129
1130         /* ideal gas contribution */
1131         ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
1132         printf("  ideal gas contribution: %f\n",
1133                ighc/moldyn->mass*KILOGRAM/JOULE);
1134
1135         /* specific heat for nvt ensemble */
1136         moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc;
1137         moldyn->c_v_nvt/=moldyn->mass;
1138
1139         /* specific heat for nve ensemble */
1140         moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2)));
1141         moldyn->c_v_nve/=moldyn->mass;
1142
1143         printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
1144         printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
1145 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)));
1146
1147         return 0;
1148 }
1149
1150 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
1151
1152         t_3dvec dim;
1153         //t_3dvec *tp;
1154         double h,dv;
1155         double y0,y1;
1156         double su,sd;
1157         t_atom *store;
1158
1159         /*
1160          * dU = - p dV
1161          *
1162          * => p = - dU/dV
1163          *
1164          */
1165
1166         /* store atomic configuration + dimension */
1167         store=malloc(moldyn->count*sizeof(t_atom));
1168         if(store==NULL) {
1169                 printf("[moldyn] allocating store mem failed\n");
1170                 return -1;
1171         }
1172         memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
1173         dim=moldyn->dim;
1174
1175         /* x1, y1 */
1176         sd=0.00001;
1177         h=(1.0-sd)*(1.0-sd)*(1.0-sd);
1178         su=pow(2.0-h,ONE_THIRD)-1.0;
1179         dv=(1.0-h)*moldyn->volume;
1180
1181         /* scale up dimension and atom positions */
1182         scale_dim(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
1183         scale_atoms(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
1184         link_cell_shutdown(moldyn);
1185         link_cell_init(moldyn,QUIET);
1186         potential_force_calc(moldyn);
1187         y1=moldyn->energy;
1188
1189         /* restore atomic configuration + dim */
1190         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
1191         moldyn->dim=dim;
1192
1193         /* scale down dimension and atom positions */
1194         scale_dim(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
1195         scale_atoms(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
1196         link_cell_shutdown(moldyn);
1197         link_cell_init(moldyn,QUIET);
1198         potential_force_calc(moldyn);
1199         y0=moldyn->energy;
1200
1201         /* calculate pressure */
1202         moldyn->tp=-(y1-y0)/(2.0*dv);
1203
1204         /* restore atomic configuration */
1205         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
1206         moldyn->dim=dim;
1207         link_cell_shutdown(moldyn);
1208         link_cell_init(moldyn,QUIET);
1209         //potential_force_calc(moldyn);
1210
1211         /* free store buffer */
1212         if(store)
1213                 free(store);
1214
1215         return moldyn->tp;
1216 }
1217
1218 double get_pressure(t_moldyn *moldyn) {
1219
1220         return moldyn->p;
1221
1222 }
1223
1224 int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1225
1226         t_3dvec *dim;
1227
1228         dim=&(moldyn->dim);
1229
1230         if(dir==SCALE_UP)
1231                 scale=1.0+scale;
1232
1233         if(dir==SCALE_DOWN)
1234                 scale=1.0-scale;
1235
1236         if(x) dim->x*=scale;
1237         if(y) dim->y*=scale;
1238         if(z) dim->z*=scale;
1239
1240         return 0;
1241 }
1242
1243 int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1244
1245         int i;
1246         t_3dvec *r;
1247
1248         if(dir==SCALE_UP)
1249                 scale=1.0+scale;
1250
1251         if(dir==SCALE_DOWN)
1252                 scale=1.0-scale;
1253
1254         for(i=0;i<moldyn->count;i++) {
1255                 r=&(moldyn->atom[i].r);
1256                 if(x) r->x*=scale;
1257                 if(y) r->y*=scale;
1258                 if(z) r->z*=scale;
1259         }
1260
1261         return 0;
1262 }
1263
1264 int scale_volume(t_moldyn *moldyn) {
1265
1266         t_3dvec *dim,*vdim;
1267         double scale;
1268         t_linkcell *lc;
1269
1270         vdim=&(moldyn->vis.dim);
1271         dim=&(moldyn->dim);
1272         lc=&(moldyn->lc);
1273
1274         /* scaling factor */
1275         if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
1276                 scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc;
1277                 scale=pow(scale,ONE_THIRD);
1278         }
1279         else {
1280                 scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
1281         }
1282
1283         /* scale the atoms and dimensions */
1284         scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1285         scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1286
1287         /* visualize dimensions */
1288         if(vdim->x!=0) {
1289                 vdim->x=dim->x;
1290                 vdim->y=dim->y;
1291                 vdim->z=dim->z;
1292         }
1293
1294         /* recalculate scaled volume */
1295         moldyn->volume=dim->x*dim->y*dim->z;
1296
1297         /* adjust/reinit linkcell */
1298         if(((int)(dim->x/moldyn->cutoff)!=lc->nx)||
1299            ((int)(dim->y/moldyn->cutoff)!=lc->ny)||
1300            ((int)(dim->z/moldyn->cutoff)!=lc->nx)) {
1301                 link_cell_shutdown(moldyn);
1302                 link_cell_init(moldyn,QUIET);
1303         } else {
1304                 lc->x*=scale;
1305                 lc->y*=scale;
1306                 lc->z*=scale;
1307         }
1308
1309         return 0;
1310
1311 }
1312
1313 double e_kin_calc(t_moldyn *moldyn) {
1314
1315         int i;
1316         t_atom *atom;
1317
1318         atom=moldyn->atom;
1319         moldyn->ekin=0.0;
1320
1321         for(i=0;i<moldyn->count;i++) {
1322                 atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
1323                 moldyn->ekin+=atom[i].ekin;
1324         }
1325
1326         return moldyn->ekin;
1327 }
1328
1329 double get_total_energy(t_moldyn *moldyn) {
1330
1331         return(moldyn->ekin+moldyn->energy);
1332 }
1333
1334 t_3dvec get_total_p(t_moldyn *moldyn) {
1335
1336         t_3dvec p,p_total;
1337         int i;
1338         t_atom *atom;
1339
1340         atom=moldyn->atom;
1341
1342         v3_zero(&p_total);
1343         for(i=0;i<moldyn->count;i++) {
1344                 v3_scale(&p,&(atom[i].v),atom[i].mass);
1345                 v3_add(&p_total,&p_total,&p);
1346         }
1347
1348         return p_total;
1349 }
1350
1351 double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
1352
1353         double tau;
1354
1355         /* nn_dist is the nearest neighbour distance */
1356
1357         tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
1358
1359         return tau;
1360 }
1361
1362 /*
1363  * numerical tricks
1364  */
1365
1366 /* linked list / cell method */
1367
1368 int link_cell_init(t_moldyn *moldyn,u8 vol) {
1369
1370         t_linkcell *lc;
1371         int i;
1372
1373         lc=&(moldyn->lc);
1374
1375         /* partitioning the md cell */
1376         lc->nx=moldyn->dim.x/moldyn->cutoff;
1377         lc->x=moldyn->dim.x/lc->nx;
1378         lc->ny=moldyn->dim.y/moldyn->cutoff;
1379         lc->y=moldyn->dim.y/lc->ny;
1380         lc->nz=moldyn->dim.z/moldyn->cutoff;
1381         lc->z=moldyn->dim.z/lc->nz;
1382         lc->cells=lc->nx*lc->ny*lc->nz;
1383
1384 #ifdef STATIC_LISTS
1385         lc->subcell=malloc(lc->cells*sizeof(int*));
1386 #else
1387         lc->subcell=malloc(lc->cells*sizeof(t_list));
1388 #endif
1389
1390         if(lc->subcell==NULL) {
1391                 perror("[moldyn] cell init (malloc)");
1392                 return -1;
1393         }
1394
1395         if(lc->cells<27)
1396                 printf("[moldyn] FATAL: less then 27 subcells!\n");
1397
1398         if(vol) {
1399 #ifdef STATIC_LISTS
1400                 printf("[moldyn] initializing 'static' linked cells (%d)\n",
1401                        lc->cells);
1402 #else
1403                 printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
1404                        lc->cells);
1405 #endif
1406                 printf("  x: %d x %f A\n",lc->nx,lc->x);
1407                 printf("  y: %d x %f A\n",lc->ny,lc->y);
1408                 printf("  z: %d x %f A\n",lc->nz,lc->z);
1409         }
1410
1411 #ifdef STATIC_LISTS
1412         /* list init */
1413         for(i=0;i<lc->cells;i++) {
1414                 lc->subcell[i]=malloc((MAX_ATOMS_PER_LIST+1)*sizeof(int));
1415                 if(lc->subcell[i]==NULL) {
1416                         perror("[moldyn] list init (malloc)");
1417                         return -1;
1418                 }
1419                 /*
1420                 if(i==0)
1421                         printf(" ---> %d malloc %p (%p)\n",
1422                                i,lc->subcell[0],lc->subcell);
1423                 */
1424         }
1425 #else
1426         for(i=0;i<lc->cells;i++)
1427                 list_init_f(&(lc->subcell[i]));
1428 #endif
1429
1430         /* update the list */
1431         link_cell_update(moldyn);
1432
1433         return 0;
1434 }
1435
1436 int link_cell_update(t_moldyn *moldyn) {
1437
1438         int count,i,j,k;
1439         int nx,ny;
1440         t_atom *atom;
1441         t_linkcell *lc;
1442 #ifdef STATIC_LISTS
1443         int p;
1444 #endif
1445
1446         atom=moldyn->atom;
1447         lc=&(moldyn->lc);
1448
1449         nx=lc->nx;
1450         ny=lc->ny;
1451
1452         for(i=0;i<lc->cells;i++)
1453 #ifdef STATIC_LISTS
1454                 memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
1455 #else
1456                 list_destroy_f(&(lc->subcell[i]));
1457 #endif
1458
1459         for(count=0;count<moldyn->count;count++) {
1460                 i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
1461                 j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
1462                 k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
1463
1464 #ifdef STATIC_LISTS
1465                 p=0;
1466                 while(lc->subcell[i+j*nx+k*nx*ny][p]!=0)
1467                         p++;
1468
1469                 if(p>=MAX_ATOMS_PER_LIST) {
1470                         printf("[moldyn] FATAL: amount of atoms too high!\n");
1471                         return -1;
1472                 }
1473
1474                 lc->subcell[i+j*nx+k*nx*ny][p]=count;
1475 #else
1476                 list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
1477                                      &(atom[count]));
1478                 /*
1479                 if(j==0&&k==0)
1480                         printf(" ---> %d %d malloc %p (%p)\n",
1481                                i,count,lc->subcell[i].current,lc->subcell);
1482                 */
1483 #endif
1484         }
1485
1486         return 0;
1487 }
1488
1489 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
1490 #ifdef STATIC_LISTS
1491                               int **cell
1492 #else
1493                               t_list *cell
1494 #endif
1495                              ) {
1496
1497         t_linkcell *lc;
1498         int a;
1499         int count1,count2;
1500         int ci,cj,ck;
1501         int nx,ny,nz;
1502         int x,y,z;
1503         u8 bx,by,bz;
1504
1505         lc=&(moldyn->lc);
1506         nx=lc->nx;
1507         ny=lc->ny;
1508         nz=lc->nz;
1509         count1=1;
1510         count2=27;
1511         a=nx*ny;
1512
1513         if(i>=nx||j>=ny||k>=nz)
1514                 printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n",
1515                        i,nx,j,ny,k,nz);
1516
1517         cell[0]=lc->subcell[i+j*nx+k*a];
1518         for(ci=-1;ci<=1;ci++) {
1519                 bx=0;
1520                 x=i+ci;
1521                 if((x<0)||(x>=nx)) {
1522                         x=(x+nx)%nx;
1523                         bx=1;
1524                 }
1525                 for(cj=-1;cj<=1;cj++) {
1526                         by=0;
1527                         y=j+cj;
1528                         if((y<0)||(y>=ny)) {
1529                                 y=(y+ny)%ny;
1530                                 by=1;
1531                         }
1532                         for(ck=-1;ck<=1;ck++) {
1533                                 bz=0;
1534                                 z=k+ck;
1535                                 if((z<0)||(z>=nz)) {
1536                                         z=(z+nz)%nz;
1537                                         bz=1;
1538                                 }
1539                                 if(!(ci|cj|ck)) continue;
1540                                 if(bx|by|bz) {
1541                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
1542                                 }
1543                                 else {
1544                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
1545                                 }
1546                         }
1547                 }
1548         }
1549
1550         lc->dnlc=count1;
1551
1552         return count1;
1553 }
1554
1555 int link_cell_shutdown(t_moldyn *moldyn) {
1556
1557         int i;
1558         t_linkcell *lc;
1559
1560         lc=&(moldyn->lc);
1561
1562         for(i=0;i<lc->cells;i++) {
1563 #ifdef STATIC_LISTS
1564                 free(lc->subcell[i]);
1565 #else
1566                 //printf(" ---> %d free %p\n",i,lc->subcell[i].start);
1567                 list_destroy_f(&(lc->subcell[i]));
1568 #endif
1569         }
1570
1571         free(lc->subcell);
1572
1573         return 0;
1574 }
1575
1576 int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
1577
1578         int count;
1579         void *ptr;
1580         t_moldyn_schedule *schedule;
1581
1582         schedule=&(moldyn->schedule);
1583         count=++(schedule->total_sched);
1584
1585         ptr=realloc(schedule->runs,count*sizeof(int));
1586         if(!ptr) {
1587                 perror("[moldyn] realloc (runs)");
1588                 return -1;
1589         }
1590         schedule->runs=ptr;
1591         schedule->runs[count-1]=runs;
1592
1593         ptr=realloc(schedule->tau,count*sizeof(double));
1594         if(!ptr) {
1595                 perror("[moldyn] realloc (tau)");
1596                 return -1;
1597         }
1598         schedule->tau=ptr;
1599         schedule->tau[count-1]=tau;
1600
1601         printf("[moldyn] schedule added:\n");
1602         printf("  number: %d | runs: %d | tau: %f\n",count-1,runs,tau);
1603
1604
1605         return 0;
1606 }
1607
1608 int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
1609
1610         moldyn->schedule.hook=hook;
1611         moldyn->schedule.hook_params=hook_params;
1612
1613         return 0;
1614 }
1615
1616 /*
1617  *
1618  * 'integration of newtons equation' - algorithms
1619  *
1620  */
1621
1622 /* start the integration */
1623
1624 int moldyn_integrate(t_moldyn *moldyn) {
1625
1626         int i;
1627         unsigned int e,m,s,v,p,t,a;
1628         t_3dvec momentum;
1629         t_moldyn_schedule *sched;
1630         t_atom *atom;
1631         int fd;
1632         char dir[128];
1633         double ds;
1634         double energy_scale;
1635         struct timeval t1,t2;
1636         //double tp;
1637
1638         sched=&(moldyn->schedule);
1639         atom=moldyn->atom;
1640
1641         /* initialize linked cell method */
1642         link_cell_init(moldyn,VERBOSE);
1643
1644         /* logging & visualization */
1645         e=moldyn->ewrite;
1646         m=moldyn->mwrite;
1647         s=moldyn->swrite;
1648         v=moldyn->vwrite;
1649         a=moldyn->awrite;
1650         p=moldyn->pwrite;
1651         t=moldyn->twrite;
1652
1653         /* sqaure of some variables */
1654         moldyn->tau_square=moldyn->tau*moldyn->tau;
1655
1656         /* get current time */
1657         gettimeofday(&t1,NULL);
1658
1659         /* calculate initial forces */
1660         potential_force_calc(moldyn);
1661 #ifdef DEBUG
1662 //return 0;
1663 #endif
1664
1665         /* some stupid checks before we actually start calculating bullshit */
1666         if(moldyn->cutoff>0.5*moldyn->dim.x)
1667                 printf("[moldyn] WARNING: cutoff > 0.5 x dim.x\n");
1668         if(moldyn->cutoff>0.5*moldyn->dim.y)
1669                 printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n");
1670         if(moldyn->cutoff>0.5*moldyn->dim.z)
1671                 printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n");
1672         ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
1673         if(ds>0.05*moldyn->nnd)
1674                 printf("[moldyn] WARNING: forces too high / tau too small!\n");
1675
1676         /* zero absolute time */
1677         moldyn->time=0.0;
1678         moldyn->total_steps=0;
1679
1680         /* debugging, ignore */
1681         moldyn->debug=0;
1682
1683         /* tell the world */
1684         printf("[moldyn] integration start, go get a coffee ...\n");
1685
1686         /* executing the schedule */
1687         sched->count=0;
1688         while(sched->count<sched->total_sched) {
1689
1690                 /* setting amount of runs and finite time step size */
1691                 moldyn->tau=sched->tau[sched->count];
1692                 moldyn->tau_square=moldyn->tau*moldyn->tau;
1693                 moldyn->time_steps=sched->runs[sched->count];
1694
1695                 /* energy scaling factor (might change!) */
1696                 energy_scale=moldyn->count*EV;
1697
1698         /* integration according to schedule */
1699
1700         for(i=0;i<moldyn->time_steps;i++) {
1701
1702                 /* integration step */
1703                 moldyn->integrate(moldyn);
1704
1705                 /* calculate kinetic energy, temperature and pressure */
1706                 e_kin_calc(moldyn);
1707                 temperature_calc(moldyn);
1708                 virial_sum(moldyn);
1709                 pressure_calc(moldyn);
1710                 /*
1711                 thermodynamic_pressure_calc(moldyn);
1712                 printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
1713                        moldyn->tp/BAR);
1714                 */
1715
1716                 /* calculate fluctuations + averages */
1717                 average_and_fluctuation_calc(moldyn);
1718
1719                 /* p/t scaling */
1720                 if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
1721                         scale_velocity(moldyn,FALSE);
1722                 if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
1723                         scale_volume(moldyn);
1724
1725                 /* check for log & visualization */
1726                 if(e) {
1727                         if(!(moldyn->total_steps%e))
1728                                 dprintf(moldyn->efd,
1729                                         "%f %f %f %f\n",
1730                                         moldyn->time,moldyn->ekin/energy_scale,
1731                                         moldyn->energy/energy_scale,
1732                                         get_total_energy(moldyn)/energy_scale);
1733                 }
1734                 if(m) {
1735                         if(!(moldyn->total_steps%m)) {
1736                                 momentum=get_total_p(moldyn);
1737                                 dprintf(moldyn->mfd,
1738                                         "%f %f %f %f %f\n",moldyn->time,
1739                                         momentum.x,momentum.y,momentum.z,
1740                                         v3_norm(&momentum));
1741                         }
1742                 }
1743                 if(p) {
1744                         if(!(moldyn->total_steps%p)) {
1745                                 dprintf(moldyn->pfd,
1746                                         "%f %f %f %f %f %f %f\n",moldyn->time,
1747                                          moldyn->p/BAR,moldyn->p_avg/BAR,
1748                                          moldyn->gp/BAR,moldyn->gp_avg/BAR,
1749                                          moldyn->tp/BAR,moldyn->tp_avg/BAR);
1750                         }
1751                 }
1752                 if(t) {
1753                         if(!(moldyn->total_steps%t)) {
1754                                 dprintf(moldyn->tfd,
1755                                         "%f %f %f\n",
1756                                         moldyn->time,moldyn->t,moldyn->t_avg);
1757                         }
1758                 }
1759                 if(v) {
1760                         if(!(moldyn->total_steps%v)) {
1761                                 dprintf(moldyn->vfd,
1762                                         "%f %f\n",moldyn->time,moldyn->volume);
1763                         }
1764                 }
1765                 if(s) {
1766                         if(!(moldyn->total_steps%s)) {
1767                                 snprintf(dir,128,"%s/s-%07.f.save",
1768                                          moldyn->vlsdir,moldyn->time);
1769                                 fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
1770                                         S_IRUSR|S_IWUSR);
1771                                 if(fd<0) perror("[moldyn] save fd open");
1772                                 else {
1773                                         write(fd,moldyn,sizeof(t_moldyn));
1774                                         write(fd,moldyn->atom,
1775                                               moldyn->count*sizeof(t_atom));
1776                                 }
1777                                 close(fd);
1778                         }
1779                 }
1780                 if(a) {
1781                         if(!(moldyn->total_steps%a)) {
1782                                 visual_atoms(moldyn);
1783                         }
1784                 }
1785
1786                 /* display progress */
1787                 //if(!(moldyn->total_steps%10)) {
1788                         /* get current time */
1789                         gettimeofday(&t2,NULL);
1790
1791 printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
1792        sched->count,i,moldyn->total_steps,
1793        moldyn->t,moldyn->t_avg,
1794        moldyn->p/BAR,moldyn->p_avg/BAR,
1795        moldyn->volume,
1796        (int)(t2.tv_sec-t1.tv_sec));
1797
1798                         fflush(stdout);
1799
1800                         /* copy over time */
1801                         t1=t2;
1802                 //}
1803
1804                 /* increase absolute time */
1805                 moldyn->time+=moldyn->tau;
1806                 moldyn->total_steps+=1;
1807
1808         }
1809
1810                 /* check for hooks */
1811                 if(sched->hook) {
1812                         printf("\n ## schedule hook %d start ##\n",
1813                                sched->count);
1814                         sched->hook(moldyn,sched->hook_params);
1815                         printf(" ## schedule hook end ##\n");
1816                 }
1817
1818                 /* increase the schedule counter */
1819                 sched->count+=1;
1820
1821         }
1822
1823         return 0;
1824 }
1825
1826 /* velocity verlet */
1827
1828 int velocity_verlet(t_moldyn *moldyn) {
1829
1830         int i,count;
1831         double tau,tau_square,h;
1832         t_3dvec delta;
1833         t_atom *atom;
1834
1835         atom=moldyn->atom;
1836         count=moldyn->count;
1837         tau=moldyn->tau;
1838         tau_square=moldyn->tau_square;
1839
1840         for(i=0;i<count;i++) {
1841                 /* check whether fixed atom */
1842                 if(atom[i].attr&ATOM_ATTR_FP)
1843                         continue;
1844                 /* new positions */
1845                 h=0.5/atom[i].mass;
1846                 v3_scale(&delta,&(atom[i].v),tau);
1847                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1848                 v3_scale(&delta,&(atom[i].f),h*tau_square);
1849                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1850                 check_per_bound(moldyn,&(atom[i].r));
1851
1852                 /* velocities [actually v(t+tau/2)] */
1853                 v3_scale(&delta,&(atom[i].f),h*tau);
1854                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1855         }
1856
1857         /* criticial check */
1858         moldyn_bc_check(moldyn);
1859
1860         /* neighbour list update */
1861         link_cell_update(moldyn);
1862
1863         /* forces depending on chosen potential */
1864         potential_force_calc(moldyn);
1865
1866         for(i=0;i<count;i++) {
1867                 /* check whether fixed atom */
1868                 if(atom[i].attr&ATOM_ATTR_FP)
1869                         continue;
1870                 /* again velocities [actually v(t+tau)] */
1871                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
1872                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1873         }
1874
1875         return 0;
1876 }
1877
1878
1879 /*
1880  *
1881  * potentials & corresponding forces & virial routine
1882  *
1883  */
1884
1885 /* generic potential and force calculation */
1886
1887 int potential_force_calc(t_moldyn *moldyn) {
1888
1889         int i,j,k,count;
1890         t_atom *itom,*jtom,*ktom;
1891         t_virial *virial;
1892         t_linkcell *lc;
1893 #ifdef STATIC_LISTS
1894         int *neighbour_i[27];
1895         int p,q;
1896         t_atom *atom;
1897 #else
1898         t_list neighbour_i[27];
1899         t_list neighbour_i2[27];
1900         t_list *this,*that;
1901 #endif
1902         u8 bc_ij,bc_ik;
1903         int dnlc;
1904
1905         count=moldyn->count;
1906         itom=moldyn->atom;
1907         lc=&(moldyn->lc);
1908 #ifdef STATIC_LISTS
1909         atom=moldyn->atom;
1910 #endif
1911
1912         /* reset energy */
1913         moldyn->energy=0.0;
1914
1915         /* reset global virial */
1916         memset(&(moldyn->gvir),0,sizeof(t_virial));
1917
1918         /* reset force, site energy and virial of every atom */
1919         for(i=0;i<count;i++) {
1920
1921                 /* reset force */
1922                 v3_zero(&(itom[i].f));
1923
1924                 /* reset virial */
1925                 virial=(&(itom[i].virial));
1926                 virial->xx=0.0;
1927                 virial->yy=0.0;
1928                 virial->zz=0.0;
1929                 virial->xy=0.0;
1930                 virial->xz=0.0;
1931                 virial->yz=0.0;
1932
1933                 /* reset site energy */
1934                 itom[i].e=0.0;
1935
1936         }
1937
1938         /* get energy, force and virial of every atom */
1939
1940         /* first (and only) loop over atoms i */
1941         for(i=0;i<count;i++) {
1942
1943                 /* single particle potential/force */
1944                 if(itom[i].attr&ATOM_ATTR_1BP)
1945                         if(moldyn->func1b)
1946                                 moldyn->func1b(moldyn,&(itom[i]));
1947
1948                 if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
1949                         continue;
1950
1951                 /* 2 body pair potential/force */
1952
1953                 link_cell_neighbour_index(moldyn,
1954                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
1955                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
1956                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
1957                                           neighbour_i);
1958
1959                 dnlc=lc->dnlc;
1960
1961                 /* first loop over atoms j */
1962                 if(moldyn->func2b) {
1963                         for(j=0;j<27;j++) {
1964
1965                                 bc_ij=(j<dnlc)?0:1;
1966 #ifdef STATIC_LISTS
1967                                 p=0;
1968
1969                                 while(neighbour_i[j][p]!=0) {
1970
1971                                         jtom=&(atom[neighbour_i[j][p]]);
1972                                         p++;
1973
1974                                         if(jtom==&(itom[i]))
1975                                                 continue;
1976
1977                                         if((jtom->attr&ATOM_ATTR_2BP)&
1978                                            (itom[i].attr&ATOM_ATTR_2BP)) {
1979                                                 moldyn->func2b(moldyn,
1980                                                                &(itom[i]),
1981                                                                jtom,
1982                                                                bc_ij);
1983                                         }
1984                                 }
1985 #else
1986                                 this=&(neighbour_i[j]);
1987                                 list_reset_f(this);
1988
1989                                 if(this->start==NULL)
1990                                         continue;
1991
1992                                 do {
1993                                         jtom=this->current->data;
1994
1995                                         if(jtom==&(itom[i]))
1996                                                 continue;
1997
1998                                         if((jtom->attr&ATOM_ATTR_2BP)&
1999                                            (itom[i].attr&ATOM_ATTR_2BP)) {
2000                                                 moldyn->func2b(moldyn,
2001                                                                &(itom[i]),
2002                                                                jtom,
2003                                                                bc_ij);
2004                                         }
2005                                 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2006 #endif
2007
2008                         }
2009                 }
2010
2011                 /* 3 body potential/force */
2012
2013                 if(!(itom[i].attr&ATOM_ATTR_3BP))
2014                         continue;
2015
2016                 /* copy the neighbour lists */
2017 #ifdef STATIC_LISTS
2018                 /* no copy needed for static lists */
2019 #else
2020                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
2021 #endif
2022
2023                 /* second loop over atoms j */
2024                 for(j=0;j<27;j++) {
2025
2026                         bc_ij=(j<dnlc)?0:1;
2027 #ifdef STATIC_LISTS
2028                         p=0;
2029
2030                         while(neighbour_i[j][p]!=0) {
2031
2032                                 jtom=&(atom[neighbour_i[j][p]]);
2033                                 p++;
2034 #else
2035                         this=&(neighbour_i[j]);
2036                         list_reset_f(this);
2037
2038                         if(this->start==NULL)
2039                                 continue;
2040
2041                         do {
2042
2043                                 jtom=this->current->data;
2044 #endif
2045
2046                                 if(jtom==&(itom[i]))
2047                                         continue;
2048
2049                                 if(!(jtom->attr&ATOM_ATTR_3BP))
2050                                         continue;
2051
2052                                 /* reset 3bp run */
2053                                 moldyn->run3bp=1;
2054
2055                                 if(moldyn->func3b_j1)
2056                                         moldyn->func3b_j1(moldyn,
2057                                                           &(itom[i]),
2058                                                           jtom,
2059                                                           bc_ij);
2060
2061                                 /* in first j loop, 3bp run can be skipped */
2062                                 if(!(moldyn->run3bp))
2063                                         continue;
2064
2065                                 /* first loop over atoms k */
2066                                 if(moldyn->func3b_k1) {
2067
2068                                 for(k=0;k<27;k++) {
2069
2070                                         bc_ik=(k<dnlc)?0:1;
2071 #ifdef STATIC_LISTS
2072                                         q=0;
2073
2074                                         while(neighbour_i[j][q]!=0) {
2075
2076                                                 ktom=&(atom[neighbour_i[k][q]]);
2077                                                 q++;
2078 #else
2079                                         that=&(neighbour_i2[k]);
2080                                         list_reset_f(that);
2081
2082                                         if(that->start==NULL)
2083                                                 continue;
2084
2085                                         do {
2086                                                 ktom=that->current->data;
2087 #endif
2088
2089                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
2090                                                         continue;
2091
2092                                                 if(ktom==jtom)
2093                                                         continue;
2094
2095                                                 if(ktom==&(itom[i]))
2096                                                         continue;
2097
2098                                                 moldyn->func3b_k1(moldyn,
2099                                                                   &(itom[i]),
2100                                                                   jtom,
2101                                                                   ktom,
2102                                                                   bc_ik|bc_ij);
2103 #ifdef STATIC_LISTS
2104                                         }
2105 #else
2106                                         } while(list_next_f(that)!=\
2107                                                 L_NO_NEXT_ELEMENT);
2108 #endif
2109
2110                                 }
2111
2112                                 }
2113
2114                                 if(moldyn->func3b_j2)
2115                                         moldyn->func3b_j2(moldyn,
2116                                                           &(itom[i]),
2117                                                           jtom,
2118                                                           bc_ij);
2119
2120                                 /* second loop over atoms k */
2121                                 if(moldyn->func3b_k2) {
2122
2123                                 for(k=0;k<27;k++) {
2124
2125                                         bc_ik=(k<dnlc)?0:1;
2126 #ifdef STATIC_LISTS
2127                                         q=0;
2128
2129                                         while(neighbour_i[j][q]!=0) {
2130
2131                                                 ktom=&(atom[neighbour_i[k][q]]);
2132                                                 q++;
2133 #else
2134                                         that=&(neighbour_i2[k]);
2135                                         list_reset_f(that);
2136
2137                                         if(that->start==NULL)
2138                                                 continue;
2139
2140                                         do {
2141                                                 ktom=that->current->data;
2142 #endif
2143
2144                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
2145                                                         continue;
2146
2147                                                 if(ktom==jtom)
2148                                                         continue;
2149
2150                                                 if(ktom==&(itom[i]))
2151                                                         continue;
2152
2153                                                 moldyn->func3b_k2(moldyn,
2154                                                                   &(itom[i]),
2155                                                                   jtom,
2156                                                                   ktom,
2157                                                                   bc_ik|bc_ij);
2158
2159 #ifdef STATIC_LISTS
2160                                         }
2161 #else
2162                                         } while(list_next_f(that)!=\
2163                                                 L_NO_NEXT_ELEMENT);
2164 #endif
2165
2166                                 }
2167
2168                                 }
2169
2170                                 /* 2bp post function */
2171                                 if(moldyn->func3b_j3) {
2172                                         moldyn->func3b_j3(moldyn,
2173                                                           &(itom[i]),
2174                                                           jtom,bc_ij);
2175                                 }
2176 #ifdef STATIC_LISTS
2177                         }
2178 #else
2179                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2180 #endif
2181
2182                 }
2183
2184 #ifdef DEBUG
2185         //printf("\n\n");
2186 #endif
2187 #ifdef VDEBUG
2188         printf("\n\n");
2189 #endif
2190
2191         }
2192
2193 #ifdef DEBUG
2194         //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
2195         if(moldyn->time>DSTART&&moldyn->time<DEND) {
2196                 printf("force:\n");
2197                 printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
2198                 printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
2199                 printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
2200         }
2201 #endif
2202
2203         /* some postprocessing */
2204         for(i=0;i<count;i++) {
2205                 /* calculate global virial */
2206                 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
2207                 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
2208                 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
2209                 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
2210                 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
2211                 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
2212
2213                 /* check forces regarding the given timestep */
2214                 if(v3_norm(&(itom[i].f))>\
2215                    0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
2216                         printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
2217                                i);
2218         }
2219
2220         return 0;
2221 }
2222
2223 /*
2224  * virial calculation
2225  */
2226
2227 //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
2228 int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
2229
2230         a->virial.xx+=f->x*d->x;
2231         a->virial.yy+=f->y*d->y;
2232         a->virial.zz+=f->z*d->z;
2233         a->virial.xy+=f->x*d->y;
2234         a->virial.xz+=f->x*d->z;
2235         a->virial.yz+=f->y*d->z;
2236
2237         return 0;
2238 }
2239
2240 /*
2241  * periodic boundary checking
2242  */
2243
2244 //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
2245 int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
2246
2247         double x,y,z;
2248         t_3dvec *dim;
2249
2250         dim=&(moldyn->dim);
2251
2252         x=dim->x/2;
2253         y=dim->y/2;
2254         z=dim->z/2;
2255
2256         if(moldyn->status&MOLDYN_STAT_PBX) {
2257                 if(a->x>=x) a->x-=dim->x;
2258                 else if(-a->x>x) a->x+=dim->x;
2259         }
2260         if(moldyn->status&MOLDYN_STAT_PBY) {
2261                 if(a->y>=y) a->y-=dim->y;
2262                 else if(-a->y>y) a->y+=dim->y;
2263         }
2264         if(moldyn->status&MOLDYN_STAT_PBZ) {
2265                 if(a->z>=z) a->z-=dim->z;
2266                 else if(-a->z>z) a->z+=dim->z;
2267         }
2268
2269         return 0;
2270 }
2271
2272 /*
2273  * debugging / critical check functions
2274  */
2275
2276 int moldyn_bc_check(t_moldyn *moldyn) {
2277
2278         t_atom *atom;
2279         t_3dvec *dim;
2280         int i;
2281         double x;
2282         u8 byte;
2283         int j,k;
2284
2285         atom=moldyn->atom;
2286         dim=&(moldyn->dim);
2287         x=dim->x/2;
2288
2289         for(i=0;i<moldyn->count;i++) {
2290                 if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {
2291                         printf("FATAL: atom %d: x: %.20f (%.20f)\n",
2292                                i,atom[i].r.x,dim->x/2);
2293                         printf("diagnostic:\n");
2294                         printf("-----------\natom.r.x:\n");
2295                         for(j=0;j<8;j++) {
2296                                 memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
2297                                 for(k=0;k<8;k++)
2298                                         printf("%d%c",
2299                                         ((byte)&(1<<k))?1:0,
2300                                         (k==7)?'\n':'|');
2301                         }
2302                         printf("---------------\nx=dim.x/2:\n");
2303                         for(j=0;j<8;j++) {
2304                                 memcpy(&byte,(u8 *)(&x)+j,1);
2305                                 for(k=0;k<8;k++)
2306                                         printf("%d%c",
2307                                         ((byte)&(1<<k))?1:0,
2308                                         (k==7)?'\n':'|');
2309                         }
2310                         if(atom[i].r.x==x) printf("the same!\n");
2311                         else printf("different!\n");
2312                 }
2313                 if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
2314                         printf("FATAL: atom %d: y: %.20f (%.20f)\n",
2315                                i,atom[i].r.y,dim->y/2);
2316                 if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
2317                         printf("FATAL: atom %d: z: %.20f (%.20f)\n",
2318                                i,atom[i].r.z,dim->z/2);
2319         }
2320
2321         return 0;
2322 }
2323
2324 /*
2325  * restore function
2326  */
2327
2328 int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
2329
2330         int fd;
2331         int cnt,size;
2332         int fsize;
2333         int corr;
2334
2335         fd=open(file,O_RDONLY);
2336         if(fd<0) {
2337                 perror("[moldyn] load save file open");
2338                 return fd;
2339         }
2340
2341         fsize=lseek(fd,0,SEEK_END);
2342         lseek(fd,0,SEEK_SET);
2343
2344         size=sizeof(t_moldyn);
2345
2346         while(size) {
2347                 cnt=read(fd,moldyn,size);
2348                 if(cnt<0) {
2349                         perror("[moldyn] load save file read (moldyn)");
2350                         return cnt;
2351                 }
2352                 size-=cnt;
2353         }
2354
2355         size=moldyn->count*sizeof(t_atom);
2356
2357         /* correcting possible atom data offset */
2358         corr=0;
2359         if(fsize!=sizeof(t_moldyn)+size) {
2360                 corr=fsize-sizeof(t_moldyn)-size;
2361                 printf("[moldyn] WARNING: lsf (illegal file size)\n");
2362                 printf("  moifying offset:\n");
2363                 printf("  - current pos: %d\n",sizeof(t_moldyn));
2364                 printf("  - atom size: %d\n",size);
2365                 printf("  - file size: %d\n",fsize);
2366                 printf("  => correction: %d\n",corr);
2367                 lseek(fd,corr,SEEK_CUR);
2368         }
2369
2370         moldyn->atom=(t_atom *)malloc(size);
2371         if(moldyn->atom==NULL) {
2372                 perror("[moldyn] load save file malloc (atoms)");
2373                 return -1;
2374         }
2375
2376         while(size) {
2377                 cnt=read(fd,moldyn->atom,size);
2378                 if(cnt<0) {
2379                         perror("[moldyn] load save file read (atoms)");
2380                         return cnt;
2381                 }
2382                 size-=cnt;
2383         }
2384
2385         // hooks etc ...
2386
2387         return 0;
2388 }
2389
2390 int moldyn_free_save_file(t_moldyn *moldyn) {
2391
2392         free(moldyn->atom);
2393
2394         return 0;
2395 }
2396
2397 int moldyn_load(t_moldyn *moldyn) {
2398
2399         // later ...
2400
2401         return 0;
2402 }
2403
2404 /*
2405  * function to find/callback all combinations of 2 body bonds
2406  */
2407
2408 int process_2b_bonds(t_moldyn *moldyn,void *data,
2409                      int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
2410                                     void *data,u8 bc)) {
2411
2412         t_linkcell *lc;
2413 #ifdef STATIC_LISTS
2414         int *neighbour[27];
2415         int p;
2416 #else
2417         t_list neighbour[27];
2418 #endif
2419         u8 bc;
2420         t_atom *itom,*jtom;
2421         int i,j;
2422         t_list *this;
2423
2424         lc=&(moldyn->lc);
2425         itom=moldyn->atom;
2426
2427         for(i=0;i<moldyn->count;i++) {
2428                 /* neighbour indexing */
2429                 link_cell_neighbour_index(moldyn,
2430                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
2431                                           (itom[i].r.y+moldyn->dim.y/2)/lc->x,
2432                                           (itom[i].r.z+moldyn->dim.z/2)/lc->x,
2433                                           neighbour);
2434
2435                 for(j=0;j<27;j++) {
2436
2437                         bc=(j<lc->dnlc)?0:1;
2438
2439 #ifdef STATIC_LISTS
2440                         p=0;
2441
2442                         while(neighbour[j][p]!=0) {
2443
2444                                 jtom=&(moldyn->atom[neighbour[j][p]]);
2445                                 p++;
2446 #else
2447                         this=&(neighbour[j]);
2448                         list_reset_f(this);
2449
2450                         if(this->start==NULL)
2451                                 continue;
2452
2453                         do {
2454
2455                                 jtom=this->current->data;
2456 #endif
2457
2458                                 /* process bond */
2459                                 process(moldyn,&(itom[i]),jtom,data,bc);
2460
2461 #ifdef STATIC_LISTS
2462                         }
2463 #else
2464                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2465 #endif
2466                 }
2467         }
2468
2469         return 0;
2470
2471 }
2472
2473 /*
2474  * post processing functions
2475  */
2476
2477 int get_line(int fd,char *line,int max) {
2478
2479         int count,ret;
2480
2481         count=0;
2482
2483         while(1) {
2484                 if(count==max) return count;
2485                 ret=read(fd,line+count,1);
2486                 if(ret<=0) return ret;
2487                 if(line[count]=='\n') {
2488                         memset(line+count,0,max-count-1);
2489                         //line[count]='\0';
2490                         return count+1;
2491                 }
2492                 count+=1;
2493         }
2494 }
2495
2496 int pair_correlation_init(t_moldyn *moldyn,double dr) {
2497
2498
2499         return 0;
2500 }
2501
2502 int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
2503
2504         int i;
2505         t_atom *atom;
2506         t_3dvec dist;
2507         double d2;
2508         int a_cnt;
2509         int b_cnt;
2510
2511         atom=moldyn->atom;
2512         dc[0]=0;
2513         dc[1]=0;
2514         dc[2]=0;
2515         a_cnt=0;
2516         b_cnt=0;
2517
2518         for(i=0;i<moldyn->count;i++) {
2519
2520                 v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
2521                 check_per_bound(moldyn,&dist);
2522                 d2=v3_absolute_square(&dist);
2523
2524                 if(atom[i].brand) {
2525                         b_cnt+=1;
2526                         dc[1]+=d2;
2527                 }
2528                 else {
2529                         a_cnt+=1;
2530                         dc[0]+=d2;
2531                 }
2532
2533                 dc[2]+=d2;
2534         }
2535
2536         dc[0]*=(1.0/(6.0*moldyn->time*a_cnt));
2537         dc[1]*=(1.0/(6.0*moldyn->time*b_cnt));
2538         dc[2]*=(1.0/(6.0*moldyn->time*moldyn->count));
2539
2540         return 0;
2541 }
2542
2543 int bonding_analyze(t_moldyn *moldyn,double *cnt) {
2544
2545         return 0;
2546 }
2547
2548 int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
2549                                        t_atom *jtom,void *data,u8 bc) {
2550
2551         t_3dvec dist;
2552         double d;
2553         int s;
2554         t_pcc *pcc;
2555
2556         /* only count pairs once,
2557          * skip same atoms */
2558         if(itom->tag>=jtom->tag)
2559                 return 0;
2560
2561         /*
2562          * pair correlation calc
2563          */
2564
2565         /* get pcc data */
2566         pcc=data;
2567
2568         /* distance */
2569         v3_sub(&dist,&(jtom->r),&(itom->r));
2570         if(bc) check_per_bound(moldyn,&dist);
2571         d=v3_absolute_square(&dist);
2572
2573         /* ignore if greater cutoff */
2574         if(d>moldyn->cutoff_square)
2575                 return 0;
2576
2577         /* fill the slots */
2578         d=sqrt(d);
2579         s=(int)(d/pcc->dr);
2580
2581         /* should never happen but it does 8) -
2582          * related to -ffloat-store problem! */
2583         if(s>=pcc->o1) {
2584                 printf("[moldyn] WARNING: pcc (%d/%d)",
2585                        s,pcc->o1);
2586                 printf("\n");
2587                 s=pcc->o1-1;
2588         }
2589
2590         if(itom->brand!=jtom->brand) {
2591                 /* mixed */
2592                 pcc->stat[s]+=1;
2593         }
2594         else {
2595                 /* type a - type a bonds */
2596                 if(itom->brand==0)
2597                         pcc->stat[s+pcc->o1]+=1;
2598                 else
2599                 /* type b - type b bonds */
2600                         pcc->stat[s+pcc->o2]+=1;
2601         }
2602
2603         return 0;
2604 }
2605
2606 int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
2607
2608         t_pcc pcc;
2609         double norm;
2610         int i;
2611
2612         pcc.dr=dr;
2613         pcc.o1=moldyn->cutoff/dr;
2614         pcc.o2=2*pcc.o1;
2615
2616         if(pcc.o1*dr<=moldyn->cutoff)
2617                 printf("[moldyn] WARNING: pcc (low #slots)\n");
2618
2619         printf("[moldyn] pair correlation calc info:\n");
2620         printf("  time: %f\n",moldyn->time);
2621         printf("  count: %d\n",moldyn->count);
2622         printf("  cutoff: %f\n",moldyn->cutoff);
2623         printf("  temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg);
2624
2625         if(ptr!=NULL) {
2626                 pcc.stat=(double *)ptr;
2627         }
2628         else {
2629                 pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double));
2630                 if(pcc.stat==NULL) {
2631                         perror("[moldyn] pair correlation malloc");
2632                         return -1;
2633                 }
2634         }
2635
2636         memset(pcc.stat,0,3*pcc.o1*sizeof(double));
2637
2638         /* process */
2639         process_2b_bonds(moldyn,&pcc,calculate_pair_correlation_process);
2640
2641         /* normalization */
2642         for(i=1;i<pcc.o1;i++) {
2643                  // normalization: 4 pi r^2 dr
2644                  // here: not double counting pairs -> 2 pi r r dr
2645                  // ... and actually it's a constant times r^2
2646                 norm=i*i*dr*dr;
2647                 pcc.stat[i]/=norm;
2648                 pcc.stat[pcc.o1+i]/=norm;
2649                 pcc.stat[pcc.o2+i]/=norm;
2650         }
2651         /* */
2652
2653         if(ptr==NULL) {
2654                 /* todo: store/print pair correlation function */
2655                 free(pcc.stat);
2656         }
2657
2658         return 0;
2659 }
2660
2661 int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
2662                          void *data,u8 bc) {
2663
2664         t_ba *ba;
2665         t_3dvec dist;
2666         double d;
2667
2668         if(itom->tag>=jtom->tag)
2669                 return 0;
2670
2671         /* distance */
2672         v3_sub(&dist,&(jtom->r),&(itom->r));
2673         if(bc) check_per_bound(moldyn,&dist);
2674         d=v3_absolute_square(&dist);
2675
2676         /* ignore if greater or equal cutoff */
2677         if(d>moldyn->cutoff_square)
2678                 return 0;
2679
2680         /* check for potential bond */
2681         if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
2682                 return 0;
2683
2684         d=sqrt(d);
2685
2686         /* now count this bonding ... */
2687         ba=data;
2688
2689         /* increase total bond counter
2690          * ... double counting!
2691          */
2692         ba->tcnt+=2;
2693
2694         if(itom->brand==0)
2695                 ba->acnt[jtom->tag]+=1;
2696         else
2697                 ba->bcnt[jtom->tag]+=1;
2698
2699         if(jtom->brand==0)
2700                 ba->acnt[itom->tag]+=1;
2701         else
2702                 ba->bcnt[itom->tag]+=1;
2703
2704         return 0;
2705 }
2706
2707 int bond_analyze(t_moldyn *moldyn,double *quality) {
2708
2709         // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
2710
2711         int qcnt;
2712         int ccnt,cset;
2713         t_ba ba;
2714         int i;
2715         t_atom *atom;
2716
2717         ba.acnt=malloc(moldyn->count*sizeof(int));
2718         if(ba.acnt==NULL) {
2719                 perror("[moldyn] bond analyze malloc (a)");
2720                 return -1;
2721         }
2722         memset(ba.acnt,0,moldyn->count*sizeof(int));
2723
2724         ba.bcnt=malloc(moldyn->count*sizeof(int));
2725         if(ba.bcnt==NULL) {
2726                 perror("[moldyn] bond analyze malloc (b)");
2727                 return -1;
2728         }
2729         memset(ba.bcnt,0,moldyn->count*sizeof(int));
2730
2731         ba.tcnt=0;
2732         qcnt=0;
2733         ccnt=0;
2734         cset=0;
2735
2736         atom=moldyn->atom;
2737
2738         process_2b_bonds(moldyn,&ba,bond_analyze_process);
2739
2740         for(i=0;i<moldyn->count;i++) {
2741                 if(atom[i].brand==0) {
2742                         if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
2743                                 qcnt+=4;
2744                 }
2745                 else {
2746                         if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
2747                                 qcnt+=4;
2748                                 ccnt+=1;
2749                         }
2750                         cset+=1;
2751                 }
2752         }
2753
2754         printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
2755         printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
2756
2757         if(quality) {
2758                 quality[0]=1.0*ccnt/cset;
2759                 quality[1]=1.0*qcnt/ba.tcnt;
2760         }
2761         else {
2762                 printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
2763                 printf("[moldyn] bond analyze:   tot_q=%f\n",1.0*qcnt/ba.tcnt);
2764         }
2765
2766         return 0;
2767 }
2768
2769 /*
2770  * visualization code
2771  */
2772
2773 int visual_init(t_moldyn *moldyn,char *filebase) {
2774
2775         strncpy(moldyn->vis.fb,filebase,128);
2776
2777         return 0;
2778 }
2779
2780 int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
2781                          void *data,u8 bc) {
2782
2783         t_vb *vb;
2784
2785         vb=data;
2786
2787         if(itom->tag>=jtom->tag)
2788                 return 0;
2789
2790         if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
2791                 return 0;
2792
2793         if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB))
2794                 dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n",
2795                         itom->r.x,itom->r.y,itom->r.z,
2796                         jtom->r.x,jtom->r.y,jtom->r.z);
2797
2798         return 0;
2799 }
2800
2801 int visual_atoms(t_moldyn *moldyn) {
2802
2803         int i;
2804         char file[128+64];
2805         t_3dvec dim;
2806         double help;
2807         t_visual *v;
2808         t_atom *atom;
2809         t_vb vb;
2810
2811         v=&(moldyn->vis);
2812         dim.x=v->dim.x;
2813         dim.y=v->dim.y;
2814         dim.z=v->dim.z;
2815         atom=moldyn->atom;
2816
2817         help=(dim.x+dim.y);
2818
2819         sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
2820         vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
2821         if(vb.fd<0) {
2822                 perror("open visual save file fd");
2823                 return -1;
2824         }
2825
2826         /* write the actual data file */
2827
2828         // povray header
2829         dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n",
2830                 moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
2831
2832         // atomic configuration
2833         for(i=0;i<moldyn->count;i++)
2834                 // atom type, positions, color and kinetic energy
2835                 dprintf(vb.fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
2836                                                     atom[i].r.x,
2837                                                     atom[i].r.y,
2838                                                     atom[i].r.z,
2839                                                     pse_col[atom[i].element],
2840                                                     atom[i].ekin);
2841
2842         // bonds between atoms
2843         process_2b_bonds(moldyn,&vb,visual_bonds_process);
2844
2845         // boundaries
2846         if(dim.x) {
2847                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2848                         -dim.x/2,-dim.y/2,-dim.z/2,
2849                         dim.x/2,-dim.y/2,-dim.z/2);
2850                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2851                         -dim.x/2,-dim.y/2,-dim.z/2,
2852                         -dim.x/2,dim.y/2,-dim.z/2);
2853                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2854                         dim.x/2,dim.y/2,-dim.z/2,
2855                         dim.x/2,-dim.y/2,-dim.z/2);
2856                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2857                         -dim.x/2,dim.y/2,-dim.z/2,
2858                         dim.x/2,dim.y/2,-dim.z/2);
2859
2860                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2861                         -dim.x/2,-dim.y/2,dim.z/2,
2862                         dim.x/2,-dim.y/2,dim.z/2);
2863                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2864                         -dim.x/2,-dim.y/2,dim.z/2,
2865                         -dim.x/2,dim.y/2,dim.z/2);
2866                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2867                         dim.x/2,dim.y/2,dim.z/2,
2868                         dim.x/2,-dim.y/2,dim.z/2);
2869                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2870                         -dim.x/2,dim.y/2,dim.z/2,
2871                         dim.x/2,dim.y/2,dim.z/2);
2872
2873                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2874                         -dim.x/2,-dim.y/2,dim.z/2,
2875                         -dim.x/2,-dim.y/2,-dim.z/2);
2876                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2877                         -dim.x/2,dim.y/2,dim.z/2,
2878                         -dim.x/2,dim.y/2,-dim.z/2);
2879                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2880                         dim.x/2,-dim.y/2,dim.z/2,
2881                         dim.x/2,-dim.y/2,-dim.z/2);
2882                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2883                         dim.x/2,dim.y/2,dim.z/2,
2884                         dim.x/2,dim.y/2,-dim.z/2);
2885         }
2886
2887         close(vb.fd);
2888
2889         return 0;
2890 }
2891