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