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