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