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