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