moldyn.c

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