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