moldyn.c

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