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