moldyn.c

   1 /*
   2  * moldyn.c - molecular dynamics library main file
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <sys/types.h>
  13 #include <sys/stat.h>
  14 #include <fcntl.h>
  15 #include <unistd.h>
  16 #include <sys/time.h>
  17 #include <time.h>
  18 #include <math.h>
  19
  20 #include "moldyn.h"
  21 #include "report/report.h"
  22
  23 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
  24
  25         printf("[moldyn] init\n");
  26
  27         memset(moldyn,0,sizeof(t_moldyn));
  28
  29         moldyn->argc=argc;
  30         moldyn->args=argv;
  31
  32         rand_init(&(moldyn->random),NULL,1);
  33         moldyn->random.status|=RAND_STAT_VERBOSE;
  34
  35         return 0;
  36 }
  37
  38 int moldyn_shutdown(t_moldyn *moldyn) {
  39
  40         printf("[moldyn] shutdown\n");
  41
  42         moldyn_log_shutdown(moldyn);
  43         link_cell_shutdown(moldyn);
  44         rand_close(&(moldyn->random));
  45         free(moldyn->atom);
  46
  47         return 0;
  48 }
  49
  50 int set_int_alg(t_moldyn *moldyn,u8 algo) {
  51
  52         printf("[moldyn] integration algorithm: ");
  53
  54         switch(algo) {
  55                 case MOLDYN_INTEGRATE_VERLET:
  56                         moldyn->integrate=velocity_verlet;
  57                         printf("velocity verlet\n");
  58                         break;
  59                 default:
  60                         printf("unknown integration algorithm: %02x\n",algo);
  61                         printf("unknown\n");
  62                         return -1;
  63         }
  64
  65         return 0;
  66 }
  67
  68 int set_cutoff(t_moldyn *moldyn,double cutoff) {
  69
  70         moldyn->cutoff=cutoff;
  71
  72         printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff);
  73
  74         return 0;
  75 }
  76
  77 int set_temperature(t_moldyn *moldyn,double t_ref) {
  78
  79         moldyn->t_ref=t_ref;
  80
  81         printf("[moldyn] temperature [K]: %f\n",moldyn->t_ref);
  82
  83         return 0;
  84 }
  85
  86 int set_pressure(t_moldyn *moldyn,double p_ref) {
  87
  88         moldyn->p_ref=p_ref;
  89
  90         printf("[moldyn] pressure [bar]: %f\n",moldyn->p_ref/BAR);
  91
  92         return 0;
  93 }
  94
  95 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
  96
  97         moldyn->pt_scale=(ptype|ttype);
  98         moldyn->t_tc=ttc;
  99         moldyn->p_tc=ptc;
 100
 101         printf("[moldyn] p/t scaling:\n");
 102
 103         printf("  p: %s",ptype?"yes":"no ");
 104         if(ptype)
 105                 printf(" | type: %02x | factor: %f",ptype,ptc);
 106         printf("\n");
 107
 108         printf("  t: %s",ttype?"yes":"no ");
 109         if(ttype)
 110                 printf(" | type: %02x | factor: %f",ttype,ttc);
 111         printf("\n");
 112
 113         return 0;
 114 }
 115
 116 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
 117
 118         moldyn->dim.x=x;
 119         moldyn->dim.y=y;
 120         moldyn->dim.z=z;
 121
 122         moldyn->volume=x*y*z;
 123
 124         if(visualize) {
 125                 moldyn->vis.dim.x=x;
 126                 moldyn->vis.dim.y=y;
 127                 moldyn->vis.dim.z=z;
 128         }
 129
 130         moldyn->dv=0.000001*moldyn->volume;
 131
 132         printf("[moldyn] dimensions in A and A^3 respectively:\n");
 133         printf("  x: %f\n",moldyn->dim.x);
 134         printf("  y: %f\n",moldyn->dim.y);
 135         printf("  z: %f\n",moldyn->dim.z);
 136         printf("  volume: %f\n",moldyn->volume);
 137         printf("  visualize simulation box: %s\n",visualize?"yes":"no");
 138         printf("  delta volume (pressure calc): %f\n",moldyn->dv);
 139
 140         return 0;
 141 }
 142
 143 int set_nn_dist(t_moldyn *moldyn,double dist) {
 144
 145         moldyn->nnd=dist;
 146
 147         return 0;
 148 }
 149
 150 int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
 151
 152         printf("[moldyn] periodic boundary conditions:\n");
 153
 154         if(x)
 155                 moldyn->status|=MOLDYN_STAT_PBX;
 156
 157         if(y)
 158                 moldyn->status|=MOLDYN_STAT_PBY;
 159
 160         if(z)
 161                 moldyn->status|=MOLDYN_STAT_PBZ;
 162
 163         printf("  x: %s\n",x?"yes":"no");
 164         printf("  y: %s\n",y?"yes":"no");
 165         printf("  z: %s\n",z?"yes":"no");
 166
 167         return 0;
 168 }
 169
 170 int set_potential1b(t_moldyn *moldyn,pf_func1b func) {
 171
 172         moldyn->func1b=func;
 173
 174         return 0;
 175 }
 176
 177 int set_potential2b(t_moldyn *moldyn,pf_func2b func) {
 178
 179         moldyn->func2b=func;
 180
 181         return 0;
 182 }
 183
 184 int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) {
 185
 186         moldyn->func3b_j1=func;
 187
 188         return 0;
 189 }
 190
 191 int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) {
 192
 193         moldyn->func3b_j2=func;
 194
 195         return 0;
 196 }
 197
 198 int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) {
 199
 200         moldyn->func3b_j3=func;
 201
 202         return 0;
 203 }
 204
 205 int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) {
 206
 207         moldyn->func3b_k1=func;
 208
 209         return 0;
 210 }
 211
 212 int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) {
 213
 214         moldyn->func3b_k2=func;
 215
 216         return 0;
 217 }
 218
 219 int set_potential_params(t_moldyn *moldyn,void *params) {
 220
 221         moldyn->pot_params=params;
 222
 223         return 0;
 224 }
 225
 226 int set_avg_skip(t_moldyn *moldyn,int skip) {
 227
 228         printf("[moldyn] skip %d steps before starting average calc\n",skip);
 229         moldyn->avg_skip=skip;
 230
 231         return 0;
 232 }
 233
 234 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
 235
 236         strncpy(moldyn->vlsdir,dir,127);
 237
 238         return 0;
 239 }
 240
 241 int moldyn_set_report(t_moldyn *moldyn,char *author,char *title) {
 242
 243         strncpy(moldyn->rauthor,author,63);
 244         strncpy(moldyn->rtitle,title,63);
 245
 246         return 0;
 247 }
 248
 249 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 250
 251         char filename[128];
 252         int ret;
 253
 254         printf("[moldyn] set log: ");
 255
 256         switch(type) {
 257                 case LOG_TOTAL_ENERGY:
 258                         moldyn->ewrite=timer;
 259                         snprintf(filename,127,"%s/energy",moldyn->vlsdir);
 260                         moldyn->efd=open(filename,
 261                                          O_WRONLY|O_CREAT|O_EXCL,
 262                                          S_IRUSR|S_IWUSR);
 263                         if(moldyn->efd<0) {
 264                                 perror("[moldyn] energy log fd open");
 265                                 return moldyn->efd;
 266                         }
 267                         dprintf(moldyn->efd,"# total energy log file\n");
 268                         printf("total energy (%d)\n",timer);
 269                         break;
 270                 case LOG_TOTAL_MOMENTUM:
 271                         moldyn->mwrite=timer;
 272                         snprintf(filename,127,"%s/momentum",moldyn->vlsdir);
 273                         moldyn->mfd=open(filename,
 274                                          O_WRONLY|O_CREAT|O_EXCL,
 275                                          S_IRUSR|S_IWUSR);
 276                         if(moldyn->mfd<0) {
 277                                 perror("[moldyn] momentum log fd open");
 278                                 return moldyn->mfd;
 279                         }
 280                         dprintf(moldyn->efd,"# total momentum log file\n");
 281                         printf("total momentum (%d)\n",timer);
 282                         break;
 283                 case LOG_PRESSURE:
 284                         moldyn->pwrite=timer;
 285                         snprintf(filename,127,"%s/pressure",moldyn->vlsdir);
 286                         moldyn->pfd=open(filename,
 287                                          O_WRONLY|O_CREAT|O_EXCL,
 288                                          S_IRUSR|S_IWUSR);
 289                         if(moldyn->pfd<0) {
 290                                 perror("[moldyn] pressure log file\n");
 291                                 return moldyn->pfd;
 292                         }
 293                         dprintf(moldyn->pfd,"# pressure log file\n");
 294                         printf("pressure (%d)\n",timer);
 295                         break;
 296                 case LOG_TEMPERATURE:
 297                         moldyn->twrite=timer;
 298                         snprintf(filename,127,"%s/temperature",moldyn->vlsdir);
 299                         moldyn->tfd=open(filename,
 300                                          O_WRONLY|O_CREAT|O_EXCL,
 301                                          S_IRUSR|S_IWUSR);
 302                         if(moldyn->tfd<0) {
 303                                 perror("[moldyn] temperature log file\n");
 304                                 return moldyn->tfd;
 305                         }
 306                         dprintf(moldyn->tfd,"# temperature log file\n");
 307                         printf("temperature (%d)\n",timer);
 308                         break;
 309                 case SAVE_STEP:
 310                         moldyn->swrite=timer;
 311                         printf("save file (%d)\n",timer);
 312                         break;
 313                 case VISUAL_STEP:
 314                         moldyn->vwrite=timer;
 315                         ret=visual_init(&(moldyn->vis),moldyn->vlsdir);
 316                         if(ret<0) {
 317                                 printf("[moldyn] visual init failure\n");
 318                                 return ret;
 319                         }
 320                         printf("visual file (%d)\n",timer);
 321                         break;
 322                 case CREATE_REPORT:
 323                         snprintf(filename,127,"%s/report.tex",moldyn->vlsdir);
 324                         moldyn->rfd=open(filename,
 325                                          O_WRONLY|O_CREAT|O_EXCL,
 326                                          S_IRUSR|S_IWUSR);
 327                         if(moldyn->rfd<0) {
 328                                 perror("[moldyn] report fd open");
 329                                 return moldyn->rfd;
 330                         }
 331                         printf("report -> ");
 332                         if(moldyn->efd) {
 333                                 snprintf(filename,127,"%s/e_plot.scr",
 334                                          moldyn->vlsdir);
 335                                 moldyn->epfd=open(filename,
 336                                                  O_WRONLY|O_CREAT|O_EXCL,
 337                                                  S_IRUSR|S_IWUSR);
 338                                 if(moldyn->epfd<0) {
 339                                         perror("[moldyn] energy plot fd open");
 340                                         return moldyn->epfd;
 341                                 }
 342                                 dprintf(moldyn->epfd,e_plot_script);
 343                                 close(moldyn->epfd);
 344                                 printf("energy ");
 345                         }
 346                         if(moldyn->pfd) {
 347                                 snprintf(filename,127,"%s/pressure_plot.scr",
 348                                          moldyn->vlsdir);
 349                                 moldyn->ppfd=open(filename,
 350                                                   O_WRONLY|O_CREAT|O_EXCL,
 351                                                   S_IRUSR|S_IWUSR);
 352                                 if(moldyn->ppfd<0) {
 353                                         perror("[moldyn] p plot fd open");
 354                                         return moldyn->ppfd;
 355                                 }
 356                                 dprintf(moldyn->ppfd,pressure_plot_script);
 357                                 close(moldyn->ppfd);
 358                                 printf("pressure ");
 359                         }
 360                         if(moldyn->tfd) {
 361                                 snprintf(filename,127,"%s/temperature_plot.scr",
 362                                          moldyn->vlsdir);
 363                                 moldyn->tpfd=open(filename,
 364                                                   O_WRONLY|O_CREAT|O_EXCL,
 365                                                   S_IRUSR|S_IWUSR);
 366                                 if(moldyn->tpfd<0) {
 367                                         perror("[moldyn] t plot fd open");
 368                                         return moldyn->tpfd;
 369                                 }
 370                                 dprintf(moldyn->tpfd,temperature_plot_script);
 371                                 close(moldyn->tpfd);
 372                                 printf("temperature ");
 373                         }
 374                         dprintf(moldyn->rfd,report_start,
 375                                 moldyn->rauthor,moldyn->rtitle);
 376                         printf("\n");
 377                         break;
 378                 default:
 379                         printf("unknown log type: %02x\n",type);
 380                         return -1;
 381         }
 382
 383         return 0;
 384 }
 385
 386 int moldyn_log_shutdown(t_moldyn *moldyn) {
 387
 388         char sc[256];
 389
 390         printf("[moldyn] log shutdown\n");
 391         if(moldyn->efd) {
 392                 close(moldyn->efd);
 393                 if(moldyn->rfd) {
 394                         dprintf(moldyn->rfd,report_energy);
 395                         snprintf(sc,255,"cd %s && gnuplot e_plot.scr",
 396                                  moldyn->vlsdir);
 397                         system(sc);
 398                 }
 399         }
 400         if(moldyn->mfd) close(moldyn->mfd);
 401         if(moldyn->pfd) {
 402                 close(moldyn->pfd);
 403                 if(moldyn->rfd)
 404                         dprintf(moldyn->rfd,report_pressure);
 405                         snprintf(sc,255,"cd %s && gnuplot pressure_plot.scr",
 406                                  moldyn->vlsdir);
 407                         system(sc);
 408         }
 409         if(moldyn->tfd) {
 410                 close(moldyn->tfd);
 411                 if(moldyn->rfd)
 412                         dprintf(moldyn->rfd,report_temperature);
 413                         snprintf(sc,255,"cd %s && gnuplot temperature_plot.scr",
 414                                  moldyn->vlsdir);
 415                         system(sc);
 416         }
 417         if(moldyn->rfd) {
 418                 dprintf(moldyn->rfd,report_end);
 419                 close(moldyn->rfd);
 420                 snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
 421                          moldyn->vlsdir);
 422                 system(sc);
 423                 snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1",
 424                          moldyn->vlsdir);
 425                 system(sc);
 426                 snprintf(sc,255,"cd %s && dvipdf report >/dev/null 2>&1",
 427                          moldyn->vlsdir);
 428                 system(sc);
 429         }
 430         if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
 431
 432         return 0;
 433 }
 434
 435 /*
 436  * creating lattice functions
 437  */
 438
 439 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 440                    u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) {
 441
 442         int new,count;
 443         int ret;
 444         t_3dvec orig;
 445         void *ptr;
 446         t_atom *atom;
 447
 448         new=a*b*c;
 449         count=moldyn->count;
 450
 451         /* how many atoms do we expect */
 452         if(type==CUBIC) new*=1;
 453         if(type==FCC) new*=4;
 454         if(type==DIAMOND) new*=8;
 455
 456         /* allocate space for atoms */
 457         ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom));
 458         if(!ptr) {
 459                 perror("[moldyn] realloc (create lattice)");
 460                 return -1;
 461         }
 462         moldyn->atom=ptr;
 463         atom=&(moldyn->atom[count]);
 464
 465         /* no atoms on the boundaries (only reason: it looks better!) */
 466         if(!origin) {
 467                 orig.x=0.5*lc;
 468                 orig.y=0.5*lc;
 469                 orig.z=0.5*lc;
 470         }
 471         else {
 472                 orig.x=origin->x;
 473                 orig.y=origin->y;
 474                 orig.z=origin->z;
 475         }
 476
 477         switch(type) {
 478                 case CUBIC:
 479                         set_nn_dist(moldyn,lc);
 480                         ret=cubic_init(a,b,c,lc,atom,&orig);
 481                         break;
 482                 case FCC:
 483                         if(!origin)
 484                                 v3_scale(&orig,&orig,0.5);
 485                         set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
 486                         ret=fcc_init(a,b,c,lc,atom,&orig);
 487                         break;
 488                 case DIAMOND:
 489                         if(!origin)
 490                                 v3_scale(&orig,&orig,0.25);
 491                         set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
 492                         ret=diamond_init(a,b,c,lc,atom,&orig);
 493                         break;
 494                 default:
 495                         printf("unknown lattice type (%02x)\n",type);
 496                         return -1;
 497         }
 498
 499         /* debug */
 500         if(ret!=new) {
 501                 printf("[moldyn] creating lattice failed\n");
 502                 printf("  amount of atoms\n");
 503                 printf("  - expected: %d\n",new);
 504                 printf("  - created: %d\n",ret);
 505                 return -1;
 506         }
 507
 508         moldyn->count+=new;
 509         printf("[moldyn] created lattice with %d atoms\n",new);
 510
 511         for(ret=0;ret<new;ret++) {
 512                 atom[ret].element=element;
 513                 atom[ret].mass=mass;
 514                 atom[ret].attr=attr;
 515                 atom[ret].brand=brand;
 516                 atom[ret].tag=count+ret;
 517                 check_per_bound(moldyn,&(atom[ret].r));
 518                 atom[ret].r_0=atom[ret].r;
 519         }
 520
 521         /* update total system mass */
 522         total_mass_calc(moldyn);
 523
 524         return ret;
 525 }
 526
 527 int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 528              t_3dvec *r,t_3dvec *v) {
 529
 530         t_atom *atom;
 531         void *ptr;
 532         int count;
 533
 534         atom=moldyn->atom;
 535         count=(moldyn->count)++;
 536
 537         ptr=realloc(atom,(count+1)*sizeof(t_atom));
 538         if(!ptr) {
 539                 perror("[moldyn] realloc (add atom)");
 540                 return -1;
 541         }
 542         moldyn->atom=ptr;
 543
 544         atom=moldyn->atom;
 545         atom[count].r=*r;
 546         atom[count].v=*v;
 547         atom[count].element=element;
 548         atom[count].mass=mass;
 549         atom[count].brand=brand;
 550         atom[count].tag=count;
 551         atom[count].attr=attr;
 552         check_per_bound(moldyn,&(atom[count].r));
 553         atom[count].r_0=atom[count].r;
 554
 555         /* update total system mass */
 556         total_mass_calc(moldyn);
 557
 558         return 0;
 559 }
 560
 561 int del_atom(t_moldyn *moldyn,int tag) {
 562
 563         t_atom *new,*old;
 564         int cnt;
 565
 566         old=moldyn->atom;
 567
 568         new=(t_atom *)malloc((moldyn->count-1)*sizeof(t_atom));
 569         if(!new) {
 570                 perror("[moldyn]malloc (del atom)");
 571                 return -1;
 572         }
 573
 574         for(cnt=0;cnt<tag;cnt++)
 575                 new[cnt]=old[cnt];
 576
 577         for(cnt=tag+1;cnt<moldyn->count;cnt++) {
 578                 new[cnt-1]=old[cnt];
 579                 new[cnt-1].tag=cnt-1;
 580         }
 581
 582         moldyn->count-=1;
 583         moldyn->atom=new;
 584
 585         free(old);
 586
 587         return 0;
 588 }
 589
 590 /* cubic init */
 591 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 592
 593         int count;
 594         t_3dvec r;
 595         int i,j,k;
 596         t_3dvec o;
 597
 598         count=0;
 599         if(origin)
 600                 v3_copy(&o,origin);
 601         else
 602                 v3_zero(&o);
 603
 604         r.x=o.x;
 605         for(i=0;i<a;i++) {
 606                 r.y=o.y;
 607                 for(j=0;j<b;j++) {
 608                         r.z=o.z;
 609                         for(k=0;k<c;k++) {
 610                                 v3_copy(&(atom[count].r),&r);
 611                                 count+=1;
 612                                 r.z+=lc;
 613                         }
 614                         r.y+=lc;
 615                 }
 616                 r.x+=lc;
 617         }
 618
 619         for(i=0;i<count;i++) {
 620                 atom[i].r.x-=(a*lc)/2.0;
 621                 atom[i].r.y-=(b*lc)/2.0;
 622                 atom[i].r.z-=(c*lc)/2.0;
 623         }
 624
 625         return count;
 626 }
 627
 628 /* fcc lattice init */
 629 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 630
 631         int count;
 632         int i,j,k,l;
 633         t_3dvec o,r,n;
 634         t_3dvec basis[3];
 635
 636         count=0;
 637         if(origin)
 638                 v3_copy(&o,origin);
 639         else
 640                 v3_zero(&o);
 641
 642         /* construct the basis */
 643         memset(basis,0,3*sizeof(t_3dvec));
 644         basis[0].x=0.5*lc;
 645         basis[0].y=0.5*lc;
 646         basis[1].x=0.5*lc;
 647         basis[1].z=0.5*lc;
 648         basis[2].y=0.5*lc;
 649         basis[2].z=0.5*lc;
 650
 651         /* fill up the room */
 652         r.x=o.x;
 653         for(i=0;i<a;i++) {
 654                 r.y=o.y;
 655                 for(j=0;j<b;j++) {
 656                         r.z=o.z;
 657                         for(k=0;k<c;k++) {
 658                                 /* first atom */
 659                                 v3_copy(&(atom[count].r),&r);
 660                                 count+=1;
 661                                 r.z+=lc;
 662                                 /* the three face centered atoms */
 663                                 for(l=0;l<3;l++) {
 664                                         v3_add(&n,&r,&basis[l]);
 665                                         v3_copy(&(atom[count].r),&n);
 666                                         count+=1;
 667                                 }
 668                         }
 669                         r.y+=lc;
 670                 }
 671                 r.x+=lc;
 672         }
 673
 674         /* coordinate transformation */
 675         for(i=0;i<count;i++) {
 676                 atom[i].r.x-=(a*lc)/2.0;
 677                 atom[i].r.y-=(b*lc)/2.0;
 678                 atom[i].r.z-=(c*lc)/2.0;
 679         }
 680
 681         return count;
 682 }
 683
 684 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 685
 686         int count;
 687         t_3dvec o;
 688
 689         count=fcc_init(a,b,c,lc,atom,origin);
 690
 691         o.x=0.25*lc;
 692         o.y=0.25*lc;
 693         o.z=0.25*lc;
 694
 695         if(origin) v3_add(&o,&o,origin);
 696
 697         count+=fcc_init(a,b,c,lc,&atom[count],&o);
 698
 699         return count;
 700 }
 701
 702 int destroy_atoms(t_moldyn *moldyn) {
 703
 704         if(moldyn->atom) free(moldyn->atom);
 705
 706         return 0;
 707 }
 708
 709 int thermal_init(t_moldyn *moldyn,u8 equi_init) {
 710
 711         /*
 712          * - gaussian distribution of velocities
 713          * - zero total momentum
 714          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 715          */
 716
 717         int i;
 718         double v,sigma;
 719         t_3dvec p_total,delta;
 720         t_atom *atom;
 721         t_random *random;
 722
 723         atom=moldyn->atom;
 724         random=&(moldyn->random);
 725
 726         printf("[moldyn] thermal init (equi init: %s)\n",equi_init?"yes":"no");
 727
 728         /* gaussian distribution of velocities */
 729         v3_zero(&p_total);
 730         for(i=0;i<moldyn->count;i++) {
 731                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
 732                 /* x direction */
 733                 v=sigma*rand_get_gauss(random);
 734                 atom[i].v.x=v;
 735                 p_total.x+=atom[i].mass*v;
 736                 /* y direction */
 737                 v=sigma*rand_get_gauss(random);
 738                 atom[i].v.y=v;
 739                 p_total.y+=atom[i].mass*v;
 740                 /* z direction */
 741                 v=sigma*rand_get_gauss(random);
 742                 atom[i].v.z=v;
 743                 p_total.z+=atom[i].mass*v;
 744         }
 745
 746         /* zero total momentum */
 747         v3_scale(&p_total,&p_total,1.0/moldyn->count);
 748         for(i=0;i<moldyn->count;i++) {
 749                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 750                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 751         }
 752
 753         /* velocity scaling */
 754         scale_velocity(moldyn,equi_init);
 755
 756         return 0;
 757 }
 758
 759 double total_mass_calc(t_moldyn *moldyn) {
 760
 761         int i;
 762
 763         moldyn->mass=0.0;
 764
 765         for(i=0;i<moldyn->count;i++)
 766                 moldyn->mass+=moldyn->atom[i].mass;
 767
 768         return moldyn->mass;
 769 }
 770
 771 double temperature_calc(t_moldyn *moldyn) {
 772
 773         /* assume up to date kinetic energy, which is 3/2 N k_B T */
 774
 775         moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
 776
 777         return moldyn->t;
 778 }
 779
 780 double get_temperature(t_moldyn *moldyn) {
 781
 782         return moldyn->t;
 783 }
 784
 785 int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 786
 787         int i;
 788         double e,scale;
 789         t_atom *atom;
 790         int count;
 791
 792         atom=moldyn->atom;
 793
 794         /*
 795          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 796          */
 797
 798         /* get kinetic energy / temperature & count involved atoms */
 799         e=0.0;
 800         count=0;
 801         for(i=0;i<moldyn->count;i++) {
 802                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) {
 803                         e+=atom[i].mass*v3_absolute_square(&(atom[i].v));
 804                         count+=1;
 805                 }
 806         }
 807         e*=0.5;
 808         if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN);
 809         else return 0;  /* no atoms involved in scaling! */
 810
 811         /* (temporary) hack for e,t = 0 */
 812         if(e==0.0) {
 813         moldyn->t=0.0;
 814                 if(moldyn->t_ref!=0.0) {
 815                         thermal_init(moldyn,equi_init);
 816                         return 0;
 817                 }
 818                 else
 819                         return 0; /* no scaling needed */
 820         }
 821
 822
 823         /* get scaling factor */
 824         scale=moldyn->t_ref/moldyn->t;
 825         if(equi_init&TRUE)
 826                 scale*=2.0;
 827         else
 828                 if(moldyn->pt_scale&T_SCALE_BERENDSEN)
 829                         scale=1.0+(scale-1.0)/moldyn->t_tc;
 830         scale=sqrt(scale);
 831
 832         /* velocity scaling */
 833         for(i=0;i<moldyn->count;i++) {
 834                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
 835                         v3_scale(&(atom[i].v),&(atom[i].v),scale);
 836         }
 837
 838         return 0;
 839 }
 840
 841 double ideal_gas_law_pressure(t_moldyn *moldyn) {
 842
 843         double p;
 844
 845         p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume;
 846
 847         return p;
 848 }
 849
 850 double virial_sum(t_moldyn *moldyn) {
 851
 852         int i;
 853         double v;
 854         t_virial *virial;
 855
 856         /* virial (sum over atom virials) */
 857         v=0.0;
 858         for(i=0;i<moldyn->count;i++) {
 859                 virial=&(moldyn->atom[i].virial);
 860                 v+=(virial->xx+virial->yy+virial->zz);
 861         }
 862         moldyn->virial=v;
 863
 864         /* global virial (absolute coordinates) */
 865         virial=&(moldyn->gvir);
 866         moldyn->gv=virial->xx+virial->yy+virial->zz;
 867
 868         return moldyn->virial;
 869 }
 870
 871 double pressure_calc(t_moldyn *moldyn) {
 872
 873         /*
 874          * PV = NkT + <W>
 875          * with W = 1/3 sum_i f_i r_i (- skipped!)
 876          * virial = sum_i f_i r_i
 877          *
 878          * => P = (2 Ekin + virial) / (3V)
 879          */
 880
 881         /* assume up to date virial & up to date kinetic energy */
 882
 883         /* pressure (atom virials) */
 884         moldyn->p=2.0*moldyn->ekin+moldyn->virial;
 885         moldyn->p/=(3.0*moldyn->volume);
 886
 887         /* pressure (absolute coordinates) */
 888         moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
 889         moldyn->gp/=(3.0*moldyn->volume);
 890
 891         return moldyn->p;
 892 }
 893
 894 int average_and_fluctuation_calc(t_moldyn *moldyn) {
 895
 896         if(moldyn->total_steps<moldyn->avg_skip)
 897                 return 0;
 898
 899         int denom=moldyn->total_steps+1-moldyn->avg_skip;
 900
 901         /* assume up to date energies, temperature, pressure etc */
 902
 903         /* kinetic energy */
 904         moldyn->k_sum+=moldyn->ekin;
 905         moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
 906         moldyn->k_avg=moldyn->k_sum/denom;
 907         moldyn->k2_avg=moldyn->k2_sum/denom;
 908         moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg);
 909
 910         /* potential energy */
 911         moldyn->v_sum+=moldyn->energy;
 912         moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
 913         moldyn->v_avg=moldyn->v_sum/denom;
 914         moldyn->v2_avg=moldyn->v2_sum/denom;
 915         moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg);
 916
 917         /* temperature */
 918         moldyn->t_sum+=moldyn->t;
 919         moldyn->t_avg=moldyn->t_sum/denom;
 920
 921         /* virial */
 922         moldyn->virial_sum+=moldyn->virial;
 923         moldyn->virial_avg=moldyn->virial_sum/denom;
 924         moldyn->gv_sum+=moldyn->gv;
 925         moldyn->gv_avg=moldyn->gv_sum/denom;
 926
 927         /* pressure */
 928         moldyn->p_sum+=moldyn->p;
 929         moldyn->p_avg=moldyn->p_sum/denom;
 930         moldyn->gp_sum+=moldyn->gp;
 931         moldyn->gp_avg=moldyn->gp_sum/denom;
 932
 933         return 0;
 934 }
 935
 936 int get_heat_capacity(t_moldyn *moldyn) {
 937
 938         double temp2,ighc;
 939
 940         /* averages needed for heat capacity calc */
 941         if(moldyn->total_steps<moldyn->avg_skip)
 942                 return 0;
 943
 944         /* (temperature average)^2 */
 945         temp2=moldyn->t_avg*moldyn->t_avg;
 946         printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
 947                moldyn->t_avg);
 948
 949         /* ideal gas contribution */
 950         ighc=3.0*moldyn->count*K_BOLTZMANN/2.0;
 951         printf("  ideal gas contribution: %f\n",
 952                ighc/moldyn->mass*KILOGRAM/JOULE);
 953
 954         /* specific heat for nvt ensemble */
 955         moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc;
 956         moldyn->c_v_nvt/=moldyn->mass;
 957
 958         /* specific heat for nve ensemble */
 959         moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2)));
 960         moldyn->c_v_nve/=moldyn->mass;
 961
 962         printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
 963         printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
 964 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)));
 965
 966         return 0;
 967 }
 968
 969 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 970
 971         t_3dvec dim,*tp;
 972         double u_up,u_down,dv;
 973         double scale,p;
 974         t_atom *store;
 975
 976         /*
 977          * dU = - p dV
 978          *
 979          * => p = - dU/dV
 980          *
 981          */
 982
 983         scale=0.00001;
 984         dv=8*scale*scale*scale*moldyn->volume;
 985
 986         store=malloc(moldyn->count*sizeof(t_atom));
 987         if(store==NULL) {
 988                 printf("[moldyn] allocating store mem failed\n");
 989                 return -1;
 990         }
 991
 992         /* save unscaled potential energy + atom/dim configuration */
 993         memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
 994         dim=moldyn->dim;
 995
 996         /* scale up dimension and atom positions */
 997         scale_dim(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
 998         scale_atoms(moldyn,SCALE_UP,scale,TRUE,TRUE,TRUE);
 999         link_cell_shutdown(moldyn);
1000         link_cell_init(moldyn,QUIET);
1001         potential_force_calc(moldyn);
1002         u_up=moldyn->energy;
1003
1004         /* restore atomic configuration + dim */
1005         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
1006         moldyn->dim=dim;
1007
1008         /* scale down dimension and atom positions */
1009         scale_dim(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
1010         scale_atoms(moldyn,SCALE_DOWN,scale,TRUE,TRUE,TRUE);
1011         link_cell_shutdown(moldyn);
1012         link_cell_init(moldyn,QUIET);
1013         potential_force_calc(moldyn);
1014         u_down=moldyn->energy;
1015
1016         /* calculate pressure */
1017         p=-(u_up-u_down)/dv;
1018 printf("-------> %.10f %.10f %f\n",u_up/EV/moldyn->count,u_down/EV/moldyn->count,p/BAR);
1019
1020         /* restore atomic configuration + dim */
1021         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
1022         moldyn->dim=dim;
1023
1024         /* restore energy */
1025         potential_force_calc(moldyn);
1026
1027         link_cell_shutdown(moldyn);
1028         link_cell_init(moldyn,QUIET);
1029
1030         return p;
1031 }
1032
1033 double get_pressure(t_moldyn *moldyn) {
1034
1035         return moldyn->p;
1036
1037 }
1038
1039 int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1040
1041         t_3dvec *dim;
1042
1043         dim=&(moldyn->dim);
1044
1045         if(dir==SCALE_UP)
1046                 scale=1.0+scale;
1047
1048         if(dir==SCALE_DOWN)
1049                 scale=1.0-scale;
1050
1051         if(x) dim->x*=scale;
1052         if(y) dim->y*=scale;
1053         if(z) dim->z*=scale;
1054
1055         return 0;
1056 }
1057
1058 int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
1059
1060         int i;
1061         t_3dvec *r;
1062
1063         if(dir==SCALE_UP)
1064                 scale=1.0+scale;
1065
1066         if(dir==SCALE_DOWN)
1067                 scale=1.0-scale;
1068
1069         for(i=0;i<moldyn->count;i++) {
1070                 r=&(moldyn->atom[i].r);
1071                 if(x) r->x*=scale;
1072                 if(y) r->y*=scale;
1073                 if(z) r->z*=scale;
1074         }
1075
1076         return 0;
1077 }
1078
1079 int scale_volume(t_moldyn *moldyn) {
1080
1081         t_3dvec *dim,*vdim;
1082         double scale;
1083         t_linkcell *lc;
1084
1085         vdim=&(moldyn->vis.dim);
1086         dim=&(moldyn->dim);
1087         lc=&(moldyn->lc);
1088
1089         /* scaling factor */
1090         if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
1091                 scale=1.0-(moldyn->p_ref-moldyn->p)/moldyn->p_tc;
1092                 scale=pow(scale,ONE_THIRD);
1093         }
1094         else {
1095                 scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
1096         }
1097 moldyn->debug=scale;
1098
1099         /* scale the atoms and dimensions */
1100         scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1101         scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
1102
1103         /* visualize dimensions */
1104         if(vdim->x!=0) {
1105                 vdim->x=dim->x;
1106                 vdim->y=dim->y;
1107                 vdim->z=dim->z;
1108         }
1109
1110         /* recalculate scaled volume */
1111         moldyn->volume=dim->x*dim->y*dim->z;
1112
1113         /* adjust/reinit linkcell */
1114         if(((int)(dim->x/moldyn->cutoff)!=lc->nx)||
1115            ((int)(dim->y/moldyn->cutoff)!=lc->ny)||
1116            ((int)(dim->z/moldyn->cutoff)!=lc->nx)) {
1117                 link_cell_shutdown(moldyn);
1118                 link_cell_init(moldyn,QUIET);
1119         } else {
1120                 lc->x*=scale;
1121                 lc->y*=scale;
1122                 lc->z*=scale;
1123         }
1124
1125         return 0;
1126
1127 }
1128
1129 double e_kin_calc(t_moldyn *moldyn) {
1130
1131         int i;
1132         t_atom *atom;
1133
1134         atom=moldyn->atom;
1135         moldyn->ekin=0.0;
1136
1137         for(i=0;i<moldyn->count;i++) {
1138                 atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
1139                 moldyn->ekin+=atom[i].ekin;
1140         }
1141
1142         return moldyn->ekin;
1143 }
1144
1145 double get_total_energy(t_moldyn *moldyn) {
1146
1147         return(moldyn->ekin+moldyn->energy);
1148 }
1149
1150 t_3dvec get_total_p(t_moldyn *moldyn) {
1151
1152         t_3dvec p,p_total;
1153         int i;
1154         t_atom *atom;
1155
1156         atom=moldyn->atom;
1157
1158         v3_zero(&p_total);
1159         for(i=0;i<moldyn->count;i++) {
1160                 v3_scale(&p,&(atom[i].v),atom[i].mass);
1161                 v3_add(&p_total,&p_total,&p);
1162         }
1163
1164         return p_total;
1165 }
1166
1167 double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
1168
1169         double tau;
1170
1171         /* nn_dist is the nearest neighbour distance */
1172
1173         tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
1174
1175         return tau;
1176 }
1177
1178 /*
1179  * numerical tricks
1180  */
1181
1182 /* linked list / cell method */
1183
1184 int link_cell_init(t_moldyn *moldyn,u8 vol) {
1185
1186         t_linkcell *lc;
1187         int i;
1188
1189         lc=&(moldyn->lc);
1190
1191         /* partitioning the md cell */
1192         lc->nx=moldyn->dim.x/moldyn->cutoff;
1193         lc->x=moldyn->dim.x/lc->nx;
1194         lc->ny=moldyn->dim.y/moldyn->cutoff;
1195         lc->y=moldyn->dim.y/lc->ny;
1196         lc->nz=moldyn->dim.z/moldyn->cutoff;
1197         lc->z=moldyn->dim.z/lc->nz;
1198
1199         lc->cells=lc->nx*lc->ny*lc->nz;
1200         lc->subcell=malloc(lc->cells*sizeof(t_list));
1201
1202         if(lc->cells<27)
1203                 printf("[moldyn] FATAL: less then 27 subcells!\n");
1204
1205         if(vol) {
1206                 printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
1207                 printf("  x: %d x %f A\n",lc->nx,lc->x);
1208                 printf("  y: %d x %f A\n",lc->ny,lc->y);
1209                 printf("  z: %d x %f A\n",lc->nz,lc->z);
1210         }
1211
1212         for(i=0;i<lc->cells;i++)
1213                 list_init_f(&(lc->subcell[i]));
1214
1215         link_cell_update(moldyn);
1216
1217         return 0;
1218 }
1219
1220 int link_cell_update(t_moldyn *moldyn) {
1221
1222         int count,i,j,k;
1223         int nx,ny;
1224         t_atom *atom;
1225         t_linkcell *lc;
1226         double x,y,z;
1227
1228         atom=moldyn->atom;
1229         lc=&(moldyn->lc);
1230
1231         nx=lc->nx;
1232         ny=lc->ny;
1233
1234         x=moldyn->dim.x/2;
1235         y=moldyn->dim.y/2;
1236         z=moldyn->dim.z/2;
1237
1238         for(i=0;i<lc->cells;i++)
1239                 list_destroy_f(&(lc->subcell[i]));
1240
1241         for(count=0;count<moldyn->count;count++) {
1242                 i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
1243                 j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
1244                 k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
1245                 list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
1246                                      &(atom[count]));
1247         }
1248
1249         return 0;
1250 }
1251
1252 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
1253
1254         t_linkcell *lc;
1255         int a;
1256         int count1,count2;
1257         int ci,cj,ck;
1258         int nx,ny,nz;
1259         int x,y,z;
1260         u8 bx,by,bz;
1261
1262         lc=&(moldyn->lc);
1263         nx=lc->nx;
1264         ny=lc->ny;
1265         nz=lc->nz;
1266         count1=1;
1267         count2=27;
1268         a=nx*ny;
1269
1270         cell[0]=lc->subcell[i+j*nx+k*a];
1271         for(ci=-1;ci<=1;ci++) {
1272                 bx=0;
1273                 x=i+ci;
1274                 if((x<0)||(x>=nx)) {
1275                         x=(x+nx)%nx;
1276                         bx=1;
1277                 }
1278                 for(cj=-1;cj<=1;cj++) {
1279                         by=0;
1280                         y=j+cj;
1281                         if((y<0)||(y>=ny)) {
1282                                 y=(y+ny)%ny;
1283                                 by=1;
1284                         }
1285                         for(ck=-1;ck<=1;ck++) {
1286                                 bz=0;
1287                                 z=k+ck;
1288                                 if((z<0)||(z>=nz)) {
1289                                         z=(z+nz)%nz;
1290                                         bz=1;
1291                                 }
1292                                 if(!(ci|cj|ck)) continue;
1293                                 if(bx|by|bz) {
1294                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
1295                                 }
1296                                 else {
1297                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
1298                                 }
1299                         }
1300                 }
1301         }
1302
1303         lc->dnlc=count1;
1304
1305         return count1;
1306 }
1307
1308 int link_cell_shutdown(t_moldyn *moldyn) {
1309
1310         int i;
1311         t_linkcell *lc;
1312
1313         lc=&(moldyn->lc);
1314
1315         for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
1316                 list_destroy_f(&(moldyn->lc.subcell[i]));
1317
1318         free(lc->subcell);
1319
1320         return 0;
1321 }
1322
1323 int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
1324
1325         int count;
1326         void *ptr;
1327         t_moldyn_schedule *schedule;
1328
1329         schedule=&(moldyn->schedule);
1330         count=++(schedule->total_sched);
1331
1332         ptr=realloc(schedule->runs,count*sizeof(int));
1333         if(!ptr) {
1334                 perror("[moldyn] realloc (runs)");
1335                 return -1;
1336         }
1337         schedule->runs=ptr;
1338         schedule->runs[count-1]=runs;
1339
1340         ptr=realloc(schedule->tau,count*sizeof(double));
1341         if(!ptr) {
1342                 perror("[moldyn] realloc (tau)");
1343                 return -1;
1344         }
1345         schedule->tau=ptr;
1346         schedule->tau[count-1]=tau;
1347
1348         printf("[moldyn] schedule added:\n");
1349         printf("  number: %d | runs: %d | tau: %f\n",count-1,runs,tau);
1350
1351
1352         return 0;
1353 }
1354
1355 int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
1356
1357         moldyn->schedule.hook=hook;
1358         moldyn->schedule.hook_params=hook_params;
1359
1360         return 0;
1361 }
1362
1363 /*
1364  *
1365  * 'integration of newtons equation' - algorithms
1366  *
1367  */
1368
1369 /* start the integration */
1370
1371 int moldyn_integrate(t_moldyn *moldyn) {
1372
1373         int i;
1374         unsigned int e,m,s,v,p,t;
1375         t_3dvec momentum;
1376         t_moldyn_schedule *sched;
1377         t_atom *atom;
1378         int fd;
1379         char dir[128];
1380         double ds;
1381         double energy_scale;
1382         struct timeval t1,t2;
1383         //double tp;
1384
1385         sched=&(moldyn->schedule);
1386         atom=moldyn->atom;
1387
1388         /* initialize linked cell method */
1389         link_cell_init(moldyn,VERBOSE);
1390
1391         /* logging & visualization */
1392         e=moldyn->ewrite;
1393         m=moldyn->mwrite;
1394         s=moldyn->swrite;
1395         v=moldyn->vwrite;
1396         p=moldyn->pwrite;
1397         t=moldyn->twrite;
1398
1399         /* sqaure of some variables */
1400         moldyn->tau_square=moldyn->tau*moldyn->tau;
1401         moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
1402
1403         /* get current time */
1404         gettimeofday(&t1,NULL);
1405
1406         /* calculate initial forces */
1407         potential_force_calc(moldyn);
1408 #ifdef DEBUG
1409 return 0;
1410 #endif
1411
1412         /* some stupid checks before we actually start calculating bullshit */
1413         if(moldyn->cutoff>0.5*moldyn->dim.x)
1414                 printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
1415         if(moldyn->cutoff>0.5*moldyn->dim.y)
1416                 printf("[moldyn] warning: cutoff > 0.5 x dim.y\n");
1417         if(moldyn->cutoff>0.5*moldyn->dim.z)
1418                 printf("[moldyn] warning: cutoff > 0.5 x dim.z\n");
1419         ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
1420         if(ds>0.05*moldyn->nnd)
1421                 printf("[moldyn] warning: forces too high / tau too small!\n");
1422
1423         /* zero absolute time */
1424         moldyn->time=0.0;
1425         moldyn->total_steps=0;
1426
1427         /* debugging, ignore */
1428         moldyn->debug=0;
1429
1430         /* tell the world */
1431         printf("[moldyn] integration start, go get a coffee ...\n");
1432
1433         /* executing the schedule */
1434         sched->count=0;
1435         while(sched->count<sched->total_sched) {
1436
1437                 /* setting amount of runs and finite time step size */
1438                 moldyn->tau=sched->tau[sched->count];
1439                 moldyn->tau_square=moldyn->tau*moldyn->tau;
1440                 moldyn->time_steps=sched->runs[sched->count];
1441
1442                 /* energy scaling factor (might change!) */
1443                 energy_scale=moldyn->count*EV;
1444
1445         /* integration according to schedule */
1446
1447         for(i=0;i<moldyn->time_steps;i++) {
1448
1449                 /* integration step */
1450                 moldyn->integrate(moldyn);
1451
1452                 /* calculate kinetic energy, temperature and pressure */
1453                 e_kin_calc(moldyn);
1454                 temperature_calc(moldyn);
1455                 virial_sum(moldyn);
1456                 pressure_calc(moldyn);
1457                 average_and_fluctuation_calc(moldyn);
1458
1459                 /* p/t scaling */
1460                 if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
1461                         scale_velocity(moldyn,FALSE);
1462                 if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
1463                         scale_volume(moldyn);
1464
1465                 /* check for log & visualization */
1466                 if(e) {
1467                         if(!(moldyn->total_steps%e))
1468                                 dprintf(moldyn->efd,
1469                                         "%f %f %f %f\n",
1470                                         moldyn->time,moldyn->ekin/energy_scale,
1471                                         moldyn->energy/energy_scale,
1472                                         get_total_energy(moldyn)/energy_scale);
1473                 }
1474                 if(m) {
1475                         if(!(moldyn->total_steps%m)) {
1476                                 momentum=get_total_p(moldyn);
1477                                 dprintf(moldyn->mfd,
1478                                         "%f %f %f %f %f\n",moldyn->time,
1479                                         momentum.x,momentum.y,momentum.z,
1480                                         v3_norm(&momentum));
1481                         }
1482                 }
1483                 if(p) {
1484                         if(!(moldyn->total_steps%p)) {
1485                                 dprintf(moldyn->pfd,
1486                                         "%f %f %f %f %f\n",moldyn->time,
1487                                          moldyn->p/BAR,moldyn->p_avg/BAR,
1488                                          moldyn->gp/BAR,moldyn->gp_avg/BAR);
1489                         }
1490                 }
1491                 if(t) {
1492                         if(!(moldyn->total_steps%t)) {
1493                                 dprintf(moldyn->tfd,
1494                                         "%f %f %f\n",
1495                                         moldyn->time,moldyn->t,moldyn->t_avg);
1496                         }
1497                 }
1498                 if(s) {
1499                         if(!(moldyn->total_steps%s)) {
1500                                 snprintf(dir,128,"%s/s-%07.f.save",
1501                                          moldyn->vlsdir,moldyn->time);
1502                                 fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
1503                                         S_IRUSR|S_IWUSR);
1504                                 if(fd<0) perror("[moldyn] save fd open");
1505                                 else {
1506                                         write(fd,moldyn,sizeof(t_moldyn));
1507                                         write(fd,moldyn->atom,
1508                                               moldyn->count*sizeof(t_atom));
1509                                 }
1510                                 close(fd);
1511                         }
1512                 }
1513                 if(v) {
1514                         if(!(moldyn->total_steps%v)) {
1515                                 visual_atoms(&(moldyn->vis),moldyn->time,
1516                                              moldyn->atom,moldyn->count);
1517                         }
1518                 }
1519
1520                 /* display progress */
1521                 if(!(moldyn->total_steps%10)) {
1522                         /* get current time */
1523                         gettimeofday(&t2,NULL);
1524
1525         printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
1526                sched->count,i,
1527                moldyn->t,moldyn->t_avg,
1528                moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
1529                moldyn->volume,
1530                (int)(t2.tv_sec-t1.tv_sec));
1531         fflush(stdout);
1532
1533                         /* copy over time */
1534                         t1=t2;
1535                 }
1536
1537                 /* increase absolute time */
1538                 moldyn->time+=moldyn->tau;
1539                 moldyn->total_steps+=1;
1540
1541         }
1542
1543                 /* check for hooks */
1544                 if(sched->hook) {
1545                         printf("\n ## schedule hook %d/%d start ##\n",
1546                                sched->count+1,sched->total_sched-1);
1547                         sched->hook(moldyn,sched->hook_params);
1548                         printf(" ## schedule hook end ##\n");
1549                 }
1550
1551                 /* increase the schedule counter */
1552                 sched->count+=1;
1553
1554         }
1555
1556         return 0;
1557 }
1558
1559 /* velocity verlet */
1560
1561 int velocity_verlet(t_moldyn *moldyn) {
1562
1563         int i,count;
1564         double tau,tau_square,h;
1565         t_3dvec delta;
1566         t_atom *atom;
1567
1568         atom=moldyn->atom;
1569         count=moldyn->count;
1570         tau=moldyn->tau;
1571         tau_square=moldyn->tau_square;
1572
1573         for(i=0;i<count;i++) {
1574                 /* new positions */
1575                 h=0.5/atom[i].mass;
1576                 v3_scale(&delta,&(atom[i].v),tau);
1577                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1578                 v3_scale(&delta,&(atom[i].f),h*tau_square);
1579                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1580                 check_per_bound(moldyn,&(atom[i].r));
1581
1582                 /* velocities [actually v(t+tau/2)] */
1583                 v3_scale(&delta,&(atom[i].f),h*tau);
1584                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1585         }
1586
1587         /* neighbour list update */
1588         link_cell_update(moldyn);
1589
1590         /* forces depending on chosen potential */
1591         potential_force_calc(moldyn);
1592
1593         for(i=0;i<count;i++) {
1594                 /* again velocities [actually v(t+tau)] */
1595                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
1596                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1597         }
1598
1599         return 0;
1600 }
1601
1602
1603 /*
1604  *
1605  * potentials & corresponding forces & virial routine
1606  *
1607  */
1608
1609 /* generic potential and force calculation */
1610
1611 int potential_force_calc(t_moldyn *moldyn) {
1612
1613         int i,j,k,count;
1614         t_atom *itom,*jtom,*ktom;
1615         t_virial *virial;
1616         t_linkcell *lc;
1617         t_list neighbour_i[27];
1618         t_list neighbour_i2[27];
1619         t_list *this,*that;
1620         u8 bc_ij,bc_ik;
1621         int dnlc;
1622
1623         count=moldyn->count;
1624         itom=moldyn->atom;
1625         lc=&(moldyn->lc);
1626
1627         /* reset energy */
1628         moldyn->energy=0.0;
1629
1630         /* reset global virial */
1631         memset(&(moldyn->gvir),0,sizeof(t_virial));
1632
1633         /* reset force, site energy and virial of every atom */
1634         for(i=0;i<count;i++) {
1635
1636                 /* reset force */
1637                 v3_zero(&(itom[i].f));
1638
1639                 /* reset virial */
1640                 virial=(&(itom[i].virial));
1641                 virial->xx=0.0;
1642                 virial->yy=0.0;
1643                 virial->zz=0.0;
1644                 virial->xy=0.0;
1645                 virial->xz=0.0;
1646                 virial->yz=0.0;
1647
1648                 /* reset site energy */
1649                 itom[i].e=0.0;
1650
1651         }
1652
1653         /* get energy, force and virial of every atom */
1654
1655         /* first (and only) loop over atoms i */
1656         for(i=0;i<count;i++) {
1657
1658                 /* single particle potential/force */
1659                 if(itom[i].attr&ATOM_ATTR_1BP)
1660                         if(moldyn->func1b)
1661                                 moldyn->func1b(moldyn,&(itom[i]));
1662
1663                 if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
1664                         continue;
1665
1666                 /* 2 body pair potential/force */
1667
1668                 link_cell_neighbour_index(moldyn,
1669                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
1670                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
1671                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
1672                                           neighbour_i);
1673
1674                 dnlc=lc->dnlc;
1675
1676                 /* first loop over atoms j */
1677                 if(moldyn->func2b) {
1678                         for(j=0;j<27;j++) {
1679
1680                                 this=&(neighbour_i[j]);
1681                                 list_reset_f(this);
1682
1683                                 if(this->start==NULL)
1684                                         continue;
1685
1686                                 bc_ij=(j<dnlc)?0:1;
1687
1688                                 do {
1689                                         jtom=this->current->data;
1690
1691                                         if(jtom==&(itom[i]))
1692                                                 continue;
1693
1694                                         if((jtom->attr&ATOM_ATTR_2BP)&
1695                                            (itom[i].attr&ATOM_ATTR_2BP)) {
1696                                                 moldyn->func2b(moldyn,
1697                                                                &(itom[i]),
1698                                                                jtom,
1699                                                                bc_ij);
1700                                         }
1701                                 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1702
1703                         }
1704                 }
1705
1706                 /* 3 body potential/force */
1707
1708                 if(!(itom[i].attr&ATOM_ATTR_3BP))
1709                         continue;
1710
1711                 /* copy the neighbour lists */
1712                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
1713
1714                 /* second loop over atoms j */
1715                 for(j=0;j<27;j++) {
1716
1717                         this=&(neighbour_i[j]);
1718                         list_reset_f(this);
1719
1720                         if(this->start==NULL)
1721                                 continue;
1722
1723                         bc_ij=(j<dnlc)?0:1;
1724
1725                         do {
1726                                 jtom=this->current->data;
1727
1728                                 if(jtom==&(itom[i]))
1729                                         continue;
1730
1731                                 if(!(jtom->attr&ATOM_ATTR_3BP))
1732                                         continue;
1733
1734                                 /* reset 3bp run */
1735                                 moldyn->run3bp=1;
1736
1737                                 if(moldyn->func3b_j1)
1738                                         moldyn->func3b_j1(moldyn,
1739                                                           &(itom[i]),
1740                                                           jtom,
1741                                                           bc_ij);
1742
1743                                 /* in first j loop, 3bp run can be skipped */
1744                                 if(!(moldyn->run3bp))
1745                                         continue;
1746
1747                                 /* first loop over atoms k */
1748                                 if(moldyn->func3b_k1) {
1749
1750                                 for(k=0;k<27;k++) {
1751
1752                                         that=&(neighbour_i2[k]);
1753                                         list_reset_f(that);
1754
1755                                         if(that->start==NULL)
1756                                                 continue;
1757
1758                                         bc_ik=(k<dnlc)?0:1;
1759
1760                                         do {
1761
1762                                                 ktom=that->current->data;
1763
1764                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
1765                                                         continue;
1766
1767                                                 if(ktom==jtom)
1768                                                         continue;
1769
1770                                                 if(ktom==&(itom[i]))
1771                                                         continue;
1772
1773                                                 moldyn->func3b_k1(moldyn,
1774                                                                   &(itom[i]),
1775                                                                   jtom,
1776                                                                   ktom,
1777                                                                   bc_ik|bc_ij);
1778
1779                                         } while(list_next_f(that)!=\
1780                                                 L_NO_NEXT_ELEMENT);
1781
1782                                 }
1783
1784                                 }
1785
1786                                 if(moldyn->func3b_j2)
1787                                         moldyn->func3b_j2(moldyn,
1788                                                           &(itom[i]),
1789                                                           jtom,
1790                                                           bc_ij);
1791
1792                                 /* second loop over atoms k */
1793                                 if(moldyn->func3b_k2) {
1794
1795                                 for(k=0;k<27;k++) {
1796
1797                                         that=&(neighbour_i2[k]);
1798                                         list_reset_f(that);
1799
1800                                         if(that->start==NULL)
1801                                                 continue;
1802
1803                                         bc_ik=(k<dnlc)?0:1;
1804
1805                                         do {
1806
1807                                                 ktom=that->current->data;
1808
1809                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
1810                                                         continue;
1811
1812                                                 if(ktom==jtom)
1813                                                         continue;
1814
1815                                                 if(ktom==&(itom[i]))
1816                                                         continue;
1817
1818                                                 moldyn->func3b_k2(moldyn,
1819                                                                   &(itom[i]),
1820                                                                   jtom,
1821                                                                   ktom,
1822                                                                   bc_ik|bc_ij);
1823
1824                                         } while(list_next_f(that)!=\
1825                                                 L_NO_NEXT_ELEMENT);
1826
1827                                 }
1828
1829                                 }
1830
1831                                 /* 2bp post function */
1832                                 if(moldyn->func3b_j3) {
1833                                         moldyn->func3b_j3(moldyn,
1834                                                           &(itom[i]),
1835                                                           jtom,bc_ij);
1836                                 }
1837
1838                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1839
1840                 }
1841
1842 #ifdef DEBUG
1843         //printf("\n\n");
1844 #endif
1845 #ifdef VDEBUG
1846         printf("\n\n");
1847 #endif
1848
1849         }
1850
1851 #ifdef DEBUG
1852         printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
1853 #endif
1854
1855         /* calculate global virial */
1856         for(i=0;i<count;i++) {
1857                 moldyn->gvir.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x;
1858                 moldyn->gvir.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y;
1859                 moldyn->gvir.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z;
1860                 moldyn->gvir.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x;
1861                 moldyn->gvir.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x;
1862                 moldyn->gvir.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y;
1863         }
1864
1865         return 0;
1866 }
1867
1868 /*
1869  * virial calculation
1870  */
1871
1872 //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
1873 int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
1874
1875         a->virial.xx+=f->x*d->x;
1876         a->virial.yy+=f->y*d->y;
1877         a->virial.zz+=f->z*d->z;
1878         a->virial.xy+=f->x*d->y;
1879         a->virial.xz+=f->x*d->z;
1880         a->virial.yz+=f->y*d->z;
1881
1882         return 0;
1883 }
1884
1885 /*
1886  * periodic boundary checking
1887  */
1888
1889 //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
1890 int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
1891
1892         double x,y,z;
1893         t_3dvec *dim;
1894
1895         dim=&(moldyn->dim);
1896
1897         x=dim->x/2;
1898         y=dim->y/2;
1899         z=dim->z/2;
1900
1901         if(moldyn->status&MOLDYN_STAT_PBX) {
1902                 if(a->x>=x) a->x-=dim->x;
1903                 else if(-a->x>x) a->x+=dim->x;
1904         }
1905         if(moldyn->status&MOLDYN_STAT_PBY) {
1906                 if(a->y>=y) a->y-=dim->y;
1907                 else if(-a->y>y) a->y+=dim->y;
1908         }
1909         if(moldyn->status&MOLDYN_STAT_PBZ) {
1910                 if(a->z>=z) a->z-=dim->z;
1911                 else if(-a->z>z) a->z+=dim->z;
1912         }
1913
1914         return 0;
1915 }
1916
1917 /*
1918  * debugging / critical check functions
1919  */
1920
1921 int moldyn_bc_check(t_moldyn *moldyn) {
1922
1923         t_atom *atom;
1924         t_3dvec *dim;
1925         int i;
1926         double x;
1927         u8 byte;
1928         int j,k;
1929
1930         atom=moldyn->atom;
1931         dim=&(moldyn->dim);
1932         x=dim->x/2;
1933
1934         for(i=0;i<moldyn->count;i++) {
1935                 if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {
1936                         printf("FATAL: atom %d: x: %.20f (%.20f)\n",
1937                                i,atom[i].r.x,dim->x/2);
1938                         printf("diagnostic:\n");
1939                         printf("-----------\natom.r.x:\n");
1940                         for(j=0;j<8;j++) {
1941                                 memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
1942                                 for(k=0;k<8;k++)
1943                                         printf("%d%c",
1944                                         ((byte)&(1<<k))?1:0,
1945                                         (k==7)?'\n':'|');
1946                         }
1947                         printf("---------------\nx=dim.x/2:\n");
1948                         for(j=0;j<8;j++) {
1949                                 memcpy(&byte,(u8 *)(&x)+j,1);
1950                                 for(k=0;k<8;k++)
1951                                         printf("%d%c",
1952                                         ((byte)&(1<<k))?1:0,
1953                                         (k==7)?'\n':'|');
1954                         }
1955                         if(atom[i].r.x==x) printf("the same!\n");
1956                         else printf("different!\n");
1957                 }
1958                 if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
1959                         printf("FATAL: atom %d: y: %.20f (%.20f)\n",
1960                                i,atom[i].r.y,dim->y/2);
1961                 if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
1962                         printf("FATAL: atom %d: z: %.20f (%.20f)\n",
1963                                i,atom[i].r.z,dim->z/2);
1964         }
1965
1966         return 0;
1967 }
1968
1969 /*
1970  * restore function
1971  */
1972
1973 int moldyn_load(t_moldyn *moldyn) {
1974
1975         // later ...
1976
1977         return 0;
1978 }
1979
1980 /*
1981  * post processing functions
1982  */
1983
1984 int get_line(int fd,char *line,int max) {
1985
1986         int count,ret;
1987
1988         count=0;
1989
1990         while(1) {
1991                 if(count==max) return count;
1992                 ret=read(fd,line+count,1);
1993                 if(ret<=0) return ret;
1994                 if(line[count]=='\n') {
1995                         line[count]='\0';
1996                         return count+1;
1997                 }
1998                 count+=1;
1999         }
2000 }
2001
2002 int analyze_bonds(t_moldyn *moldyn) {
2003
2004
2005
2006
2007         return 0;
2008 }
2009