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