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