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