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