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                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
2287                 v3_scale(&delta,&(atom[i].f),h*tau_square);
2288 #ifdef CONSTRAINT_110_5832
2289                 if(i==5832) {
2290                         delta.y=-delta.x;
2291                 }
2292 #endif
2293                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
2294                 //check_per_bound_and_care_for_pbc(moldyn,&(atom[i]));
2295                 check_per_bound(moldyn,&(atom[i].r));
2296
2297                 /* velocities [actually v(t+tau/2)] */
2298                 v3_scale(&delta,&(atom[i].f),h*tau);
2299                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
2300         }
2301
2302         /* criticial check */
2303         moldyn_bc_check(moldyn);
2304
2305         /* neighbour list update */
2306         link_cell_update(moldyn);
2307
2308         /* forces depending on chosen potential */
2309 #ifndef ALBE_FAST
2310         // if albe, use albe force calc routine
2311         //if(moldyn->func3b_j1==albe_mult_3bp_j1)
2312         //      albe_potential_force_calc(moldyn);
2313         //else
2314                 potential_force_calc(moldyn);
2315 #else
2316         albe_potential_force_calc(moldyn);
2317 #endif
2318
2319 #ifdef CONSTRAINT_110_5832
2320         if(count==5833) {
2321                 atom[5832].f.x=0.5*(atom[5832].f.x-atom[5832].f.y);
2322                 atom[5832].f.y=-atom[5832].f.x;
2323         }
2324 #endif
2325         for(i=0;i<count;i++) {
2326                 /* check whether fixed atom */
2327                 if(atom[i].attr&ATOM_ATTR_FP)
2328                         continue;
2329                 /* again velocities [actually v(t+tau)] */
2330                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
2331                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
2332         }
2333
2334         return 0;
2335 }
2336
2337
2338 /*
2339  *
2340  * potentials & corresponding forces & virial routine
2341  *
2342  */
2343
2344 /* generic potential and force calculation */
2345
2346 int potential_force_calc(t_moldyn *moldyn) {
2347
2348         int i,j,k,count;
2349         t_atom *itom,*jtom,*ktom;
2350         t_virial *virial;
2351         t_linkcell *lc;
2352 #ifdef STATIC_LISTS
2353         int *neighbour_i[27];
2354         int p,q;
2355         t_atom *atom;
2356 #elif LOWMEM_LISTS
2357         int neighbour_i[27];
2358         int p,q;
2359 #else
2360         t_list neighbour_i[27];
2361         t_list neighbour_i2[27];
2362         t_list *this,*that;
2363 #endif
2364         u8 bc_ij,bc_ik;
2365         int dnlc;
2366
2367         count=moldyn->count;
2368         itom=moldyn->atom;
2369         lc=&(moldyn->lc);
2370 #ifdef STATIC_LISTS
2371         atom=moldyn->atom;
2372 #endif
2373
2374         /* reset energy */
2375         moldyn->energy=0.0;
2376
2377         /* reset global virial */
2378         memset(&(moldyn->gvir),0,sizeof(t_virial));
2379
2380         /* reset force, site energy and virial of every atom */
2381 #ifdef PARALLEL
2382         i=omp_get_thread_num();
2383         #pragma omp parallel for private(virial)
2384 #endif
2385         for(i=0;i<count;i++) {
2386
2387                 /* reset force */
2388                 v3_zero(&(itom[i].f));
2389
2390                 /* reset virial */
2391                 virial=(&(itom[i].virial));
2392                 virial->xx=0.0;
2393                 virial->yy=0.0;
2394                 virial->zz=0.0;
2395                 virial->xy=0.0;
2396                 virial->xz=0.0;
2397                 virial->yz=0.0;
2398
2399                 /* reset site energy */
2400                 itom[i].e=0.0;
2401
2402         }
2403
2404         /* get energy, force and virial of every atom */
2405
2406         /* first (and only) loop over atoms i */
2407         for(i=0;i<count;i++) {
2408
2409                 /* single particle potential/force */
2410                 if(itom[i].attr&ATOM_ATTR_1BP)
2411                         if(moldyn->func1b)
2412                                 moldyn->func1b(moldyn,&(itom[i]));
2413
2414                 if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
2415                         continue;
2416
2417                 /* 2 body pair potential/force */
2418
2419                 link_cell_neighbour_index(moldyn,
2420                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
2421                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
2422                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
2423                                           neighbour_i);
2424
2425                 dnlc=lc->dnlc;
2426
2427                 /* first loop over atoms j */
2428                 if(moldyn->func2b) {
2429                         for(j=0;j<27;j++) {
2430
2431                                 bc_ij=(j<dnlc)?0:1;
2432 #ifdef STATIC_LISTS
2433                                 p=0;
2434
2435                                 while(neighbour_i[j][p]!=-1) {
2436
2437                                         jtom=&(atom[neighbour_i[j][p]]);
2438                                         p++;
2439 #elif LOWMEM_LISTS
2440                                 p=neighbour_i[j];
2441
2442                                 while(p!=-1) {
2443
2444                                         jtom=&(itom[p]);
2445                                         p=lc->subcell->list[p];
2446 #else
2447                                 this=&(neighbour_i[j]);
2448                                 list_reset_f(this);
2449
2450                                 if(this->start==NULL)
2451                                         continue;
2452
2453                                 do {
2454                                         jtom=this->current->data;
2455 #endif
2456
2457                                         if(jtom==&(itom[i]))
2458                                                 continue;
2459
2460                                         if((jtom->attr&ATOM_ATTR_2BP)&
2461                                            (itom[i].attr&ATOM_ATTR_2BP)) {
2462                                                 moldyn->func2b(moldyn,
2463                                                                &(itom[i]),
2464                                                                jtom,
2465                                                                bc_ij);
2466                                         }
2467 #ifdef STATIC_LISTS
2468                                 }
2469 #elif LOWMEM_LISTS
2470                                 }
2471 #else
2472                                 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2473 #endif
2474
2475                         }
2476                 }
2477
2478                 /* 3 body potential/force */
2479
2480                 if(!(itom[i].attr&ATOM_ATTR_3BP))
2481                         continue;
2482
2483                 /* copy the neighbour lists */
2484 #ifdef STATIC_LISTS
2485                 /* no copy needed for static lists */
2486 #elif LOWMEM_LISTS
2487                 /* no copy needed for lowmem lists */
2488 #else
2489                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
2490 #endif
2491
2492                 /* second loop over atoms j */
2493                 for(j=0;j<27;j++) {
2494
2495                         bc_ij=(j<dnlc)?0:1;
2496 #ifdef STATIC_LISTS
2497                         p=0;
2498
2499                         while(neighbour_i[j][p]!=-1) {
2500
2501                                 jtom=&(atom[neighbour_i[j][p]]);
2502                                 p++;
2503 #elif LOWMEM_LISTS
2504                                 p=neighbour_i[j];
2505
2506                                 while(p!=-1) {
2507
2508                                         jtom=&(itom[p]);
2509                                         p=lc->subcell->list[p];
2510 #else
2511                         this=&(neighbour_i[j]);
2512                         list_reset_f(this);
2513
2514                         if(this->start==NULL)
2515                                 continue;
2516
2517                         do {
2518
2519                                 jtom=this->current->data;
2520 #endif
2521
2522                                 if(jtom==&(itom[i]))
2523                                         continue;
2524
2525                                 if(!(jtom->attr&ATOM_ATTR_3BP))
2526                                         continue;
2527
2528                                 /* reset 3bp run */
2529                                 moldyn->run3bp=1;
2530
2531                                 if(moldyn->func3b_j1)
2532                                         moldyn->func3b_j1(moldyn,
2533                                                           &(itom[i]),
2534                                                           jtom,
2535                                                           bc_ij);
2536
2537                                 /* in first j loop, 3bp run can be skipped */
2538                                 if(!(moldyn->run3bp))
2539                                         continue;
2540
2541                                 /* first loop over atoms k */
2542                                 if(moldyn->func3b_k1) {
2543
2544                                 for(k=0;k<27;k++) {
2545
2546                                         bc_ik=(k<dnlc)?0:1;
2547 #ifdef STATIC_LISTS
2548                                         q=0;
2549
2550                                         while(neighbour_i[k][q]!=-1) {
2551
2552                                                 ktom=&(atom[neighbour_i[k][q]]);
2553                                                 q++;
2554 #elif LOWMEM_LISTS
2555                                         q=neighbour_i[k];
2556
2557                                         while(q!=-1) {
2558
2559                                                 ktom=&(itom[q]);
2560                                                 q=lc->subcell->list[q];
2561 #else
2562                                         that=&(neighbour_i2[k]);
2563                                         list_reset_f(that);
2564
2565                                         if(that->start==NULL)
2566                                                 continue;
2567
2568                                         do {
2569                                                 ktom=that->current->data;
2570 #endif
2571
2572                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
2573                                                         continue;
2574
2575                                                 if(ktom==jtom)
2576                                                         continue;
2577
2578                                                 if(ktom==&(itom[i]))
2579                                                         continue;
2580
2581                                                 moldyn->func3b_k1(moldyn,
2582                                                                   &(itom[i]),
2583                                                                   jtom,
2584                                                                   ktom,
2585                                                                   bc_ik|bc_ij);
2586
2587 #ifdef STATIC_LISTS
2588                                         }
2589 #elif LOWMEM_LISTS
2590                                         }
2591 #else
2592                                         } while(list_next_f(that)!=\
2593                                                 L_NO_NEXT_ELEMENT);
2594 #endif
2595
2596                                 }
2597
2598                                 }
2599
2600                                 if(moldyn->func3b_j2)
2601                                         moldyn->func3b_j2(moldyn,
2602                                                           &(itom[i]),
2603                                                           jtom,
2604                                                           bc_ij);
2605
2606                                 /* second loop over atoms k */
2607                                 if(moldyn->func3b_k2) {
2608
2609                                 for(k=0;k<27;k++) {
2610
2611                                         bc_ik=(k<dnlc)?0:1;
2612 #ifdef STATIC_LISTS
2613                                         q=0;
2614
2615                                         while(neighbour_i[k][q]!=-1) {
2616
2617                                                 ktom=&(atom[neighbour_i[k][q]]);
2618                                                 q++;
2619 #elif LOWMEM_LISTS
2620                                         q=neighbour_i[k];
2621
2622                                         while(q!=-1) {
2623
2624                                                 ktom=&(itom[q]);
2625                                                 q=lc->subcell->list[q];
2626 #else
2627                                         that=&(neighbour_i2[k]);
2628                                         list_reset_f(that);
2629
2630                                         if(that->start==NULL)
2631                                                 continue;
2632
2633                                         do {
2634                                                 ktom=that->current->data;
2635 #endif
2636
2637                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
2638                                                         continue;
2639
2640                                                 if(ktom==jtom)
2641                                                         continue;
2642
2643                                                 if(ktom==&(itom[i]))
2644                                                         continue;
2645
2646                                                 moldyn->func3b_k2(moldyn,
2647                                                                   &(itom[i]),
2648                                                                   jtom,
2649                                                                   ktom,
2650                                                                   bc_ik|bc_ij);
2651
2652 #ifdef STATIC_LISTS
2653                                         }
2654 #elif LOWMEM_LISTS
2655                                         }
2656 #else
2657                                         } while(list_next_f(that)!=\
2658                                                 L_NO_NEXT_ELEMENT);
2659 #endif
2660
2661                                 }
2662
2663                                 }
2664
2665                                 /* 2bp post function */
2666                                 if(moldyn->func3b_j3) {
2667                                         moldyn->func3b_j3(moldyn,
2668                                                           &(itom[i]),
2669                                                           jtom,bc_ij);
2670                                 }
2671 #ifdef STATIC_LISTS
2672                         }
2673 #elif LOWMEM_LISTS
2674                         }
2675 #else
2676                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
2677 #endif
2678
2679                 }
2680
2681 #ifdef DEBUG
2682         //printf("\n\n");
2683 #endif
2684 #ifdef VDEBUG
2685         printf("\n\n");
2686 #endif
2687
2688         }
2689
2690 #ifdef DEBUG
2691         //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
2692         if(moldyn->time>DSTART&&moldyn->time<DEND) {
2693                 printf("force:\n");
2694                 printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
2695                 printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
2696                 printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
2697         }
2698 #endif
2699
2700         /* some postprocessing */
2701 #ifdef PARALLEL
2702         #pragma omp parallel for
2703 #endif
2704         for(i=0;i<count;i++) {
2705                 /* calculate global virial */
2706                 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
2707                 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
2708                 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
2709                 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
2710                 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
2711                 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
2712
2713                 /* check forces regarding the given timestep */
2714                 if(v3_norm(&(itom[i].f))>\
2715                    0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
2716                         printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
2717                                i);
2718         }
2719
2720         return 0;
2721 }
2722
2723 /*
2724  * virial calculation
2725  */
2726
2727 //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
2728 int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
2729
2730         a->virial.xx+=f->x*d->x;
2731         a->virial.yy+=f->y*d->y;
2732         a->virial.zz+=f->z*d->z;
2733         a->virial.xy+=f->x*d->y;
2734         a->virial.xz+=f->x*d->z;
2735         a->virial.yz+=f->y*d->z;
2736
2737         return 0;
2738 }
2739
2740 /*
2741  * periodic boundary checking
2742  */
2743
2744 //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
2745 int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
2746
2747         double x,y,z;
2748         t_3dvec *dim;
2749
2750         dim=&(moldyn->dim);
2751
2752         x=dim->x/2;
2753         y=dim->y/2;
2754         z=dim->z/2;
2755
2756         if(moldyn->status&MOLDYN_STAT_PBX) {
2757                 if(a->x>=x) a->x-=dim->x;
2758                 else if(-a->x>x) a->x+=dim->x;
2759         }
2760         if(moldyn->status&MOLDYN_STAT_PBY) {
2761                 if(a->y>=y) a->y-=dim->y;
2762                 else if(-a->y>y) a->y+=dim->y;
2763         }
2764         if(moldyn->status&MOLDYN_STAT_PBZ) {
2765                 if(a->z>=z) a->z-=dim->z;
2766                 else if(-a->z>z) a->z+=dim->z;
2767         }
2768
2769         return 0;
2770 }
2771
2772 int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a) {
2773
2774         double x,y,z;
2775         t_3dvec *dim;
2776
2777         dim=&(moldyn->dim);
2778
2779         x=dim->x/2;
2780         y=dim->y/2;
2781         z=dim->z/2;
2782
2783         if(moldyn->status&MOLDYN_STAT_PBX) {
2784                 if(a->r.x>=x) {
2785                         a->pbc[0]+=1;
2786                         a->r.x-=dim->x;
2787                 }
2788                 else if(-a->r.x>x) {
2789                         a->pbc[0]-=1;
2790                         a->r.x+=dim->x;
2791                 }
2792         }
2793         if(moldyn->status&MOLDYN_STAT_PBY) {
2794                 if(a->r.y>=y) {
2795                         a->pbc[1]+=1;
2796                         a->r.y-=dim->y;
2797                 }
2798                 else if(-a->r.y>y) {
2799                         a->pbc[1]-=1;
2800                         a->r.y+=dim->y;
2801                 }
2802         }
2803         if(moldyn->status&MOLDYN_STAT_PBZ) {
2804                 if(a->r.z>=z) {
2805                         a->pbc[2]+=1;
2806                         a->r.z-=dim->z;
2807                 }
2808                 else if(-a->r.z>z) {
2809                         a->pbc[2]-=1;
2810                         a->r.z+=dim->z;
2811                 }
2812         }
2813
2814         return 0;
2815 }
2816
2817 /*
2818  * debugging / critical check functions
2819  */
2820
2821 int moldyn_bc_check(t_moldyn *moldyn) {
2822
2823         t_atom *atom;
2824         t_3dvec *dim;
2825         int i;
2826         double x;
2827         u8 byte;
2828         int j,k;
2829
2830         atom=moldyn->atom;
2831         dim=&(moldyn->dim);
2832         x=dim->x/2;
2833
2834         for(i=0;i<moldyn->count;i++) {
2835                 if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {
2836                         printf("FATAL: atom %d: x: %.20f (%.20f)\n",
2837                                i,atom[i].r.x,dim->x/2);
2838                         printf("diagnostic:\n");
2839                         printf("-----------\natom.r.x:\n");
2840                         for(j=0;j<8;j++) {
2841                                 memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
2842                                 for(k=0;k<8;k++)
2843                                         printf("%d%c",
2844                                         ((byte)&(1<<k))?1:0,
2845                                         (k==7)?'\n':'|');
2846                         }
2847                         printf("---------------\nx=dim.x/2:\n");
2848                         for(j=0;j<8;j++) {
2849                                 memcpy(&byte,(u8 *)(&x)+j,1);
2850                                 for(k=0;k<8;k++)
2851                                         printf("%d%c",
2852                                         ((byte)&(1<<k))?1:0,
2853                                         (k==7)?'\n':'|');
2854                         }
2855                         if(atom[i].r.x==x) printf("the same!\n");
2856                         else printf("different!\n");
2857                 }
2858                 if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
2859                         printf("FATAL: atom %d: y: %.20f (%.20f)\n",
2860                                i,atom[i].r.y,dim->y/2);
2861                 if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
2862                         printf("FATAL: atom %d: z: %.20f (%.20f)\n",
2863                                i,atom[i].r.z,dim->z/2);
2864         }
2865
2866         return 0;
2867 }
2868
2869 /*
2870  * restore function
2871  */
2872
2873 int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
2874
2875         int fd;
2876         int cnt,size;
2877         int fsize;
2878         int corr;
2879
2880         fd=open(file,O_RDONLY);
2881         if(fd<0) {
2882                 perror("[moldyn] load save file open");
2883                 return fd;
2884         }
2885
2886         fsize=lseek(fd,0,SEEK_END);
2887         lseek(fd,0,SEEK_SET);
2888
2889         size=sizeof(t_moldyn);
2890
2891         while(size) {
2892                 cnt=read(fd,moldyn,size);
2893                 if(cnt<0) {
2894                         perror("[moldyn] load save file read (moldyn)");
2895                         return cnt;
2896                 }
2897                 size-=cnt;
2898         }
2899
2900         size=moldyn->count*sizeof(t_atom);
2901
2902         /* correcting possible atom data offset */
2903         corr=0;
2904         if(fsize!=sizeof(t_moldyn)+size) {
2905                 corr=fsize-sizeof(t_moldyn)-size;
2906                 printf("[moldyn] WARNING: lsf (illegal file size)\n");
2907                 printf("  modifying offset:\n");
2908                 printf("  - current pos: %d\n",sizeof(t_moldyn));
2909                 printf("  - atom size: %d\n",size);
2910                 printf("  - file size: %d\n",fsize);
2911                 printf("  => correction: %d\n",corr);
2912                 lseek(fd,corr,SEEK_CUR);
2913         }
2914
2915         moldyn->atom=(t_atom *)malloc(size);
2916         if(moldyn->atom==NULL) {
2917                 perror("[moldyn] load save file malloc (atoms)");
2918                 return -1;
2919         }
2920
2921         while(size) {
2922                 cnt=read(fd,moldyn->atom,size);
2923                 if(cnt<0) {
2924                         perror("[moldyn] load save file read (atoms)");
2925                         return cnt;
2926                 }
2927                 size-=cnt;
2928         }
2929
2930 #ifdef PTHREADS
2931         amutex=malloc(moldyn->count*sizeof(pthread_mutex_t));
2932         if(amutex==NULL) {
2933                 perror("[moldyn] load save file (mutexes)");
2934                 return -1;
2935         }
2936         for(cnt=0;cnt<moldyn->count;cnt++)
2937                 pthread_mutex_init(&(amutex[cnt]),NULL);
2938 #endif
2939
2940         // hooks etc ...
2941
2942         return 0;
2943 }
2944
2945 int moldyn_free_save_file(t_moldyn *moldyn) {
2946
2947         free(moldyn->atom);
2948
2949         return 0;
2950 }
2951
2952 int moldyn_load(t_moldyn *moldyn) {
2953
2954         // later ...
2955
2956         return 0;
2957 }
2958
2959 /*
2960  * function to find/callback all combinations of 2 body bonds
2961  */
2962
2963 int process_2b_bonds(t_moldyn *moldyn,void *data,
2964                      int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
2965                                     void *data,u8 bc)) {
2966
2967         t_linkcell *lc;
2968 #ifdef STATIC_LISTS
2969         int *neighbour[27];
2970         int p;
2971 #elif LOWMEM_LISTS
2972         int neighbour[27];
2973         int p;
2974 #else
2975         t_list neighbour[27];
2976         t_list *this;
2977 #endif
2978         u8 bc;
2979         t_atom *itom,*jtom;
2980         int i,j;
2981
2982         lc=&(moldyn->lc);
2983         itom=moldyn->atom;
2984
2985         for(i=0;i<moldyn->count;i++) {
2986                 /* neighbour indexing */
2987                 link_cell_neighbour_index(moldyn,
2988                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
2989                                           (itom[i].r.y+moldyn->dim.y/2)/lc->x,
2990                                           (itom[i].r.z+moldyn->dim.z/2)/lc->x,
2991                                           neighbour);
2992
2993                 for(j=0;j<27;j++) {
2994
2995                         bc=(j<lc->dnlc)?0:1;
2996
2997 #ifdef STATIC_LISTS
2998                         p=0;
2999
3000                         while(neighbour[j][p]!=-1) {
3001
3002                                 jtom=&(moldyn->atom[neighbour[j][p]]);
3003                                 p++;
3004 #elif LOWMEM_LISTS
3005                         p=neighbour[j];
3006
3007                         while(p!=-1) {
3008
3009                                 jtom=&(itom[p]);
3010                                 p=lc->subcell->list[p];
3011 #else
3012                         this=&(neighbour[j]);
3013                         list_reset_f(this);
3014
3015                         if(this->start==NULL)
3016                                 continue;
3017
3018                         do {
3019
3020                                 jtom=this->current->data;
3021 #endif
3022
3023                                 /* process bond */
3024                                 process(moldyn,&(itom[i]),jtom,data,bc);
3025
3026 #ifdef STATIC_LISTS
3027                         }
3028 #elif LOWMEM_LISTS
3029                         }
3030 #else
3031                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
3032 #endif
3033                 }
3034         }
3035
3036         return 0;
3037
3038 }
3039
3040 /*
3041  * function to find neighboured atoms
3042  */
3043
3044 int process_neighbours(t_moldyn *moldyn,void *data,t_atom *atom,
3045                        int (*process)(t_moldyn *moldyn,t_atom *atom,t_atom *natom,
3046                                       void *data,u8 bc)) {
3047
3048         t_linkcell *lc;
3049 #ifdef STATIC_LISTS
3050         int *neighbour[27];
3051         int p;
3052 #elif LOWMEM_LISTS
3053         int neighbour[27];
3054         int p;
3055 #else
3056         t_list neighbour[27];
3057         t_list *this;
3058 #endif
3059         u8 bc;
3060         t_atom *natom;
3061         int j;
3062
3063         lc=&(moldyn->lc);
3064
3065         /* neighbour indexing */
3066         link_cell_neighbour_index(moldyn,
3067                                   (atom->r.x+moldyn->dim.x/2)/lc->x,
3068                                   (atom->r.y+moldyn->dim.y/2)/lc->x,
3069                                   (atom->r.z+moldyn->dim.z/2)/lc->x,
3070                                   neighbour);
3071
3072         for(j=0;j<27;j++) {
3073
3074                 bc=(j<lc->dnlc)?0:1;
3075
3076 #ifdef STATIC_LISTS
3077                 p=0;
3078
3079                 while(neighbour[j][p]!=-1) {
3080
3081                         natom=&(moldyn->atom[neighbour[j][p]]);
3082                         p++;
3083 #elif LOWMEM_LISTS
3084                 p=neighbour[j];
3085
3086                 while(p!=-1) {
3087
3088                         natom=&(moldyn->atom[p]);
3089                         p=lc->subcell->list[p];
3090 #else
3091                 this=&(neighbour[j]);
3092                 list_reset_f(this);
3093
3094                 if(this->start==NULL)
3095                         continue;
3096
3097                 do {
3098
3099                         natom=this->current->data;
3100 #endif
3101
3102                         /* process bond */
3103                         process(moldyn,atom,natom,data,bc);
3104
3105 #ifdef STATIC_LISTS
3106                 }
3107 #elif LOWMEM_LISTS
3108                 }
3109 #else
3110                 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
3111 #endif
3112         }
3113
3114         return 0;
3115
3116 }
3117
3118 /*
3119  * post processing functions
3120  */
3121
3122 int get_line(int fd,char *line,int max) {
3123
3124         int count,ret;
3125
3126         count=0;
3127
3128         while(1) {
3129                 if(count==max) return count;
3130                 ret=read(fd,line+count,1);
3131                 if(ret<=0) return ret;
3132                 if(line[count]=='\n') {
3133                         memset(line+count,0,max-count-1);
3134                         //line[count]='\0';
3135                         return count+1;
3136                 }
3137                 count+=1;
3138         }
3139 }
3140
3141 int pair_correlation_init(t_moldyn *moldyn,double dr) {
3142
3143
3144         return 0;
3145 }
3146
3147 int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
3148
3149         int i;
3150         t_atom *atom;
3151         t_3dvec dist;
3152         t_3dvec final_r;
3153         double d2;
3154         int a_cnt;
3155         int b_cnt;
3156
3157         atom=moldyn->atom;
3158         dc[0]=0;
3159         dc[1]=0;
3160         dc[2]=0;
3161         a_cnt=0;
3162         b_cnt=0;
3163
3164         for(i=0;i<moldyn->count;i++) {
3165
3166                 /* care for pb crossing */
3167                 final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x;
3168                 final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y;
3169                 final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z;
3170                 /* calculate distance */
3171                 v3_sub(&dist,&final_r,&(atom[i].r_0));
3172                 d2=v3_absolute_square(&dist);
3173
3174                 if(atom[i].brand) {
3175                         b_cnt+=1;
3176                         dc[1]+=d2;
3177                 }
3178                 else {
3179                         a_cnt+=1;
3180                         dc[0]+=d2;
3181                 }
3182
3183                 dc[2]+=d2;
3184         }
3185
3186         dc[0]*=(1.0/(6.0*moldyn->time*a_cnt));
3187         dc[1]*=(1.0/(6.0*moldyn->time*b_cnt));
3188         dc[2]*=(1.0/(6.0*moldyn->time*moldyn->count));
3189
3190         return 0;
3191 }
3192
3193 int bonding_analyze(t_moldyn *moldyn,double *cnt) {
3194
3195         return 0;
3196 }
3197
3198 int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
3199                                        t_atom *jtom,void *data,u8 bc) {
3200
3201         t_3dvec dist;
3202         double d;
3203         int s;
3204         t_pcc *pcc;
3205
3206         /* only count pairs once,
3207          * skip same atoms */
3208         if(itom->tag>=jtom->tag)
3209                 return 0;
3210
3211         /*
3212          * pair correlation calc
3213          */
3214
3215         /* get pcc data */
3216         pcc=data;
3217
3218         /* distance */
3219         v3_sub(&dist,&(jtom->r),&(itom->r));
3220         if(bc) check_per_bound(moldyn,&dist);
3221         d=v3_absolute_square(&dist);
3222
3223         /* ignore if greater cutoff */
3224         if(d>moldyn->cutoff_square)
3225                 return 0;
3226
3227         /* fill the slots */
3228         d=sqrt(d);
3229         s=(int)(d/pcc->dr);
3230
3231         /* should never happen but it does 8) -
3232          * related to -ffloat-store problem! */
3233         if(s>=pcc->o1) {
3234                 printf("[moldyn] WARNING: pcc (%d/%d)",
3235                        s,pcc->o1);
3236                 printf("\n");
3237                 s=pcc->o1-1;
3238         }
3239
3240         if(itom->brand!=jtom->brand) {
3241                 /* mixed */
3242                 pcc->stat[s]+=1;
3243         }
3244         else {
3245                 /* type a - type a bonds */
3246                 if(itom->brand==0)
3247                         pcc->stat[s+pcc->o1]+=1;
3248                 else
3249                 /* type b - type b bonds */
3250                         pcc->stat[s+pcc->o2]+=1;
3251         }
3252
3253         return 0;
3254 }
3255
3256 int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
3257
3258         t_pcc pcc;
3259         double norm;
3260         int i;
3261
3262         pcc.dr=dr;
3263         pcc.o1=moldyn->cutoff/dr;
3264         pcc.o2=2*pcc.o1;
3265
3266         if(pcc.o1*dr<=moldyn->cutoff)
3267                 printf("[moldyn] WARNING: pcc (low #slots)\n");
3268
3269         printf("[moldyn] pair correlation calc info:\n");
3270         printf("  time: %f\n",moldyn->time);
3271         printf("  count: %d\n",moldyn->count);
3272         printf("  cutoff: %f\n",moldyn->cutoff);
3273         printf("  temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg);
3274
3275         if(ptr!=NULL) {
3276                 pcc.stat=(double *)ptr;
3277         }
3278         else {
3279                 pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double));
3280                 if(pcc.stat==NULL) {
3281                         perror("[moldyn] pair correlation malloc");
3282                         return -1;
3283                 }
3284         }
3285
3286         memset(pcc.stat,0,3*pcc.o1*sizeof(double));
3287
3288         /* process */
3289         process_2b_bonds(moldyn,&pcc,calculate_pair_correlation_process);
3290
3291         /* normalization */
3292         for(i=1;i<pcc.o1;i++) {
3293                  // normalization: 4 pi r^2 dr
3294                  // here: not double counting pairs -> 2 pi r r dr
3295                  // ... and actually it's a constant times r^2
3296                 norm=i*i*dr*dr;
3297                 pcc.stat[i]/=norm;
3298                 pcc.stat[pcc.o1+i]/=norm;
3299                 pcc.stat[pcc.o2+i]/=norm;
3300         }
3301         /* */
3302
3303         if(ptr==NULL) {
3304                 /* todo: store/print pair correlation function */
3305                 free(pcc.stat);
3306         }
3307
3308         return 0;
3309 }
3310
3311 int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
3312                          void *data,u8 bc) {
3313
3314         t_ba *ba;
3315         t_3dvec dist;
3316         double d;
3317
3318         if(itom->tag>=jtom->tag)
3319                 return 0;
3320
3321         /* distance */
3322         v3_sub(&dist,&(jtom->r),&(itom->r));
3323         if(bc) check_per_bound(moldyn,&dist);
3324         d=v3_absolute_square(&dist);
3325
3326         /* ignore if greater or equal cutoff */
3327         if(d>moldyn->cutoff_square)
3328                 return 0;
3329
3330         /* check for potential bond */
3331         if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
3332                 return 0;
3333
3334         /* now count this bonding ... */
3335         ba=data;
3336
3337         /* increase total bond counter
3338          */
3339         ba->tcnt+=1;
3340
3341         if(itom->brand==0)
3342                 ba->acnt[jtom->tag]+=1;
3343         else
3344                 ba->bcnt[jtom->tag]+=1;
3345
3346         if(jtom->brand==0)
3347                 ba->acnt[itom->tag]+=1;
3348         else
3349                 ba->bcnt[itom->tag]+=1;
3350
3351         return 0;
3352 }
3353
3354 int bond_analyze(t_moldyn *moldyn,double *quality) {
3355
3356         int qcnt4;
3357         int qcnt3;
3358         int ccnt4;
3359         int ccnt3;
3360         int bcnt;
3361         t_ba ba;
3362         int i;
3363         t_atom *atom;
3364
3365         ba.acnt=malloc(moldyn->count*sizeof(int));
3366         if(ba.acnt==NULL) {
3367                 perror("[moldyn] bond analyze malloc (a)");
3368                 return -1;
3369         }
3370         memset(ba.acnt,0,moldyn->count*sizeof(int));
3371
3372         ba.bcnt=malloc(moldyn->count*sizeof(int));
3373         if(ba.bcnt==NULL) {
3374                 perror("[moldyn] bond analyze malloc (b)");
3375                 return -1;
3376         }
3377         memset(ba.bcnt,0,moldyn->count*sizeof(int));
3378
3379         ba.tcnt=0;
3380         qcnt4=0; qcnt3=0;
3381         ccnt4=0; ccnt3=0;
3382         bcnt=0;
3383
3384         atom=moldyn->atom;
3385
3386         process_2b_bonds(moldyn,&ba,bond_analyze_process);
3387
3388         for(i=0;i<moldyn->count;i++) {
3389                 if(atom[i].brand==0) {
3390                         if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
3391                                 qcnt4+=4;
3392                         if((ba.acnt[i]==0)&(ba.bcnt[i]==3))
3393                                 qcnt3+=3;
3394                 }
3395                 else {
3396                         if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
3397                                 qcnt4+=4;
3398                                 ccnt4+=1;
3399                         }
3400                         if((ba.acnt[i]==3)&(ba.bcnt[i]==0)) {
3401                                 qcnt3+=4;
3402                                 ccnt3+=1;
3403                         }
3404                         bcnt+=1;
3405                 }
3406         }
3407
3408         if(quality) {
3409                 quality[0]=1.0*ccnt4/bcnt;
3410                 quality[1]=1.0*ccnt3/bcnt;
3411         }
3412         else {
3413                 printf("[moldyn] bond analyze: %f %f\n",
3414                        1.0*ccnt4/bcnt,1.0*ccnt3/bcnt);
3415         }
3416
3417         return 0;
3418 }
3419
3420 /*
3421  * visualization code
3422  */
3423
3424 int visual_init(t_moldyn *moldyn,char *filebase) {
3425
3426         strncpy(moldyn->vis.fb,filebase,128);
3427
3428         return 0;
3429 }
3430
3431 int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
3432                          void *data,u8 bc) {
3433
3434         t_vb *vb;
3435
3436         vb=data;
3437
3438         if(itom->tag>=jtom->tag)
3439                 return 0;
3440
3441         if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
3442                 return 0;
3443
3444         if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB))
3445                 dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n",
3446                         itom->r.x,itom->r.y,itom->r.z,
3447                         jtom->r.x,jtom->r.y,jtom->r.z);
3448
3449         return 0;
3450 }
3451
3452 #ifdef VISUAL_THREAD
3453 void *visual_atoms(void *ptr) {
3454 #else
3455 int visual_atoms(t_moldyn *moldyn) {
3456 #endif
3457
3458         int i;
3459         char file[128+64];
3460         t_3dvec dim;
3461         double help;
3462         t_visual *v;
3463         t_atom *atom;
3464         t_vb vb;
3465         t_3dvec strain;
3466 #ifdef VISUAL_THREAD
3467         t_moldyn *moldyn;
3468
3469         moldyn=ptr;
3470 #endif
3471
3472         v=&(moldyn->vis);
3473         dim.x=v->dim.x;
3474         dim.y=v->dim.y;
3475         dim.z=v->dim.z;
3476         atom=moldyn->atom;
3477
3478         help=(dim.x+dim.y);
3479
3480         sprintf(file,"%s/atomic_conf_%08.f.xyz",v->fb,moldyn->time);
3481         vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
3482         if(vb.fd<0) {
3483                 perror("open visual save file fd");
3484 #ifndef VISUAL_THREAD
3485                 return -1;
3486 #endif
3487         }
3488
3489         /* write the actual data file */
3490
3491         // povray header
3492         dprintf(vb.fd,"# [P] %d %08.f <%f,%f,%f>\n",
3493                 moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
3494
3495         // atomic configuration
3496         for(i=0;i<moldyn->count;i++) {
3497                 v3_sub(&strain,&(atom[i].r),&(atom[i].r_0));
3498                 check_per_bound(moldyn,&strain);
3499                 // atom type, positions, color and kinetic energy
3500                 dprintf(vb.fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
3501                                                     atom[i].r.x,
3502                                                     atom[i].r.y,
3503                                                     atom[i].r.z,
3504                                                     pse_col[atom[i].element],
3505                                                     //atom[i].ekin);
3506                                                     sqrt(v3_absolute_square(&strain)));
3507         }
3508
3509         // bonds between atoms
3510 #ifndef VISUAL_THREAD
3511         process_2b_bonds(moldyn,&vb,visual_bonds_process);
3512 #endif
3513
3514         // boundaries
3515         if(dim.x) {
3516                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3517                         -dim.x/2,-dim.y/2,-dim.z/2,
3518                         dim.x/2,-dim.y/2,-dim.z/2);
3519                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3520                         -dim.x/2,-dim.y/2,-dim.z/2,
3521                         -dim.x/2,dim.y/2,-dim.z/2);
3522                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3523                         dim.x/2,dim.y/2,-dim.z/2,
3524                         dim.x/2,-dim.y/2,-dim.z/2);
3525                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3526                         -dim.x/2,dim.y/2,-dim.z/2,
3527                         dim.x/2,dim.y/2,-dim.z/2);
3528
3529                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3530                         -dim.x/2,-dim.y/2,dim.z/2,
3531                         dim.x/2,-dim.y/2,dim.z/2);
3532                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3533                         -dim.x/2,-dim.y/2,dim.z/2,
3534                         -dim.x/2,dim.y/2,dim.z/2);
3535                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3536                         dim.x/2,dim.y/2,dim.z/2,
3537                         dim.x/2,-dim.y/2,dim.z/2);
3538                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3539                         -dim.x/2,dim.y/2,dim.z/2,
3540                         dim.x/2,dim.y/2,dim.z/2);
3541
3542                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3543                         -dim.x/2,-dim.y/2,dim.z/2,
3544                         -dim.x/2,-dim.y/2,-dim.z/2);
3545                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3546                         -dim.x/2,dim.y/2,dim.z/2,
3547                         -dim.x/2,dim.y/2,-dim.z/2);
3548                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3549                         dim.x/2,-dim.y/2,dim.z/2,
3550                         dim.x/2,-dim.y/2,-dim.z/2);
3551                 dprintf(vb.fd,"# [D] %f %f %f %f %f %f\n",
3552                         dim.x/2,dim.y/2,dim.z/2,
3553                         dim.x/2,dim.y/2,-dim.z/2);
3554         }
3555
3556         close(vb.fd);
3557
3558 #ifdef VISUAL_THREAD
3559         pthread_exit(NULL);
3560
3561 }
3562 #else
3563
3564         return 0;
3565 }
3566 #endif
3567
3568 /*
3569  * fpu cntrol functions
3570  */
3571
3572 // set rounding to double (eliminates -ffloat-store!)
3573 int fpu_set_rtd(void) {
3574
3575         fpu_control_t ctrl;
3576
3577         _FPU_GETCW(ctrl);
3578
3579         ctrl&=~_FPU_EXTENDED;
3580         ctrl|=_FPU_DOUBLE;
3581
3582         _FPU_SETCW(ctrl);
3583
3584         return 0;
3585 }
3586