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