moldyn.c

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