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