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