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