moldyn.c

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