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