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 #ifdef LOWMEM_LISTS
 604         ptr=realloc(moldyn->lc.subcell->list,(count+1)*sizeof(int));
 605         if(!ptr) {
 606                 perror("[moldyn] list realloc (add atom)");
 607                 return -1;
 608         }
 609         moldyn->lc.subcell->list=ptr;
 610 #endif
 611
 612         atom=moldyn->atom;
 613
 614         /* initialize new atom */
 615         memset(&(atom[count]),0,sizeof(t_atom));
 616         atom[count].r=*r;
 617         atom[count].v=*v;
 618         atom[count].element=element;
 619         atom[count].mass=mass;
 620         atom[count].brand=brand;
 621         atom[count].tag=count;
 622         atom[count].attr=attr;
 623         check_per_bound(moldyn,&(atom[count].r));
 624         atom[count].r_0=atom[count].r;
 625
 626         /* update total system mass */
 627         total_mass_calc(moldyn);
 628
 629         return 0;
 630 }
 631
 632 int del_atom(t_moldyn *moldyn,int tag) {
 633
 634         t_atom *new,*old;
 635         int cnt;
 636
 637         old=moldyn->atom;
 638
 639         new=(t_atom *)malloc((moldyn->count-1)*sizeof(t_atom));
 640         if(!new) {
 641                 perror("[moldyn]malloc (del atom)");
 642                 return -1;
 643         }
 644
 645         for(cnt=0;cnt<tag;cnt++)
 646                 new[cnt]=old[cnt];
 647
 648         for(cnt=tag+1;cnt<moldyn->count;cnt++) {
 649                 new[cnt-1]=old[cnt];
 650                 new[cnt-1].tag=cnt-1;
 651         }
 652
 653         moldyn->count-=1;
 654         moldyn->atom=new;
 655
 656         free(old);
 657
 658         return 0;
 659 }
 660
 661 /* cubic init */
 662 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 663
 664         int count;
 665         t_3dvec r;
 666         int i,j,k;
 667         t_3dvec o;
 668
 669         count=0;
 670         if(origin)
 671                 v3_copy(&o,origin);
 672         else
 673                 v3_zero(&o);
 674
 675         r.x=o.x;
 676         for(i=0;i<a;i++) {
 677                 r.y=o.y;
 678                 for(j=0;j<b;j++) {
 679                         r.z=o.z;
 680                         for(k=0;k<c;k++) {
 681                                 v3_copy(&(atom[count].r),&r);
 682                                 count+=1;
 683                                 r.z+=lc;
 684                         }
 685                         r.y+=lc;
 686                 }
 687                 r.x+=lc;
 688         }
 689
 690         for(i=0;i<count;i++) {
 691                 atom[i].r.x-=(a*lc)/2.0;
 692                 atom[i].r.y-=(b*lc)/2.0;
 693                 atom[i].r.z-=(c*lc)/2.0;
 694         }
 695
 696         return count;
 697 }
 698
 699 /* fcc lattice init */
 700 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 701
 702         int count;
 703         int i,j,k,l;
 704         t_3dvec o,r,n;
 705         t_3dvec basis[3];
 706
 707         count=0;
 708         if(origin)
 709                 v3_copy(&o,origin);
 710         else
 711                 v3_zero(&o);
 712
 713         /* construct the basis */
 714         memset(basis,0,3*sizeof(t_3dvec));
 715         basis[0].x=0.5*lc;
 716         basis[0].y=0.5*lc;
 717         basis[1].x=0.5*lc;
 718         basis[1].z=0.5*lc;
 719         basis[2].y=0.5*lc;
 720         basis[2].z=0.5*lc;
 721
 722         /* fill up the room */
 723         r.x=o.x;
 724         for(i=0;i<a;i++) {
 725                 r.y=o.y;
 726                 for(j=0;j<b;j++) {
 727                         r.z=o.z;
 728                         for(k=0;k<c;k++) {
 729                                 /* first atom */
 730                                 v3_copy(&(atom[count].r),&r);
 731                                 count+=1;
 732                                 r.z+=lc;
 733                                 /* the three face centered atoms */
 734                                 for(l=0;l<3;l++) {
 735                                         v3_add(&n,&r,&basis[l]);
 736                                         v3_copy(&(atom[count].r),&n);
 737                                         count+=1;
 738                                 }
 739                         }
 740                         r.y+=lc;
 741                 }
 742                 r.x+=lc;
 743         }
 744
 745         /* coordinate transformation */
 746         for(i=0;i<count;i++) {
 747                 atom[i].r.x-=(a*lc)/2.0;
 748                 atom[i].r.y-=(b*lc)/2.0;
 749                 atom[i].r.z-=(c*lc)/2.0;
 750         }
 751
 752         return count;
 753 }
 754
 755 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 756
 757         int count;
 758         t_3dvec o;
 759
 760         count=fcc_init(a,b,c,lc,atom,origin);
 761
 762         o.x=0.25*lc;
 763         o.y=0.25*lc;
 764         o.z=0.25*lc;
 765
 766         if(origin) v3_add(&o,&o,origin);
 767
 768         count+=fcc_init(a,b,c,lc,&atom[count],&o);
 769
 770         return count;
 771 }
 772
 773 int destroy_atoms(t_moldyn *moldyn) {
 774
 775         if(moldyn->atom) free(moldyn->atom);
 776
 777         return 0;
 778 }
 779
 780 int thermal_init(t_moldyn *moldyn,u8 equi_init) {
 781
 782         /*
 783          * - gaussian distribution of velocities
 784          * - zero total momentum
 785          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 786          */
 787
 788         int i;
 789         double v,sigma;
 790         t_3dvec p_total,delta;
 791         t_atom *atom;
 792         t_random *random;
 793
 794         atom=moldyn->atom;
 795         random=&(moldyn->random);
 796
 797         printf("[moldyn] thermal init (equi init: %s)\n",equi_init?"yes":"no");
 798
 799         /* gaussian distribution of velocities */
 800         v3_zero(&p_total);
 801         for(i=0;i<moldyn->count;i++) {
 802                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
 803                 /* x direction */
 804                 v=sigma*rand_get_gauss(random);
 805                 atom[i].v.x=v;
 806                 p_total.x+=atom[i].mass*v;
 807                 /* y direction */
 808                 v=sigma*rand_get_gauss(random);
 809                 atom[i].v.y=v;
 810                 p_total.y+=atom[i].mass*v;
 811                 /* z direction */
 812                 v=sigma*rand_get_gauss(random);
 813                 atom[i].v.z=v;
 814                 p_total.z+=atom[i].mass*v;
 815         }
 816
 817         /* zero total momentum */
 818         v3_scale(&p_total,&p_total,1.0/moldyn->count);
 819         for(i=0;i<moldyn->count;i++) {
 820                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 821                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 822         }
 823
 824         /* velocity scaling */
 825         scale_velocity(moldyn,equi_init);
 826
 827         return 0;
 828 }
 829
 830 double total_mass_calc(t_moldyn *moldyn) {
 831
 832         int i;
 833
 834         moldyn->mass=0.0;
 835
 836         for(i=0;i<moldyn->count;i++)
 837                 moldyn->mass+=moldyn->atom[i].mass;
 838
 839         return moldyn->mass;
 840 }
 841
 842 double temperature_calc(t_moldyn *moldyn) {
 843
 844         /* assume up to date kinetic energy, which is 3/2 N k_B T */
 845
 846         moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
 847
 848         return moldyn->t;
 849 }
 850
 851 double get_temperature(t_moldyn *moldyn) {
 852
 853         return moldyn->t;
 854 }
 855
 856 int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 857
 858         int i;
 859         double e,scale;
 860         t_atom *atom;
 861         int count;
 862
 863         atom=moldyn->atom;
 864
 865         /*
 866          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 867          */
 868
 869         /* get kinetic energy / temperature & count involved atoms */
 870         e=0.0;
 871         count=0;
 872         for(i=0;i<moldyn->count;i++) {
 873                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) {
 874                         e+=atom[i].mass*v3_absolute_square(&(atom[i].v));
 875                         count+=1;
 876                 }
 877         }
 878         e*=0.5;
 879         if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN);
 880         else return 0;  /* no atoms involved in scaling! */
 881
 882         /* (temporary) hack for e,t = 0 */
 883         if(e==0.0) {
 884         moldyn->t=0.0;
 885                 if(moldyn->t_ref!=0.0) {
 886                         thermal_init(moldyn,equi_init);
 887                         return 0;
 888                 }
 889                 else
 890                         return 0; /* no scaling needed */
 891         }
 892
 893
 894         /* get scaling factor */
 895         scale=moldyn->t_ref/moldyn->t;
 896         if(equi_init&TRUE)
 897                 scale*=2.0;
 898         else
 899                 if(moldyn->pt_scale&T_SCALE_BERENDSEN)
 900                         scale=1.0+(scale-1.0)*moldyn->tau/moldyn->t_tc;
 901         scale=sqrt(scale);
 902
 903         /* velocity scaling */
 904         for(i=0;i<moldyn->count;i++) {
 905                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
 906                         v3_scale(&(atom[i].v),&(atom[i].v),scale);
 907         }
 908
 909         return 0;
 910 }
 911
 912 double ideal_gas_law_pressure(t_moldyn *moldyn) {
 913
 914         double p;
 915
 916         p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume;
 917
 918         return p;
 919 }
 920
 921 double virial_sum(t_moldyn *moldyn) {
 922
 923         int i;
 924         t_virial *virial;
 925
 926         /* virial (sum over atom virials) */
 927         moldyn->virial=0.0;
 928         moldyn->vir.xx=0.0;
 929         moldyn->vir.yy=0.0;
 930         moldyn->vir.zz=0.0;
 931         moldyn->vir.xy=0.0;
 932         moldyn->vir.xz=0.0;
 933         moldyn->vir.yz=0.0;
 934         for(i=0;i<moldyn->count;i++) {
 935                 virial=&(moldyn->atom[i].virial);
 936                 moldyn->virial+=(virial->xx+virial->yy+virial->zz);
 937                 moldyn->vir.xx+=virial->xx;
 938                 moldyn->vir.yy+=virial->yy;
 939                 moldyn->vir.zz+=virial->zz;
 940                 moldyn->vir.xy+=virial->xy;
 941                 moldyn->vir.xz+=virial->xz;
 942                 moldyn->vir.yz+=virial->yz;
 943         }
 944
 945         /* global virial (absolute coordinates) */
 946         virial=&(moldyn->gvir);
 947         moldyn->gv=virial->xx+virial->yy+virial->zz;
 948
 949         return moldyn->virial;
 950 }
 951
 952 double pressure_calc(t_moldyn *moldyn) {
 953
 954         /*
 955          * PV = NkT + <W>
 956          * with W = 1/3 sum_i f_i r_i (- skipped!)
 957          * virial = sum_i f_i r_i
 958          *
 959          * => P = (2 Ekin + virial) / (3V)
 960          */
 961
 962         /* assume up to date virial & up to date kinetic energy */
 963
 964         /* pressure (atom virials) */
 965         moldyn->p=2.0*moldyn->ekin+moldyn->virial;
 966         moldyn->p/=(3.0*moldyn->volume);
 967
 968         /* pressure (absolute coordinates) */
 969         moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
 970         moldyn->gp/=(3.0*moldyn->volume);
 971
 972         return moldyn->p;
 973 }
 974
 975 int average_reset(t_moldyn *moldyn) {
 976
 977         printf("[moldyn] average reset\n");
 978
 979         /* update skip value */
 980         moldyn->avg_skip=moldyn->total_steps;
 981
 982         /* kinetic energy */
 983         moldyn->k_sum=0.0;
 984         moldyn->k2_sum=0.0;
 985
 986         /* potential energy */
 987         moldyn->v_sum=0.0;
 988         moldyn->v2_sum=0.0;
 989
 990         /* temperature */
 991         moldyn->t_sum=0.0;
 992
 993         /* virial */
 994         moldyn->virial_sum=0.0;
 995         moldyn->gv_sum=0.0;
 996
 997         /* pressure */
 998         moldyn->p_sum=0.0;
 999         moldyn->gp_sum=0.0;
1000         moldyn->tp_sum=0.0;
1001
1002         return 0;
1003 }
1004
1005 int average_and_fluctuation_calc(t_moldyn *moldyn) {
1006
1007         int denom;
1008
1009         if(moldyn->total_steps<moldyn->avg_skip)
1010                 return 0;
1011
1012         denom=moldyn->total_steps+1-moldyn->avg_skip;
1013
1014         /* assume up to date energies, temperature, pressure etc */
1015
1016         /* kinetic energy */
1017         moldyn->k_sum+=moldyn->ekin;
1018         moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
1019         moldyn->k_avg=moldyn->k_sum/denom;
1020         moldyn->k2_avg=moldyn->k2_sum/denom;
1021         moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg);
1022
1023         /* potential energy */
1024         moldyn->v_sum+=moldyn->energy;
1025         moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
1026         moldyn->v_avg=moldyn->v_sum/denom;
1027         moldyn->v2_avg=moldyn->v2_sum/denom;
1028         moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg);
1029
1030         /* temperature */
1031         moldyn->t_sum+=moldyn->t;
1032         moldyn->t_avg=moldyn->t_sum/denom;
1033
1034         /* virial */
1035         moldyn->virial_sum+=moldyn->virial;
1036         moldyn->virial_avg=moldyn->virial_sum/denom;
1037         moldyn->gv_sum+=moldyn->gv;
1038         moldyn->gv_avg=moldyn->gv_sum/denom;
1039
1040         /* pressure */
1041         moldyn->p_sum+=moldyn->p;
1042         moldyn->p_avg=moldyn->p_sum/denom;
1043         moldyn->gp_sum+=moldyn->gp;
1044         moldyn->gp_avg=moldyn->gp_sum/denom;
1045         moldyn->tp_sum+=moldyn->tp;
1046         moldyn->tp_avg=moldyn->tp_sum/denom;
1047
1048         return 0;
1049 }
1050
1051 int get_heat_capacity(t_moldyn *moldyn) {
1052
1053         double temp2,ighc;
1054
1055         /* averages needed for heat capacity calc */
1056         if(moldyn->total_steps<moldyn->avg_skip)
1057                 return 0;
1058
1059         /* (temperature average)^2 */
1060         temp2=moldyn->t_avg*moldyn->t_avg;
1061         printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
1062                moldyn->t_avg);
1063
1064         /* ideal gas contribution */
1065         ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
1066         printf("  ideal gas contribution: %f\n",
1067                ighc/moldyn->mass*KILOGRAM/JOULE);
1068
1069         /* specific heat for nvt ensemble */
1070         moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc;
1071         moldyn->c_v_nvt/=moldyn->mass;
1072
1073         /* specific heat for nve ensemble */
1074         moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2)));
1075         moldyn->c_v_nve/=moldyn->mass;
1076
1077         printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
1078         printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
1079 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)));
1080
1081         return 0;
1082 }
1083
1084 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
1085
1086         t_3dvec dim;
1087         //t_3dvec *tp;
1088         double h,dv;
1089         double y0,y1;
1090         double su,sd;
1091         t_atom *store;
1092
1093         /*
1094          * dU = - p dV
1095          *
1096          * => p = - dU/dV
1097          *
1098          */
1099
1100         /* store atomic configuration + dimension */
1101         store=malloc(moldyn->count*sizeof(t_atom));
1102         if(store==NULL) {
1103                 printf("[moldyn] allocating store mem failed\n");
1104                 return -1;
1105         }
1106         memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
1107         dim=moldyn->dim;
1108
1109         /* x1, y1 */
1110         sd=0.00001;
1111         h=(1.0-sd)*(1.0-sd)*(1.0-sd);
1112         su=pow(2.0-h,ONE_THIRD)-1.0;
1113         dv=(1.0-h)*moldyn->volume;
1114
1115         /* scale up dimension and atom positions */
1116         scale_dim(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
1117         scale_atoms(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
1118         link_cell_shutdown(moldyn);
1119         link_cell_init(moldyn,QUIET);
1120         potential_force_calc(moldyn);
1121         y1=moldyn->energy;
1122
1123         /* restore atomic configuration + dim */
1124         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
1125         moldyn->dim=dim;
1126
1127         /* scale down dimension and atom positions */
1128         scale_dim(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
1129         scale_atoms(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
1130         link_cell_shutdown(moldyn);
1131         link_cell_init(moldyn,QUIET);
1132         potential_force_calc(moldyn);
1133         y0=moldyn->energy;
1134
1135         /* calculate pressure */
1136         moldyn->tp=-(y1-y0)/(2.0*dv);
1137
1138         /* restore atomic configuration */
1139         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
1140         moldyn->dim=dim;
1141         link_cell_shutdown(moldyn);
1142         link_cell_init(moldyn,QUIET);
1143         //potential_force_calc(moldyn);
1144
1145         /* free store buffer */
1146         if(store)
1147                 free(store);
1148
1149         return moldyn->tp;
1150 }
1151
1152 double get_pressure(t_moldyn *moldyn) {
1153
1154         return moldyn->p;
1155
1156 }
1157
1158 int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1159
1160         t_3dvec *dim;
1161
1162         dim=&(moldyn->dim);
1163
1164         if(dir==SCALE_UP)
1165                 scale=1.0+scale;
1166
1167         if(dir==SCALE_DOWN)
1168                 scale=1.0-scale;
1169
1170         if(x) dim->x*=scale;
1171         if(y) dim->y*=scale;
1172         if(z) dim->z*=scale;
1173
1174         return 0;
1175 }
1176
1177 int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1178
1179         int i;
1180         t_3dvec *r;
1181
1182         if(dir==SCALE_UP)
1183                 scale=1.0+scale;
1184
1185         if(dir==SCALE_DOWN)
1186                 scale=1.0-scale;
1187
1188         for(i=0;i<moldyn->count;i++) {
1189                 r=&(moldyn->atom[i].r);
1190                 if(x) r->x*=scale;
1191                 if(y) r->y*=scale;
1192                 if(z) r->z*=scale;
1193         }
1194
1195         return 0;
1196 }
1197
1198 int scale_volume(t_moldyn *moldyn) {
1199
1200         t_3dvec *dim,*vdim;
1201         double scale;
1202         t_linkcell *lc;
1203
1204         vdim=&(moldyn->vis.dim);
1205         dim=&(moldyn->dim);
1206         lc=&(moldyn->lc);
1207
1208         /* scaling factor */
1209         if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
1210                 scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau;
1211                 scale=pow(scale,ONE_THIRD);
1212         }
1213         else {
1214                 scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
1215         }
1216
1217         /* scale the atoms and dimensions */
1218         scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1219         scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1220
1221         /* visualize dimensions */
1222         if(vdim->x!=0) {
1223                 vdim->x=dim->x;
1224                 vdim->y=dim->y;
1225                 vdim->z=dim->z;
1226         }
1227
1228         /* recalculate scaled volume */
1229         moldyn->volume=dim->x*dim->y*dim->z;
1230
1231         /* adjust/reinit linkcell */
1232         if(((int)(dim->x/moldyn->cutoff)!=lc->nx)||
1233            ((int)(dim->y/moldyn->cutoff)!=lc->ny)||
1234            ((int)(dim->z/moldyn->cutoff)!=lc->nx)) {
1235                 link_cell_shutdown(moldyn);
1236                 link_cell_init(moldyn,QUIET);
1237         } else {
1238                 lc->x*=scale;
1239                 lc->y*=scale;
1240                 lc->z*=scale;
1241         }
1242
1243         return 0;
1244
1245 }
1246
1247 double e_kin_calc(t_moldyn *moldyn) {
1248
1249         int i;
1250         t_atom *atom;
1251
1252         atom=moldyn->atom;
1253         moldyn->ekin=0.0;
1254
1255         for(i=0;i<moldyn->count;i++) {
1256                 atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
1257                 moldyn->ekin+=atom[i].ekin;
1258         }
1259
1260         return moldyn->ekin;
1261 }
1262
1263 double get_total_energy(t_moldyn *moldyn) {
1264
1265         return(moldyn->ekin+moldyn->energy);
1266 }
1267
1268 t_3dvec get_total_p(t_moldyn *moldyn) {
1269
1270         t_3dvec p,p_total;
1271         int i;
1272         t_atom *atom;
1273
1274         atom=moldyn->atom;
1275
1276         v3_zero(&p_total);
1277         for(i=0;i<moldyn->count;i++) {
1278                 v3_scale(&p,&(atom[i].v),atom[i].mass);
1279                 v3_add(&p_total,&p_total,&p);
1280         }
1281
1282         return p_total;
1283 }
1284
1285 double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
1286
1287         double tau;
1288
1289         /* nn_dist is the nearest neighbour distance */
1290
1291         tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
1292
1293         return tau;
1294 }
1295
1296 /*
1297  * numerical tricks
1298  */
1299
1300 /* linked list / cell method */
1301
1302 int link_cell_init(t_moldyn *moldyn,u8 vol) {
1303
1304         t_linkcell *lc;
1305 #ifndef LOWMEM_LISTS
1306         int i;
1307 #endif
1308
1309         lc=&(moldyn->lc);
1310
1311         /* partitioning the md cell */
1312         lc->nx=moldyn->dim.x/moldyn->cutoff;
1313         lc->x=moldyn->dim.x/lc->nx;
1314         lc->ny=moldyn->dim.y/moldyn->cutoff;
1315         lc->y=moldyn->dim.y/lc->ny;
1316         lc->nz=moldyn->dim.z/moldyn->cutoff;
1317         lc->z=moldyn->dim.z/lc->nz;
1318         lc->cells=lc->nx*lc->ny*lc->nz;
1319
1320 #ifdef STATIC_LISTS
1321         lc->subcell=malloc(lc->cells*sizeof(int*));
1322 #elif LOWMEM_LISTS
1323         lc->subcell=malloc(sizeof(t_lowmem_list));
1324 #else
1325         lc->subcell=malloc(lc->cells*sizeof(t_list));
1326 #endif
1327
1328         if(lc->subcell==NULL) {
1329                 perror("[moldyn] cell init (malloc)");
1330                 return -1;
1331         }
1332
1333         if(lc->cells<27)
1334                 printf("[moldyn] FATAL: less then 27 subcells! (%d)\n",
1335                        lc->cells);
1336
1337         if(vol) {
1338 #ifdef STATIC_LISTS
1339                 printf("[moldyn] initializing 'static' linked cells (%d)\n",
1340                        lc->cells);
1341 #elif LOWMEM_LISTS
1342                 printf("[moldyn] initializing 'lowmem' linked cells (%d)\n",
1343                        lc->cells);
1344 #else
1345                 printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
1346                        lc->cells);
1347 #endif
1348                 printf("  x: %d x %f A\n",lc->nx,lc->x);
1349                 printf("  y: %d x %f A\n",lc->ny,lc->y);
1350                 printf("  z: %d x %f A\n",lc->nz,lc->z);
1351         }
1352
1353 #ifdef STATIC_LISTS
1354         /* list init */
1355         for(i=0;i<lc->cells;i++) {
1356                 lc->subcell[i]=malloc((MAX_ATOMS_PER_LIST+1)*sizeof(int));
1357                 if(lc->subcell[i]==NULL) {
1358                         perror("[moldyn] list init (malloc)");
1359                         return -1;
1360                 }
1361                 /*
1362                 if(i==0)
1363                         printf(" ---> %d malloc %p (%p)\n",
1364                                i,lc->subcell[0],lc->subcell);
1365                 */
1366         }
1367 #elif LOWMEM_LISTS
1368         lc->subcell->head=malloc(lc->cells*sizeof(int));
1369         if(lc->subcell->head==NULL) {
1370                 perror("[moldyn] head init (malloc)");
1371                 return -1;
1372         }
1373         lc->subcell->list=malloc(moldyn->count*sizeof(int));
1374         if(lc->subcell->list==NULL) {
1375                 perror("[moldyn] list init (malloc)");
1376                 return -1;
1377         }
1378 #else
1379         for(i=0;i<lc->cells;i++)
1380                 list_init_f(&(lc->subcell[i]));
1381 #endif
1382
1383         /* update the list */
1384         link_cell_update(moldyn);
1385
1386         return 0;
1387 }
1388
1389 int link_cell_update(t_moldyn *moldyn) {
1390
1391         int count,i,j,k;
1392         int nx,nxy;
1393         t_atom *atom;
1394         t_linkcell *lc;
1395 #ifdef STATIC_LISTS
1396         int p;
1397 #elif LOWMEM_LISTS
1398         int p;
1399 #endif
1400
1401         atom=moldyn->atom;
1402         lc=&(moldyn->lc);
1403
1404         nx=lc->nx;
1405         nxy=nx*lc->ny;
1406
1407         for(i=0;i<lc->cells;i++)
1408 #ifdef STATIC_LISTS
1409                 memset(lc->subcell[i],-1,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
1410 #elif LOWMEM_LISTS
1411                 lc->subcell->head[i]=-1;
1412 #else
1413                 list_destroy_f(&(lc->subcell[i]));
1414 #endif
1415
1416         for(count=0;count<moldyn->count;count++) {
1417                 i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
1418                 j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
1419                 k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
1420
1421 #ifdef STATIC_LISTS
1422                 p=0;
1423                 while(lc->subcell[i+j*nx+k*nxy][p]!=-1)
1424                         p++;
1425
1426                 if(p>=MAX_ATOMS_PER_LIST) {
1427                         printf("[moldyn] FATAL: amount of atoms too high!\n");
1428                         return -1;
1429                 }
1430
1431                 lc->subcell[i+j*nx+k*nxy][p]=count;
1432 #elif LOWMEM_LISTS
1433                 p=i+j*nx+k*nxy;
1434                 lc->subcell->list[count]=lc->subcell->head[p];
1435                 lc->subcell->head[p]=count;
1436 #else
1437                 list_add_immediate_f(&(lc->subcell[i+j*nx+k*nxy]),
1438                                      &(atom[count]));
1439                 /*
1440                 if(j==0&&k==0)
1441                         printf(" ---> %d %d malloc %p (%p)\n",
1442                                i,count,lc->subcell[i].current,lc->subcell);
1443                 */
1444 #endif
1445         }
1446
1447         return 0;
1448 }
1449
1450 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
1451 #ifdef STATIC_LISTS
1452                               int **cell
1453 #elif LOWMEM_LISTS
1454                               int *cell
1455 #else
1456                               t_list *cell
1457 #endif
1458                              ) {
1459
1460         t_linkcell *lc;
1461         int a;
1462         int count1,count2;
1463         int ci,cj,ck;
1464         int nx,ny,nz;
1465         int x,y,z;
1466         u8 bx,by,bz;
1467
1468         lc=&(moldyn->lc);
1469         nx=lc->nx;
1470         ny=lc->ny;
1471         nz=lc->nz;
1472         count1=1;
1473         count2=27;
1474         a=nx*ny;
1475
1476         if(i>=nx||j>=ny||k>=nz)
1477                 printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n",
1478                        i,nx,j,ny,k,nz);
1479
1480 #ifndef LOWMEM_LISTS
1481         cell[0]=lc->subcell[i+j*nx+k*a];
1482 #else
1483         cell[0]=lc->subcell->head[i+j*nx+k*a];
1484 #endif
1485         for(ci=-1;ci<=1;ci++) {
1486                 bx=0;
1487                 x=i+ci;
1488                 if((x<0)||(x>=nx)) {
1489                         x=(x+nx)%nx;
1490                         bx=1;
1491                 }
1492                 for(cj=-1;cj<=1;cj++) {
1493                         by=0;
1494                         y=j+cj;
1495                         if((y<0)||(y>=ny)) {
1496                                 y=(y+ny)%ny;
1497                                 by=1;
1498                         }
1499                         for(ck=-1;ck<=1;ck++) {
1500                                 bz=0;
1501                                 z=k+ck;
1502                                 if((z<0)||(z>=nz)) {
1503                                         z=(z+nz)%nz;
1504                                         bz=1;
1505                                 }
1506                                 if(!(ci|cj|ck)) continue;
1507                                 if(bx|by|bz) {
1508 #ifndef LOWMEM_LISTS
1509                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
1510 #else
1511                                 cell[--count2]=lc->subcell->head[x+y*nx+z*a];
1512 #endif
1513
1514                                 }
1515                                 else {
1516 #ifndef LOWMEM_LISTS
1517                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
1518 #else
1519                                 cell[count1++]=lc->subcell->head[x+y*nx+z*a];
1520 #endif
1521                                 }
1522                         }
1523                 }
1524         }
1525
1526         lc->dnlc=count1;
1527
1528         return count1;
1529 }
1530
1531 int link_cell_shutdown(t_moldyn *moldyn) {
1532
1533 #ifndef LOWMEM_LISTS
1534         int i;
1535 #endif
1536         t_linkcell *lc;
1537
1538         lc=&(moldyn->lc);
1539
1540 #if LOWMEM_LISTS
1541         free(lc->subcell->head);
1542         free(lc->subcell->list);
1543
1544 #else
1545
1546         for(i=0;i<lc->cells;i++) {
1547 #ifdef STATIC_LISTS
1548                 free(lc->subcell[i]);
1549 #else
1550                 //printf(" ---> %d free %p\n",i,lc->subcell[i].start);
1551                 list_destroy_f(&(lc->subcell[i]));
1552 #endif
1553         }
1554 #endif
1555
1556         free(lc->subcell);
1557
1558         return 0;
1559 }
1560
1561 int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
1562
1563         int count;
1564         void *ptr;
1565         t_moldyn_schedule *schedule;
1566
1567         schedule=&(moldyn->schedule);
1568         count=++(schedule->total_sched);
1569
1570         ptr=realloc(schedule->runs,count*sizeof(int));
1571         if(!ptr) {
1572                 perror("[moldyn] realloc (runs)");
1573                 return -1;
1574         }
1575         schedule->runs=ptr;
1576         schedule->runs[count-1]=runs;
1577
1578         ptr=realloc(schedule->tau,count*sizeof(double));
1579         if(!ptr) {
1580                 perror("[moldyn] realloc (tau)");
1581                 return -1;
1582         }
1583         schedule->tau=ptr;
1584         schedule->tau[count-1]=tau;
1585
1586         printf("[moldyn] schedule added:\n");
1587         printf("  number: %d | runs: %d | tau: %f\n",count-1,runs,tau);
1588
1589
1590         return 0;
1591 }
1592
1593 int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
1594
1595         moldyn->schedule.hook=hook;
1596         moldyn->schedule.hook_params=hook_params;
1597
1598         return 0;
1599 }
1600
1601 /*
1602  *
1603  * 'integration of newtons equation' - algorithms
1604  *
1605  */
1606
1607 /* start the integration */
1608
1609 int moldyn_integrate(t_moldyn *moldyn) {
1610
1611         int i;
1612         unsigned int e,m,s,v,p,t,a;
1613         t_3dvec momentum;
1614         t_moldyn_schedule *sched;
1615         t_atom *atom;
1616         int fd;
1617         char dir[128];
1618         double ds;
1619         double energy_scale;
1620         struct timeval t1,t2;
1621         //double tp;
1622
1623         sched=&(moldyn->schedule);
1624         atom=moldyn->atom;
1625
1626         /* initialize linked cell method */
1627         link_cell_init(moldyn,VERBOSE);
1628
1629         /* logging & visualization */
1630         e=moldyn->ewrite;
1631         m=moldyn->mwrite;
1632         s=moldyn->swrite;
1633         v=moldyn->vwrite;
1634         a=moldyn->awrite;
1635         p=moldyn->pwrite;
1636         t=moldyn->twrite;
1637
1638         /* sqaure of some variables */
1639         moldyn->tau_square=moldyn->tau*moldyn->tau;
1640
1641         /* get current time */
1642         gettimeofday(&t1,NULL);
1643
1644         /* calculate initial forces */
1645         potential_force_calc(moldyn);
1646 #ifdef DEBUG
1647 //return 0;
1648 #endif
1649
1650         /* some stupid checks before we actually start calculating bullshit */
1651         if(moldyn->cutoff>0.5*moldyn->dim.x)
1652                 printf("[moldyn] WARNING: cutoff > 0.5 x dim.x\n");
1653         if(moldyn->cutoff>0.5*moldyn->dim.y)
1654                 printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n");
1655         if(moldyn->cutoff>0.5*moldyn->dim.z)
1656                 printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n");
1657         if(moldyn->count) {
1658                 ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
1659                 if(ds>0.05*moldyn->nnd)
1660                 printf("[moldyn] WARNING: forces too high / tau too small!\n");
1661         }
1662
1663         /* zero absolute time */
1664         // should have right values!
1665         //moldyn->time=0.0;
1666         //moldyn->total_steps=0;
1667
1668         /* debugging, ignore */
1669         moldyn->debug=0;
1670
1671         /* tell the world */
1672         printf("[moldyn] integration start, go get a coffee ...\n");
1673
1674         /* executing the schedule */
1675         sched->count=0;
1676         while(sched->count<sched->total_sched) {
1677
1678                 /* setting amount of runs and finite time step size */
1679                 moldyn->tau=sched->tau[sched->count];
1680                 moldyn->tau_square=moldyn->tau*moldyn->tau;
1681                 moldyn->time_steps=sched->runs[sched->count];
1682
1683                 /* energy scaling factor (might change!) */
1684                 energy_scale=moldyn->count*EV;
1685
1686         /* integration according to schedule */
1687
1688         for(i=0;i<moldyn->time_steps;i++) {
1689
1690                 /* integration step */
1691                 moldyn->integrate(moldyn);
1692
1693                 /* calculate kinetic energy, temperature and pressure */
1694                 e_kin_calc(moldyn);
1695                 temperature_calc(moldyn);
1696                 virial_sum(moldyn);
1697                 pressure_calc(moldyn);
1698                 /*
1699                 thermodynamic_pressure_calc(moldyn);
1700                 printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
1701                        moldyn->tp/BAR);
1702                 */
1703
1704                 /* calculate fluctuations + averages */
1705                 average_and_fluctuation_calc(moldyn);
1706
1707                 /* p/t scaling */
1708                 if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
1709                         scale_velocity(moldyn,FALSE);
1710                 if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
1711                         scale_volume(moldyn);
1712
1713                 /* check for log & visualization */
1714                 if(e) {
1715                         if(!(moldyn->total_steps%e))
1716                                 dprintf(moldyn->efd,
1717                                         "%f %f %f %f\n",
1718                                         moldyn->time,moldyn->ekin/energy_scale,
1719                                         moldyn->energy/energy_scale,
1720                                         get_total_energy(moldyn)/energy_scale);
1721                 }
1722                 if(m) {
1723                         if(!(moldyn->total_steps%m)) {
1724                                 momentum=get_total_p(moldyn);
1725                                 dprintf(moldyn->mfd,
1726                                         "%f %f %f %f %f\n",moldyn->time,
1727                                         momentum.x,momentum.y,momentum.z,
1728                                         v3_norm(&momentum));
1729                         }
1730                 }
1731                 if(p) {
1732                         if(!(moldyn->total_steps%p)) {
1733                                 dprintf(moldyn->pfd,
1734                                         "%f %f %f %f %f %f %f\n",moldyn->time,
1735                                          moldyn->p/BAR,moldyn->p_avg/BAR,
1736                                          moldyn->gp/BAR,moldyn->gp_avg/BAR,
1737                                          moldyn->tp/BAR,moldyn->tp_avg/BAR);
1738                         }
1739                 }
1740                 if(t) {
1741                         if(!(moldyn->total_steps%t)) {
1742                                 dprintf(moldyn->tfd,
1743                                         "%f %f %f\n",
1744                                         moldyn->time,moldyn->t,moldyn->t_avg);
1745                         }
1746                 }
1747                 if(v) {
1748                         if(!(moldyn->total_steps%v)) {
1749                                 dprintf(moldyn->vfd,
1750                                         "%f %f\n",moldyn->time,moldyn->volume);
1751                         }
1752                 }
1753                 if(s) {
1754                         if(!(moldyn->total_steps%s)) {
1755                                 snprintf(dir,128,"%s/s-%07.f.save",
1756                                          moldyn->vlsdir,moldyn->time);
1757                                 fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
1758                                         S_IRUSR|S_IWUSR);
1759                                 if(fd<0) perror("[moldyn] save fd open");
1760                                 else {
1761                                         write(fd,moldyn,sizeof(t_moldyn));
1762                                         write(fd,moldyn->atom,
1763                                               moldyn->count*sizeof(t_atom));
1764                                 }
1765                                 close(fd);
1766                         }
1767                 }
1768                 if(a) {
1769                         if(!(moldyn->total_steps%a)) {
1770                                 visual_atoms(moldyn);
1771                         }
1772                 }
1773
1774                 /* display progress */
1775                 if(!(moldyn->total_steps%10)) {
1776                         /* get current time */
1777                         gettimeofday(&t2,NULL);
1778
1779 printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
1780        sched->count,i,moldyn->total_steps,
1781        moldyn->t,moldyn->t_avg,
1782        moldyn->p/BAR,moldyn->p_avg/BAR,
1783        //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
1784        moldyn->volume,
1785        (int)(t2.tv_sec-t1.tv_sec));
1786
1787                         fflush(stdout);
1788
1789                         /* copy over time */
1790                         t1=t2;
1791                 }
1792
1793                 /* increase absolute time */
1794                 moldyn->time+=moldyn->tau;
1795                 moldyn->total_steps+=1;
1796
1797         }
1798
1799                 /* check for hooks */
1800                 if(sched->hook) {
1801                         printf("\n ## schedule hook %d start ##\n",
1802                                sched->count);
1803                         sched->hook(moldyn,sched->hook_params);
1804                         printf(" ## schedule hook end ##\n");
1805                 }
1806
1807                 /* increase the schedule counter */
1808                 sched->count+=1;
1809
1810         }
1811
1812         return 0;
1813 }
1814
1815 /* velocity verlet */
1816
1817 int velocity_verlet(t_moldyn *moldyn) {
1818
1819         int i,count;
1820         double tau,tau_square,h;
1821         t_3dvec delta;
1822         t_atom *atom;
1823
1824         atom=moldyn->atom;
1825         count=moldyn->count;
1826         tau=moldyn->tau;
1827         tau_square=moldyn->tau_square;
1828
1829         for(i=0;i<count;i++) {
1830                 /* check whether fixed atom */
1831                 if(atom[i].attr&ATOM_ATTR_FP)
1832                         continue;
1833                 /* new positions */
1834                 h=0.5/atom[i].mass;
1835                 v3_scale(&delta,&(atom[i].v),tau);
1836                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1837                 v3_scale(&delta,&(atom[i].f),h*tau_square);
1838                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1839                 check_per_bound(moldyn,&(atom[i].r));
1840
1841                 /* velocities [actually v(t+tau/2)] */
1842                 v3_scale(&delta,&(atom[i].f),h*tau);
1843                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1844         }
1845
1846         /* criticial check */
1847         moldyn_bc_check(moldyn);
1848
1849         /* neighbour list update */
1850         link_cell_update(moldyn);
1851
1852         /* forces depending on chosen potential */
1853 #ifndef ALBE_FAST
1854         potential_force_calc(moldyn);
1855 #else
1856         albe_potential_force_calc(moldyn);
1857 #endif
1858
1859         for(i=0;i<count;i++) {
1860                 /* check whether fixed atom */
1861                 if(atom[i].attr&ATOM_ATTR_FP)
1862                         continue;
1863                 /* again velocities [actually v(t+tau)] */
1864                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
1865                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1866         }
1867
1868         return 0;
1869 }
1870
1871
1872 /*
1873  *
1874  * potentials & corresponding forces & virial routine
1875  *
1876  */
1877
1878 /* generic potential and force calculation */
1879
1880 int potential_force_calc(t_moldyn *moldyn) {
1881
1882         int i,j,k,count;
1883         t_atom *itom,*jtom,*ktom;
1884         t_virial *virial;
1885         t_linkcell *lc;
1886 #ifdef STATIC_LISTS
1887         int *neighbour_i[27];
1888         int p,q;
1889         t_atom *atom;
1890 #elif LOWMEM_LISTS
1891         int neighbour_i[27];
1892         int p,q;
1893 #else
1894         t_list neighbour_i[27];
1895         t_list neighbour_i2[27];
1896         t_list *this,*that;
1897 #endif
1898         u8 bc_ij,bc_ik;
1899         int dnlc;
1900
1901         count=moldyn->count;
1902         itom=moldyn->atom;
1903         lc=&(moldyn->lc);
1904 #ifdef STATIC_LISTS
1905         atom=moldyn->atom;
1906 #endif
1907
1908         /* reset energy */
1909         moldyn->energy=0.0;
1910
1911         /* reset global virial */
1912         memset(&(moldyn->gvir),0,sizeof(t_virial));
1913
1914         /* reset force, site energy and virial of every atom */
1915 #ifdef PARALLEL
1916         i=omp_get_thread_num();
1917         #pragma omp parallel for private(virial)
1918 #endif
1919         for(i=0;i<count;i++) {
1920
1921                 /* reset force */
1922                 v3_zero(&(itom[i].f));
1923
1924                 /* reset virial */
1925                 virial=(&(itom[i].virial));
1926                 virial->xx=0.0;
1927                 virial->yy=0.0;
1928                 virial->zz=0.0;
1929                 virial->xy=0.0;
1930                 virial->xz=0.0;
1931                 virial->yz=0.0;
1932
1933                 /* reset site energy */
1934                 itom[i].e=0.0;
1935
1936         }
1937
1938         /* get energy, force and virial of every atom */
1939
1940         /* first (and only) loop over atoms i */
1941         for(i=0;i<count;i++) {
1942
1943                 /* single particle potential/force */
1944                 if(itom[i].attr&ATOM_ATTR_1BP)
1945                         if(moldyn->func1b)
1946                                 moldyn->func1b(moldyn,&(itom[i]));
1947
1948                 if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
1949                         continue;
1950
1951                 /* 2 body pair potential/force */
1952
1953                 link_cell_neighbour_index(moldyn,
1954                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
1955                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
1956                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
1957                                           neighbour_i);
1958
1959                 dnlc=lc->dnlc;
1960
1961                 /* first loop over atoms j */
1962                 if(moldyn->func2b) {
1963                         for(j=0;j<27;j++) {
1964
1965                                 bc_ij=(j<dnlc)?0:1;
1966 #ifdef STATIC_LISTS
1967                                 p=0;
1968
1969                                 while(neighbour_i[j][p]!=-1) {
1970
1971                                         jtom=&(atom[neighbour_i[j][p]]);
1972                                         p++;
1973 #elif LOWMEM_LISTS
1974                                 p=neighbour_i[j];
1975
1976                                 while(p!=-1) {
1977
1978                                         jtom=&(itom[p]);
1979                                         p=lc->subcell->list[p];
1980 #else
1981                                 this=&(neighbour_i[j]);
1982                                 list_reset_f(this);
1983
1984                                 if(this->start==NULL)
1985                                         continue;
1986
1987                                 do {
1988                                         jtom=this->current->data;
1989 #endif
1990
1991                                         if(jtom==&(itom[i]))
1992                                                 continue;
1993
1994                                         if((jtom->attr&ATOM_ATTR_2BP)&
1995                                            (itom[i].attr&ATOM_ATTR_2BP)) {
1996                                                 moldyn->func2b(moldyn,
1997                                                                &(itom[i]),
1998                                                                jtom,
1999                                                                bc_ij);
2000                                         }
2001 #ifdef STATIC_LISTS
2002                                 }
2003 #elif LOWMEM_LISTS
2004                                 }
2005 #else
2006                                 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2007 #endif
2008
2009                         }
2010                 }
2011
2012                 /* 3 body potential/force */
2013
2014                 if(!(itom[i].attr&ATOM_ATTR_3BP))
2015                         continue;
2016
2017                 /* copy the neighbour lists */
2018 #ifdef STATIC_LISTS
2019                 /* no copy needed for static lists */
2020 #elif LOWMEM_LISTS
2021                 /* no copy needed for lowmem lists */
2022 #else
2023                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
2024 #endif
2025
2026                 /* second loop over atoms j */
2027                 for(j=0;j<27;j++) {
2028
2029                         bc_ij=(j<dnlc)?0:1;
2030 #ifdef STATIC_LISTS
2031                         p=0;
2032
2033                         while(neighbour_i[j][p]!=-1) {
2034
2035                                 jtom=&(atom[neighbour_i[j][p]]);
2036                                 p++;
2037 #elif LOWMEM_LISTS
2038                                 p=neighbour_i[j];
2039
2040                                 while(p!=-1) {
2041
2042                                         jtom=&(itom[p]);
2043                                         p=lc->subcell->list[p];
2044 #else
2045                         this=&(neighbour_i[j]);
2046                         list_reset_f(this);
2047
2048                         if(this->start==NULL)
2049                                 continue;
2050
2051                         do {
2052
2053                                 jtom=this->current->data;
2054 #endif
2055
2056                                 if(jtom==&(itom[i]))
2057                                         continue;
2058
2059                                 if(!(jtom->attr&ATOM_ATTR_3BP))
2060                                         continue;
2061
2062                                 /* reset 3bp run */
2063                                 moldyn->run3bp=1;
2064
2065                                 if(moldyn->func3b_j1)
2066                                         moldyn->func3b_j1(moldyn,
2067                                                           &(itom[i]),
2068                                                           jtom,
2069                                                           bc_ij);
2070
2071                                 /* in first j loop, 3bp run can be skipped */
2072                                 if(!(moldyn->run3bp))
2073                                         continue;
2074
2075                                 /* first loop over atoms k */
2076                                 if(moldyn->func3b_k1) {
2077
2078                                 for(k=0;k<27;k++) {
2079
2080                                         bc_ik=(k<dnlc)?0:1;
2081 #ifdef STATIC_LISTS
2082                                         q=0;
2083
2084                                         while(neighbour_i[k][q]!=-1) {
2085
2086                                                 ktom=&(atom[neighbour_i[k][q]]);
2087                                                 q++;
2088 #elif LOWMEM_LISTS
2089                                         q=neighbour_i[k];
2090
2091                                         while(q!=-1) {
2092
2093                                                 ktom=&(itom[q]);
2094                                                 q=lc->subcell->list[q];
2095 #else
2096                                         that=&(neighbour_i2[k]);
2097                                         list_reset_f(that);
2098
2099                                         if(that->start==NULL)
2100                                                 continue;
2101
2102                                         do {
2103                                                 ktom=that->current->data;
2104 #endif
2105
2106                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
2107                                                         continue;
2108
2109                                                 if(ktom==jtom)
2110                                                         continue;
2111
2112                                                 if(ktom==&(itom[i]))
2113                                                         continue;
2114
2115                                                 moldyn->func3b_k1(moldyn,
2116                                                                   &(itom[i]),
2117                                                                   jtom,
2118                                                                   ktom,
2119                                                                   bc_ik|bc_ij);
2120 #ifdef STATIC_LISTS
2121                                         }
2122 #elif LOWMEM_LISTS
2123                                         }
2124 #else
2125                                         } while(list_next_f(that)!=\
2126                                                 L_NO_NEXT_ELEMENT);
2127 #endif
2128
2129                                 }
2130
2131                                 }
2132
2133                                 if(moldyn->func3b_j2)
2134                                         moldyn->func3b_j2(moldyn,
2135                                                           &(itom[i]),
2136                                                           jtom,
2137                                                           bc_ij);
2138
2139                                 /* second loop over atoms k */
2140                                 if(moldyn->func3b_k2) {
2141
2142                                 for(k=0;k<27;k++) {
2143
2144                                         bc_ik=(k<dnlc)?0:1;
2145 #ifdef STATIC_LISTS
2146                                         q=0;
2147
2148                                         while(neighbour_i[k][q]!=-1) {
2149
2150                                                 ktom=&(atom[neighbour_i[k][q]]);
2151                                                 q++;
2152 #elif LOWMEM_LISTS
2153                                         q=neighbour_i[k];
2154
2155                                         while(q!=-1) {
2156
2157                                                 ktom=&(itom[q]);
2158                                                 q=lc->subcell->list[q];
2159 #else
2160                                         that=&(neighbour_i2[k]);
2161                                         list_reset_f(that);
2162
2163                                         if(that->start==NULL)
2164                                                 continue;
2165
2166                                         do {
2167                                                 ktom=that->current->data;
2168 #endif
2169
2170                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
2171                                                         continue;
2172
2173                                                 if(ktom==jtom)
2174                                                         continue;
2175
2176                                                 if(ktom==&(itom[i]))
2177                                                         continue;
2178
2179                                                 moldyn->func3b_k2(moldyn,
2180                                                                   &(itom[i]),
2181                                                                   jtom,
2182                                                                   ktom,
2183                                                                   bc_ik|bc_ij);
2184
2185 #ifdef STATIC_LISTS
2186                                         }
2187 #elif LOWMEM_LISTS
2188                                         }
2189 #else
2190                                         } while(list_next_f(that)!=\
2191                                                 L_NO_NEXT_ELEMENT);
2192 #endif
2193
2194                                 }
2195
2196                                 }
2197
2198                                 /* 2bp post function */
2199                                 if(moldyn->func3b_j3) {
2200                                         moldyn->func3b_j3(moldyn,
2201                                                           &(itom[i]),
2202                                                           jtom,bc_ij);
2203                                 }
2204 #ifdef STATIC_LISTS
2205                         }
2206 #elif LOWMEM_LISTS
2207                         }
2208 #else
2209                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2210 #endif
2211
2212                 }
2213
2214 #ifdef DEBUG
2215         //printf("\n\n");
2216 #endif
2217 #ifdef VDEBUG
2218         printf("\n\n");
2219 #endif
2220
2221         }
2222
2223 #ifdef DEBUG
2224         //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
2225         if(moldyn->time>DSTART&&moldyn->time<DEND) {
2226                 printf("force:\n");
2227                 printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
2228                 printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
2229                 printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
2230         }
2231 #endif
2232
2233         /* some postprocessing */
2234 #ifdef PARALLEL
2235         #pragma omp parallel for
2236 #endif
2237         for(i=0;i<count;i++) {
2238                 /* calculate global virial */
2239                 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
2240                 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
2241                 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
2242                 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
2243                 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
2244                 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
2245
2246                 /* check forces regarding the given timestep */
2247                 if(v3_norm(&(itom[i].f))>\
2248                    0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
2249                         printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
2250                                i);
2251         }
2252
2253         return 0;
2254 }
2255
2256 /*
2257  * virial calculation
2258  */
2259
2260 //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
2261 int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
2262
2263         a->virial.xx+=f->x*d->x;
2264         a->virial.yy+=f->y*d->y;
2265         a->virial.zz+=f->z*d->z;
2266         a->virial.xy+=f->x*d->y;
2267         a->virial.xz+=f->x*d->z;
2268         a->virial.yz+=f->y*d->z;
2269
2270         return 0;
2271 }
2272
2273 /*
2274  * periodic boundary checking
2275  */
2276
2277 //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
2278 int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
2279
2280         double x,y,z;
2281         t_3dvec *dim;
2282
2283         dim=&(moldyn->dim);
2284
2285         x=dim->x/2;
2286         y=dim->y/2;
2287         z=dim->z/2;
2288
2289         if(moldyn->status&MOLDYN_STAT_PBX) {
2290                 if(a->x>=x) a->x-=dim->x;
2291                 else if(-a->x>x) a->x+=dim->x;
2292         }
2293         if(moldyn->status&MOLDYN_STAT_PBY) {
2294                 if(a->y>=y) a->y-=dim->y;
2295                 else if(-a->y>y) a->y+=dim->y;
2296         }
2297         if(moldyn->status&MOLDYN_STAT_PBZ) {
2298                 if(a->z>=z) a->z-=dim->z;
2299                 else if(-a->z>z) a->z+=dim->z;
2300         }
2301
2302         return 0;
2303 }
2304
2305 /*
2306  * debugging / critical check functions
2307  */
2308
2309 int moldyn_bc_check(t_moldyn *moldyn) {
2310
2311         t_atom *atom;
2312         t_3dvec *dim;
2313         int i;
2314         double x;
2315         u8 byte;
2316         int j,k;
2317
2318         atom=moldyn->atom;
2319         dim=&(moldyn->dim);
2320         x=dim->x/2;
2321
2322         for(i=0;i<moldyn->count;i++) {
2323                 if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {
2324                         printf("FATAL: atom %d: x: %.20f (%.20f)\n",
2325                                i,atom[i].r.x,dim->x/2);
2326                         printf("diagnostic:\n");
2327                         printf("-----------\natom.r.x:\n");
2328                         for(j=0;j<8;j++) {
2329                                 memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
2330                                 for(k=0;k<8;k++)
2331                                         printf("%d%c",
2332                                         ((byte)&(1<<k))?1:0,
2333                                         (k==7)?'\n':'|');
2334                         }
2335                         printf("---------------\nx=dim.x/2:\n");
2336                         for(j=0;j<8;j++) {
2337                                 memcpy(&byte,(u8 *)(&x)+j,1);
2338                                 for(k=0;k<8;k++)
2339                                         printf("%d%c",
2340                                         ((byte)&(1<<k))?1:0,
2341                                         (k==7)?'\n':'|');
2342                         }
2343                         if(atom[i].r.x==x) printf("the same!\n");
2344                         else printf("different!\n");
2345                 }
2346                 if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
2347                         printf("FATAL: atom %d: y: %.20f (%.20f)\n",
2348                                i,atom[i].r.y,dim->y/2);
2349                 if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
2350                         printf("FATAL: atom %d: z: %.20f (%.20f)\n",
2351                                i,atom[i].r.z,dim->z/2);
2352         }
2353
2354         return 0;
2355 }
2356
2357 /*
2358  * restore function
2359  */
2360
2361 int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
2362
2363         int fd;
2364         int cnt,size;
2365         int fsize;
2366         int corr;
2367
2368         fd=open(file,O_RDONLY);
2369         if(fd<0) {
2370                 perror("[moldyn] load save file open");
2371                 return fd;
2372         }
2373
2374         fsize=lseek(fd,0,SEEK_END);
2375         lseek(fd,0,SEEK_SET);
2376
2377         size=sizeof(t_moldyn);
2378
2379         while(size) {
2380                 cnt=read(fd,moldyn,size);
2381                 if(cnt<0) {
2382                         perror("[moldyn] load save file read (moldyn)");
2383                         return cnt;
2384                 }
2385                 size-=cnt;
2386         }
2387
2388         size=moldyn->count*sizeof(t_atom);
2389
2390         /* correcting possible atom data offset */
2391         corr=0;
2392         if(fsize!=sizeof(t_moldyn)+size) {
2393                 corr=fsize-sizeof(t_moldyn)-size;
2394                 printf("[moldyn] WARNING: lsf (illegal file size)\n");
2395                 printf("  moifying offset:\n");
2396                 printf("  - current pos: %d\n",sizeof(t_moldyn));
2397                 printf("  - atom size: %d\n",size);
2398                 printf("  - file size: %d\n",fsize);
2399                 printf("  => correction: %d\n",corr);
2400                 lseek(fd,corr,SEEK_CUR);
2401         }
2402
2403         moldyn->atom=(t_atom *)malloc(size);
2404         if(moldyn->atom==NULL) {
2405                 perror("[moldyn] load save file malloc (atoms)");
2406                 return -1;
2407         }
2408
2409         while(size) {
2410                 cnt=read(fd,moldyn->atom,size);
2411                 if(cnt<0) {
2412                         perror("[moldyn] load save file read (atoms)");
2413                         return cnt;
2414                 }
2415                 size-=cnt;
2416         }
2417
2418         // hooks etc ...
2419
2420         return 0;
2421 }
2422
2423 int moldyn_free_save_file(t_moldyn *moldyn) {
2424
2425         free(moldyn->atom);
2426
2427         return 0;
2428 }
2429
2430 int moldyn_load(t_moldyn *moldyn) {
2431
2432         // later ...
2433
2434         return 0;
2435 }
2436
2437 /*
2438  * function to find/callback all combinations of 2 body bonds
2439  */
2440
2441 int process_2b_bonds(t_moldyn *moldyn,void *data,
2442                      int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
2443                                     void *data,u8 bc)) {
2444
2445         t_linkcell *lc;
2446 #ifdef STATIC_LISTS
2447         int *neighbour[27];
2448         int p;
2449 #elif LOWMEM_LISTS
2450         int neighbour[27];
2451         int p;
2452 #else
2453         t_list neighbour[27];
2454         t_list *this;
2455 #endif
2456         u8 bc;
2457         t_atom *itom,*jtom;
2458         int i,j;
2459
2460         lc=&(moldyn->lc);
2461         itom=moldyn->atom;
2462
2463         for(i=0;i<moldyn->count;i++) {
2464                 /* neighbour indexing */
2465                 link_cell_neighbour_index(moldyn,
2466                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
2467                                           (itom[i].r.y+moldyn->dim.y/2)/lc->x,
2468                                           (itom[i].r.z+moldyn->dim.z/2)/lc->x,
2469                                           neighbour);
2470
2471                 for(j=0;j<27;j++) {
2472
2473                         bc=(j<lc->dnlc)?0:1;
2474
2475 #ifdef STATIC_LISTS
2476                         p=0;
2477
2478                         while(neighbour[j][p]!=-1) {
2479
2480                                 jtom=&(moldyn->atom[neighbour[j][p]]);
2481                                 p++;
2482 #elif LOWMEM_LISTS
2483                         p=neighbour[j];
2484
2485                         while(p!=-1) {
2486
2487                                 jtom=&(itom[p]);
2488                                 p=lc->subcell->list[p];
2489 #else
2490                         this=&(neighbour[j]);
2491                         list_reset_f(this);
2492
2493                         if(this->start==NULL)
2494                                 continue;
2495
2496                         do {
2497
2498                                 jtom=this->current->data;
2499 #endif
2500
2501                                 /* process bond */
2502                                 process(moldyn,&(itom[i]),jtom,data,bc);
2503
2504 #ifdef STATIC_LISTS
2505                         }
2506 #elif LOWMEM_LISTS
2507                         }
2508 #else
2509                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2510 #endif
2511                 }
2512         }
2513
2514         return 0;
2515
2516 }
2517
2518 /*
2519  * function to find neighboured atoms
2520  */
2521
2522 int process_neighbours(t_moldyn *moldyn,void *data,t_atom *atom,
2523                        int (*process)(t_moldyn *moldyn,t_atom *atom,t_atom *natom,
2524                                       void *data,u8 bc)) {
2525
2526         t_linkcell *lc;
2527 #ifdef STATIC_LISTS
2528         int *neighbour[27];
2529         int p;
2530 #elif LOWMEM_LISTS
2531         int neighbour[27];
2532         int p;
2533 #else
2534         t_list neighbour[27];
2535         t_list *this;
2536 #endif
2537         u8 bc;
2538         t_atom *natom;
2539         int j;
2540
2541         lc=&(moldyn->lc);
2542
2543         /* neighbour indexing */
2544         link_cell_neighbour_index(moldyn,
2545                                   (atom->r.x+moldyn->dim.x/2)/lc->x,
2546                                   (atom->r.y+moldyn->dim.y/2)/lc->x,
2547                                   (atom->r.z+moldyn->dim.z/2)/lc->x,
2548                                   neighbour);
2549
2550         for(j=0;j<27;j++) {
2551
2552                 bc=(j<lc->dnlc)?0:1;
2553
2554 #ifdef STATIC_LISTS
2555                 p=0;
2556
2557                 while(neighbour[j][p]!=-1) {
2558
2559                         natom=&(moldyn->atom[neighbour[j][p]]);
2560                         p++;
2561 #elif LOWMEM_LISTS
2562                 p=neighbour[j];
2563
2564                 while(p!=-1) {
2565
2566                         natom=&(moldyn->atom[p]);
2567                         p=lc->subcell->list[p];
2568 #else
2569                 this=&(neighbour[j]);
2570                 list_reset_f(this);
2571
2572                 if(this->start==NULL)
2573                         continue;
2574
2575                 do {
2576
2577                         natom=this->current->data;
2578 #endif
2579
2580                         /* process bond */
2581                         process(moldyn,atom,natom,data,bc);
2582
2583 #ifdef STATIC_LISTS
2584                 }
2585 #elif LOWMEM_LISTS
2586                 }
2587 #else
2588                 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2589 #endif
2590         }
2591
2592         return 0;
2593
2594 }
2595
2596 /*
2597  * post processing functions
2598  */
2599
2600 int get_line(int fd,char *line,int max) {
2601
2602         int count,ret;
2603
2604         count=0;
2605
2606         while(1) {
2607                 if(count==max) return count;
2608                 ret=read(fd,line+count,1);
2609                 if(ret<=0) return ret;
2610                 if(line[count]=='\n') {
2611                         memset(line+count,0,max-count-1);
2612                         //line[count]='\0';
2613                         return count+1;
2614                 }
2615                 count+=1;
2616         }
2617 }
2618
2619 int pair_correlation_init(t_moldyn *moldyn,double dr) {
2620
2621
2622         return 0;
2623 }
2624
2625 int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
2626
2627         int i;
2628         t_atom *atom;
2629         t_3dvec dist;
2630         double d2;
2631         int a_cnt;
2632         int b_cnt;
2633
2634         atom=moldyn->atom;
2635         dc[0]=0;
2636         dc[1]=0;
2637         dc[2]=0;
2638         a_cnt=0;
2639         b_cnt=0;
2640
2641         for(i=0;i<moldyn->count;i++) {
2642
2643                 v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
2644                 check_per_bound(moldyn,&dist);
2645                 d2=v3_absolute_square(&dist);
2646
2647                 if(atom[i].brand) {
2648                         b_cnt+=1;
2649                         dc[1]+=d2;
2650                 }
2651                 else {
2652                         a_cnt+=1;
2653                         dc[0]+=d2;
2654                 }
2655
2656                 dc[2]+=d2;
2657         }
2658
2659         dc[0]*=(1.0/(6.0*moldyn->time*a_cnt));
2660         dc[1]*=(1.0/(6.0*moldyn->time*b_cnt));
2661         dc[2]*=(1.0/(6.0*moldyn->time*moldyn->count));
2662
2663         return 0;
2664 }
2665
2666 int bonding_analyze(t_moldyn *moldyn,double *cnt) {
2667
2668         return 0;
2669 }
2670
2671 int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
2672                                        t_atom *jtom,void *data,u8 bc) {
2673
2674         t_3dvec dist;
2675         double d;
2676         int s;
2677         t_pcc *pcc;
2678
2679         /* only count pairs once,
2680          * skip same atoms */
2681         if(itom->tag>=jtom->tag)
2682                 return 0;
2683
2684         /*
2685          * pair correlation calc
2686          */
2687
2688         /* get pcc data */
2689         pcc=data;
2690
2691         /* distance */
2692         v3_sub(&dist,&(jtom->r),&(itom->r));
2693         if(bc) check_per_bound(moldyn,&dist);
2694         d=v3_absolute_square(&dist);
2695
2696         /* ignore if greater cutoff */
2697         if(d>moldyn->cutoff_square)
2698                 return 0;
2699
2700         /* fill the slots */
2701         d=sqrt(d);
2702         s=(int)(d/pcc->dr);
2703
2704         /* should never happen but it does 8) -
2705          * related to -ffloat-store problem! */
2706         if(s>=pcc->o1) {
2707                 printf("[moldyn] WARNING: pcc (%d/%d)",
2708                        s,pcc->o1);
2709                 printf("\n");
2710                 s=pcc->o1-1;
2711         }
2712
2713         if(itom->brand!=jtom->brand) {
2714                 /* mixed */
2715                 pcc->stat[s]+=1;
2716         }
2717         else {
2718                 /* type a - type a bonds */
2719                 if(itom->brand==0)
2720                         pcc->stat[s+pcc->o1]+=1;
2721                 else
2722                 /* type b - type b bonds */
2723                         pcc->stat[s+pcc->o2]+=1;
2724         }
2725
2726         return 0;
2727 }
2728
2729 int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
2730
2731         t_pcc pcc;
2732         double norm;
2733         int i;
2734
2735         pcc.dr=dr;
2736         pcc.o1=moldyn->cutoff/dr;
2737         pcc.o2=2*pcc.o1;
2738
2739         if(pcc.o1*dr<=moldyn->cutoff)
2740                 printf("[moldyn] WARNING: pcc (low #slots)\n");
2741
2742         printf("[moldyn] pair correlation calc info:\n");
2743         printf("  time: %f\n",moldyn->time);
2744         printf("  count: %d\n",moldyn->count);
2745         printf("  cutoff: %f\n",moldyn->cutoff);
2746         printf("  temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg);
2747
2748         if(ptr!=NULL) {
2749                 pcc.stat=(double *)ptr;
2750         }
2751         else {
2752                 pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double));
2753                 if(pcc.stat==NULL) {
2754                         perror("[moldyn] pair correlation malloc");
2755                         return -1;
2756                 }
2757         }
2758
2759         memset(pcc.stat,0,3*pcc.o1*sizeof(double));
2760
2761         /* process */
2762         process_2b_bonds(moldyn,&pcc,calculate_pair_correlation_process);
2763
2764         /* normalization */
2765         for(i=1;i<pcc.o1;i++) {
2766                  // normalization: 4 pi r^2 dr
2767                  // here: not double counting pairs -> 2 pi r r dr
2768                  // ... and actually it's a constant times r^2
2769                 norm=i*i*dr*dr;
2770                 pcc.stat[i]/=norm;
2771                 pcc.stat[pcc.o1+i]/=norm;
2772                 pcc.stat[pcc.o2+i]/=norm;
2773         }
2774         /* */
2775
2776         if(ptr==NULL) {
2777                 /* todo: store/print pair correlation function */
2778                 free(pcc.stat);
2779         }
2780
2781         return 0;
2782 }
2783
2784 int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
2785                          void *data,u8 bc) {
2786
2787         t_ba *ba;
2788         t_3dvec dist;
2789         double d;
2790
2791         if(itom->tag>=jtom->tag)
2792                 return 0;
2793
2794         /* distance */
2795         v3_sub(&dist,&(jtom->r),&(itom->r));
2796         if(bc) check_per_bound(moldyn,&dist);
2797         d=v3_absolute_square(&dist);
2798
2799         /* ignore if greater or equal cutoff */
2800         if(d>moldyn->cutoff_square)
2801                 return 0;
2802
2803         /* check for potential bond */
2804         if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
2805                 return 0;
2806
2807         /* now count this bonding ... */
2808         ba=data;
2809
2810         /* increase total bond counter
2811          * ... double counting!
2812          */
2813         ba->tcnt+=2;
2814
2815         if(itom->brand==0)
2816                 ba->acnt[jtom->tag]+=1;
2817         else
2818                 ba->bcnt[jtom->tag]+=1;
2819
2820         if(jtom->brand==0)
2821                 ba->acnt[itom->tag]+=1;
2822         else
2823                 ba->bcnt[itom->tag]+=1;
2824
2825         return 0;
2826 }
2827
2828 int bond_analyze(t_moldyn *moldyn,double *quality) {
2829
2830         // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
2831
2832         int qcnt;
2833         int ccnt,cset;
2834         t_ba ba;
2835         int i;
2836         t_atom *atom;
2837
2838         ba.acnt=malloc(moldyn->count*sizeof(int));
2839         if(ba.acnt==NULL) {
2840                 perror("[moldyn] bond analyze malloc (a)");
2841                 return -1;
2842         }
2843         memset(ba.acnt,0,moldyn->count*sizeof(int));
2844
2845         ba.bcnt=malloc(moldyn->count*sizeof(int));
2846         if(ba.bcnt==NULL) {
2847                 perror("[moldyn] bond analyze malloc (b)");
2848                 return -1;
2849         }
2850         memset(ba.bcnt,0,moldyn->count*sizeof(int));
2851
2852         ba.tcnt=0;
2853         qcnt=0;
2854         ccnt=0;
2855         cset=0;
2856
2857         atom=moldyn->atom;
2858
2859         process_2b_bonds(moldyn,&ba,bond_analyze_process);
2860
2861         for(i=0;i<moldyn->count;i++) {
2862                 if(atom[i].brand==0) {
2863                         if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
2864                                 qcnt+=4;
2865                 }
2866                 else {
2867                         if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
2868                                 qcnt+=4;
2869                                 ccnt+=1;
2870                         }
2871                         cset+=1;
2872                 }
2873         }
2874
2875         printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
2876         printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
2877
2878         if(quality) {
2879                 quality[0]=1.0*ccnt/cset;
2880                 quality[1]=1.0*qcnt/ba.tcnt;
2881         }
2882         else {
2883                 printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
2884                 printf("[moldyn] bond analyze:   tot_q=%f\n",1.0*qcnt/ba.tcnt);
2885         }
2886
2887         return 0;
2888 }
2889
2890 /*
2891  * visualization code
2892  */
2893
2894 int visual_init(t_moldyn *moldyn,char *filebase) {
2895
2896         strncpy(moldyn->vis.fb,filebase,128);
2897
2898         return 0;
2899 }
2900
2901 int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
2902                          void *data,u8 bc) {
2903
2904         t_vb *vb;
2905
2906         vb=data;
2907
2908         if(itom->tag>=jtom->tag)
2909                 return 0;
2910
2911         if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
2912                 return 0;
2913
2914         if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB))
2915                 dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n",
2916                         itom->r.x,itom->r.y,itom->r.z,
2917                         jtom->r.x,jtom->r.y,jtom->r.z);
2918
2919         return 0;
2920 }
2921
2922 int visual_atoms(t_moldyn *moldyn) {
2923
2924         int i;
2925         char file[128+64];
2926         t_3dvec dim;
2927         double help;
2928         t_visual *v;
2929         t_atom *atom;
2930         t_vb vb;
2931
2932         v=&(moldyn->vis);
2933         dim.x=v->dim.x;
2934         dim.y=v->dim.y;
2935         dim.z=v->dim.z;
2936         atom=moldyn->atom;
2937
2938         help=(dim.x+dim.y);
2939
2940         sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
2941         vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
2942         if(vb.fd<0) {
2943                 perror("open visual save file fd");
2944                 return -1;
2945         }
2946
2947         /* write the actual data file */
2948
2949         // povray header
2950         dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n",
2951                 moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
2952
2953         // atomic configuration
2954         for(i=0;i<moldyn->count;i++)
2955                 // atom type, positions, color and kinetic energy
2956                 dprintf(vb.fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
2957                                                     atom[i].r.x,
2958                                                     atom[i].r.y,
2959                                                     atom[i].r.z,
2960                                                     pse_col[atom[i].element],
2961                                                     atom[i].ekin);
2962
2963         // bonds between atoms
2964         process_2b_bonds(moldyn,&vb,visual_bonds_process);
2965
2966         // boundaries
2967         if(dim.x) {
2968                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2969                         -dim.x/2,-dim.y/2,-dim.z/2,
2970                         dim.x/2,-dim.y/2,-dim.z/2);
2971                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2972                         -dim.x/2,-dim.y/2,-dim.z/2,
2973                         -dim.x/2,dim.y/2,-dim.z/2);
2974                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2975                         dim.x/2,dim.y/2,-dim.z/2,
2976                         dim.x/2,-dim.y/2,-dim.z/2);
2977                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2978                         -dim.x/2,dim.y/2,-dim.z/2,
2979                         dim.x/2,dim.y/2,-dim.z/2);
2980
2981                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2982                         -dim.x/2,-dim.y/2,dim.z/2,
2983                         dim.x/2,-dim.y/2,dim.z/2);
2984                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2985                         -dim.x/2,-dim.y/2,dim.z/2,
2986                         -dim.x/2,dim.y/2,dim.z/2);
2987                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2988                         dim.x/2,dim.y/2,dim.z/2,
2989                         dim.x/2,-dim.y/2,dim.z/2);
2990                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2991                         -dim.x/2,dim.y/2,dim.z/2,
2992                         dim.x/2,dim.y/2,dim.z/2);
2993
2994                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2995                         -dim.x/2,-dim.y/2,dim.z/2,
2996                         -dim.x/2,-dim.y/2,-dim.z/2);
2997                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
2998                         -dim.x/2,dim.y/2,dim.z/2,
2999                         -dim.x/2,dim.y/2,-dim.z/2);
3000                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3001                         dim.x/2,-dim.y/2,dim.z/2,
3002                         dim.x/2,-dim.y/2,-dim.z/2);
3003                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3004                         dim.x/2,dim.y/2,dim.z/2,
3005                         dim.x/2,dim.y/2,-dim.z/2);
3006         }
3007
3008         close(vb.fd);
3009
3010         return 0;
3011 }
3012
3013 /*
3014  * fpu cntrol functions
3015  */
3016
3017 // set rounding to double (eliminates -ffloat-store!)
3018 int fpu_set_rtd(void) {
3019
3020         fpu_control_t ctrl;
3021
3022         _FPU_GETCW(ctrl);
3023
3024         ctrl&=~_FPU_EXTENDED;
3025         ctrl|=_FPU_DOUBLE;
3026
3027         _FPU_SETCW(ctrl);
3028
3029         return 0;
3030 }
3031