moldyn.c

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