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 <math.h>
  17
  18 #include "moldyn.h"
  19 #include "report/report.h"
  20
  21 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
  22
  23         printf("[moldyn] init\n");
  24
  25         memset(moldyn,0,sizeof(t_moldyn));
  26
  27         rand_init(&(moldyn->random),NULL,1);
  28         moldyn->random.status|=RAND_STAT_VERBOSE;
  29
  30         return 0;
  31 }
  32
  33 int moldyn_shutdown(t_moldyn *moldyn) {
  34
  35         printf("[moldyn] shutdown\n");
  36
  37         moldyn_log_shutdown(moldyn);
  38         link_cell_shutdown(moldyn);
  39         rand_close(&(moldyn->random));
  40         free(moldyn->atom);
  41
  42         return 0;
  43 }
  44
  45 int set_int_alg(t_moldyn *moldyn,u8 algo) {
  46
  47         printf("[moldyn] integration algorithm: ");
  48
  49         switch(algo) {
  50                 case MOLDYN_INTEGRATE_VERLET:
  51                         moldyn->integrate=velocity_verlet;
  52                         printf("velocity verlet\n");
  53                         break;
  54                 default:
  55                         printf("unknown integration algorithm: %02x\n",algo);
  56                         printf("unknown\n");
  57                         return -1;
  58         }
  59
  60         return 0;
  61 }
  62
  63 int set_cutoff(t_moldyn *moldyn,double cutoff) {
  64
  65         moldyn->cutoff=cutoff;
  66
  67         printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff);
  68
  69         return 0;
  70 }
  71
  72 int set_temperature(t_moldyn *moldyn,double t_ref) {
  73
  74         moldyn->t_ref=t_ref;
  75
  76         printf("[moldyn] temperature [K]: %f\n",moldyn->t_ref);
  77
  78         return 0;
  79 }
  80
  81 int set_pressure(t_moldyn *moldyn,double p_ref) {
  82
  83         moldyn->p_ref=p_ref;
  84
  85         printf("[moldyn] pressure [atm]: %f\n",moldyn->p_ref/ATM);
  86
  87         return 0;
  88 }
  89
  90 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
  91
  92         moldyn->pt_scale=(ptype|ttype);
  93         moldyn->t_tc=ttc;
  94         moldyn->p_tc=ptc;
  95
  96         printf("[moldyn] p/t scaling:\n");
  97
  98         printf("  p: %s",ptype?"yes":"no ");
  99         if(ptype)
 100                 printf(" | type: %02x | factor: %f",ptype,ptc);
 101         printf("\n");
 102
 103         printf("  t: %s",ttype?"yes":"no ");
 104         if(ttype)
 105                 printf(" | type: %02x | factor: %f",ttype,ttc);
 106         printf("\n");
 107
 108         return 0;
 109 }
 110
 111 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
 112
 113         moldyn->dim.x=x;
 114         moldyn->dim.y=y;
 115         moldyn->dim.z=z;
 116
 117         moldyn->volume=x*y*z;
 118
 119         if(visualize) {
 120                 moldyn->vis.dim.x=x;
 121                 moldyn->vis.dim.y=y;
 122                 moldyn->vis.dim.z=z;
 123         }
 124
 125         moldyn->dv=0.000001*moldyn->volume;
 126
 127         printf("[moldyn] dimensions in A and A^3 respectively:\n");
 128         printf("  x: %f\n",moldyn->dim.x);
 129         printf("  y: %f\n",moldyn->dim.y);
 130         printf("  z: %f\n",moldyn->dim.z);
 131         printf("  volume: %f\n",moldyn->volume);
 132         printf("  visualize simulation box: %s\n",visualize?"yes":"no");
 133         printf("  delta volume (pressure calc): %f\n",moldyn->dv);
 134
 135         return 0;
 136 }
 137
 138 int set_nn_dist(t_moldyn *moldyn,double dist) {
 139
 140         moldyn->nnd=dist;
 141
 142         return 0;
 143 }
 144
 145 int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
 146
 147         printf("[moldyn] periodic boundary conditions:\n");
 148
 149         if(x)
 150                 moldyn->status|=MOLDYN_STAT_PBX;
 151
 152         if(y)
 153                 moldyn->status|=MOLDYN_STAT_PBY;
 154
 155         if(z)
 156                 moldyn->status|=MOLDYN_STAT_PBZ;
 157
 158         printf("  x: %s\n",x?"yes":"no");
 159         printf("  y: %s\n",y?"yes":"no");
 160         printf("  z: %s\n",z?"yes":"no");
 161
 162         return 0;
 163 }
 164
 165 int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params) {
 166
 167         moldyn->func1b=func;
 168
 169         return 0;
 170 }
 171
 172 int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) {
 173
 174         moldyn->func2b=func;
 175
 176         return 0;
 177 }
 178
 179 int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) {
 180
 181         moldyn->func2b_post=func;
 182
 183         return 0;
 184 }
 185
 186 int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) {
 187
 188         moldyn->func3b=func;
 189
 190         return 0;
 191 }
 192
 193 int set_potential_params(t_moldyn *moldyn,void *params) {
 194
 195         moldyn->pot_params=params;
 196
 197         return 0;
 198 }
 199
 200 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
 201
 202         strncpy(moldyn->vlsdir,dir,127);
 203
 204         return 0;
 205 }
 206
 207 int moldyn_set_report(t_moldyn *moldyn,char *author,char *title) {
 208
 209         strncpy(moldyn->rauthor,author,63);
 210         strncpy(moldyn->rtitle,title,63);
 211
 212         return 0;
 213 }
 214
 215 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 216
 217         char filename[128];
 218         int ret;
 219
 220         printf("[moldyn] set log: ");
 221
 222         switch(type) {
 223                 case LOG_TOTAL_ENERGY:
 224                         moldyn->ewrite=timer;
 225                         snprintf(filename,127,"%s/energy",moldyn->vlsdir);
 226                         moldyn->efd=open(filename,
 227                                          O_WRONLY|O_CREAT|O_EXCL,
 228                                          S_IRUSR|S_IWUSR);
 229                         if(moldyn->efd<0) {
 230                                 perror("[moldyn] energy log fd open");
 231                                 return moldyn->efd;
 232                         }
 233                         dprintf(moldyn->efd,"# total energy log file\n");
 234                         printf("total energy (%d)\n",timer);
 235                         break;
 236                 case LOG_TOTAL_MOMENTUM:
 237                         moldyn->mwrite=timer;
 238                         snprintf(filename,127,"%s/momentum",moldyn->vlsdir);
 239                         moldyn->mfd=open(filename,
 240                                          O_WRONLY|O_CREAT|O_EXCL,
 241                                          S_IRUSR|S_IWUSR);
 242                         if(moldyn->mfd<0) {
 243                                 perror("[moldyn] momentum log fd open");
 244                                 return moldyn->mfd;
 245                         }
 246                         dprintf(moldyn->efd,"# total momentum log file\n");
 247                         printf("total momentum (%d)\n",timer);
 248                         break;
 249                 case SAVE_STEP:
 250                         moldyn->swrite=timer;
 251                         printf("save file (%d)\n",timer);
 252                         break;
 253                 case VISUAL_STEP:
 254                         moldyn->vwrite=timer;
 255                         ret=visual_init(&(moldyn->vis),moldyn->vlsdir);
 256                         if(ret<0) {
 257                                 printf("[moldyn] visual init failure\n");
 258                                 return ret;
 259                         }
 260                         printf("visual file (%d)\n",timer);
 261                         break;
 262                 case CREATE_REPORT:
 263                         snprintf(filename,127,"%s/report.tex",moldyn->vlsdir);
 264                         moldyn->rfd=open(filename,
 265                                          O_WRONLY|O_CREAT|O_EXCL,
 266                                          S_IRUSR|S_IWUSR);
 267                         if(moldyn->rfd<0) {
 268                                 perror("[moldyn] report fd open");
 269                                 return moldyn->rfd;
 270                         }
 271                         snprintf(filename,127,"%s/plot.scr",moldyn->vlsdir);
 272                         moldyn->pfd=open(filename,
 273                                          O_WRONLY|O_CREAT|O_EXCL,
 274                                          S_IRUSR|S_IWUSR);
 275                         if(moldyn->pfd<0) {
 276                                 perror("[moldyn] plot fd open");
 277                                 return moldyn->pfd;
 278                         }
 279                         dprintf(moldyn->rfd,report_start,
 280                                 moldyn->rauthor,moldyn->rtitle);
 281                         dprintf(moldyn->pfd,plot_script);
 282                         close(moldyn->pfd);
 283                         break;
 284                 default:
 285                         printf("unknown log type: %02x\n",type);
 286                         return -1;
 287         }
 288
 289         return 0;
 290 }
 291
 292 int moldyn_log_shutdown(t_moldyn *moldyn) {
 293
 294         char sc[256];
 295
 296         printf("[moldyn] log shutdown\n");
 297         if(moldyn->efd) close(moldyn->efd);
 298         if(moldyn->mfd) close(moldyn->mfd);
 299         if(moldyn->rfd) {
 300                 dprintf(moldyn->rfd,report_end);
 301                 close(moldyn->rfd);
 302                 snprintf(sc,255,"cd %s && gnuplot plot.scr",moldyn->vlsdir);
 303                 system(sc);
 304                 snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir);
 305                 system(sc);
 306                 snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir);
 307                 system(sc);
 308                 snprintf(sc,255,"cd %s && dvipdf report",moldyn->vlsdir);
 309                 system(sc);
 310         }
 311         if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
 312
 313         return 0;
 314 }
 315
 316 /*
 317  * creating lattice functions
 318  */
 319
 320 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 321                    u8 attr,u8 brand,int a,int b,int c) {
 322
 323         int new,count;
 324         int ret;
 325         t_3dvec origin;
 326         void *ptr;
 327         t_atom *atom;
 328
 329         new=a*b*c;
 330         count=moldyn->count;
 331
 332         /* how many atoms do we expect */
 333         if(type==CUBIC) new*=1;
 334         if(type==FCC) new*=4;
 335         if(type==DIAMOND) new*=8;
 336
 337         /* allocate space for atoms */
 338         ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom));
 339         if(!ptr) {
 340                 perror("[moldyn] realloc (create lattice)");
 341                 return -1;
 342         }
 343         moldyn->atom=ptr;
 344         atom=&(moldyn->atom[count]);
 345
 346         /* no atoms on the boundaries (only reason: it looks better!) */
 347         origin.x=0.5*lc;
 348         origin.y=0.5*lc;
 349         origin.z=0.5*lc;
 350
 351         switch(type) {
 352                 case CUBIC:
 353                         set_nn_dist(moldyn,lc);
 354                         ret=cubic_init(a,b,c,lc,atom,&origin);
 355                         break;
 356                 case FCC:
 357                         v3_scale(&origin,&origin,0.5);
 358                         set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
 359                         ret=fcc_init(a,b,c,lc,atom,&origin);
 360                         break;
 361                 case DIAMOND:
 362                         v3_scale(&origin,&origin,0.25);
 363                         set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
 364                         ret=diamond_init(a,b,c,lc,atom,&origin);
 365                         break;
 366                 default:
 367                         printf("unknown lattice type (%02x)\n",type);
 368                         return -1;
 369         }
 370
 371         /* debug */
 372         if(ret!=new) {
 373                 printf("[moldyn] creating lattice failed\n");
 374                 printf("  amount of atoms\n");
 375                 printf("  - expected: %d\n",new);
 376                 printf("  - created: %d\n",ret);
 377                 return -1;
 378         }
 379
 380         moldyn->count+=new;
 381         printf("[moldyn] created lattice with %d atoms\n",new);
 382
 383         for(ret=0;ret<new;ret++) {
 384                 atom[ret].element=element;
 385                 atom[ret].mass=mass;
 386                 atom[ret].attr=attr;
 387                 atom[ret].brand=brand;
 388                 atom[ret].tag=count+ret;
 389                 check_per_bound(moldyn,&(atom[ret].r));
 390         }
 391
 392         return ret;
 393 }
 394
 395 /* cubic init */
 396 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 397
 398         int count;
 399         t_3dvec r;
 400         int i,j,k;
 401         t_3dvec o;
 402
 403         count=0;
 404         if(origin)
 405                 v3_copy(&o,origin);
 406         else
 407                 v3_zero(&o);
 408
 409         r.x=o.x;
 410         for(i=0;i<a;i++) {
 411                 r.y=o.y;
 412                 for(j=0;j<b;j++) {
 413                         r.z=o.z;
 414                         for(k=0;k<c;k++) {
 415                                 v3_copy(&(atom[count].r),&r);
 416                                 count+=1;
 417                                 r.z+=lc;
 418                         }
 419                         r.y+=lc;
 420                 }
 421                 r.x+=lc;
 422         }
 423
 424         for(i=0;i<count;i++) {
 425                 atom[i].r.x-=(a*lc)/2.0;
 426                 atom[i].r.y-=(b*lc)/2.0;
 427                 atom[i].r.z-=(c*lc)/2.0;
 428         }
 429
 430         return count;
 431 }
 432
 433 /* fcc lattice init */
 434 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 435
 436         int count;
 437         int i,j,k,l;
 438         t_3dvec o,r,n;
 439         t_3dvec basis[3];
 440
 441         count=0;
 442         if(origin)
 443                 v3_copy(&o,origin);
 444         else
 445                 v3_zero(&o);
 446
 447         /* construct the basis */
 448         memset(basis,0,3*sizeof(t_3dvec));
 449         basis[0].x=0.5*lc;
 450         basis[0].y=0.5*lc;
 451         basis[1].x=0.5*lc;
 452         basis[1].z=0.5*lc;
 453         basis[2].y=0.5*lc;
 454         basis[2].z=0.5*lc;
 455
 456         /* fill up the room */
 457         r.x=o.x;
 458         for(i=0;i<a;i++) {
 459                 r.y=o.y;
 460                 for(j=0;j<b;j++) {
 461                         r.z=o.z;
 462                         for(k=0;k<c;k++) {
 463                                 /* first atom */
 464                                 v3_copy(&(atom[count].r),&r);
 465                                 count+=1;
 466                                 r.z+=lc;
 467                                 /* the three face centered atoms */
 468                                 for(l=0;l<3;l++) {
 469                                         v3_add(&n,&r,&basis[l]);
 470                                         v3_copy(&(atom[count].r),&n);
 471                                         count+=1;
 472                                 }
 473                         }
 474                         r.y+=lc;
 475                 }
 476                 r.x+=lc;
 477         }
 478
 479         /* coordinate transformation */
 480         for(i=0;i<count;i++) {
 481                 atom[i].r.x-=(a*lc)/2.0;
 482                 atom[i].r.y-=(b*lc)/2.0;
 483                 atom[i].r.z-=(c*lc)/2.0;
 484         }
 485
 486         return count;
 487 }
 488
 489 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 490
 491         int count;
 492         t_3dvec o;
 493
 494         count=fcc_init(a,b,c,lc,atom,origin);
 495
 496         o.x=0.25*lc;
 497         o.y=0.25*lc;
 498         o.z=0.25*lc;
 499
 500         if(origin) v3_add(&o,&o,origin);
 501
 502         count+=fcc_init(a,b,c,lc,&atom[count],&o);
 503
 504         return count;
 505 }
 506
 507 int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 508              t_3dvec *r,t_3dvec *v) {
 509
 510         t_atom *atom;
 511         void *ptr;
 512         int count;
 513
 514         atom=moldyn->atom;
 515         count=(moldyn->count)++;
 516
 517         ptr=realloc(atom,(count+1)*sizeof(t_atom));
 518         if(!ptr) {
 519                 perror("[moldyn] realloc (add atom)");
 520                 return -1;
 521         }
 522         moldyn->atom=ptr;
 523
 524         atom=moldyn->atom;
 525         atom[count].r=*r;
 526         atom[count].v=*v;
 527         atom[count].element=element;
 528         atom[count].mass=mass;
 529         atom[count].brand=brand;
 530         atom[count].tag=count;
 531         atom[count].attr=attr;
 532
 533         return 0;
 534 }
 535
 536 int destroy_atoms(t_moldyn *moldyn) {
 537
 538         if(moldyn->atom) free(moldyn->atom);
 539
 540         return 0;
 541 }
 542
 543 int thermal_init(t_moldyn *moldyn,u8 equi_init) {
 544
 545         /*
 546          * - gaussian distribution of velocities
 547          * - zero total momentum
 548          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 549          */
 550
 551         int i;
 552         double v,sigma;
 553         t_3dvec p_total,delta;
 554         t_atom *atom;
 555         t_random *random;
 556
 557         atom=moldyn->atom;
 558         random=&(moldyn->random);
 559
 560         printf("[moldyn] thermal init (equi init: %s)\n",equi_init?"yes":"no");
 561
 562         /* gaussian distribution of velocities */
 563         v3_zero(&p_total);
 564         for(i=0;i<moldyn->count;i++) {
 565                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
 566                 /* x direction */
 567                 v=sigma*rand_get_gauss(random);
 568                 atom[i].v.x=v;
 569                 p_total.x+=atom[i].mass*v;
 570                 /* y direction */
 571                 v=sigma*rand_get_gauss(random);
 572                 atom[i].v.y=v;
 573                 p_total.y+=atom[i].mass*v;
 574                 /* z direction */
 575                 v=sigma*rand_get_gauss(random);
 576                 atom[i].v.z=v;
 577                 p_total.z+=atom[i].mass*v;
 578         }
 579
 580         /* zero total momentum */
 581         v3_scale(&p_total,&p_total,1.0/moldyn->count);
 582         for(i=0;i<moldyn->count;i++) {
 583                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 584                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 585         }
 586
 587         /* velocity scaling */
 588         scale_velocity(moldyn,equi_init);
 589
 590         return 0;
 591 }
 592
 593 double temperature_calc(t_moldyn *moldyn) {
 594
 595         /* assume up to date kinetic energy, which is 3/2 N k_B T */
 596
 597         moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
 598
 599         return moldyn->t;
 600 }
 601
 602 double get_temperature(t_moldyn *moldyn) {
 603
 604         return moldyn->t;
 605 }
 606
 607 int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 608
 609         int i;
 610         double e,scale;
 611         t_atom *atom;
 612         int count;
 613
 614         atom=moldyn->atom;
 615
 616         /*
 617          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 618          */
 619
 620         /* get kinetic energy / temperature & count involved atoms */
 621         e=0.0;
 622         count=0;
 623         for(i=0;i<moldyn->count;i++) {
 624                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) {
 625                         e+=atom[i].mass*v3_absolute_square(&(atom[i].v));
 626                         count+=1;
 627                 }
 628         }
 629         e*=0.5;
 630         if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN);
 631         else return 0;  /* no atoms involved in scaling! */
 632
 633         /* (temporary) hack for e,t = 0 */
 634         if(e==0.0) {
 635         moldyn->t=0.0;
 636                 if(moldyn->t_ref!=0.0) {
 637                         thermal_init(moldyn,equi_init);
 638                         return 0;
 639                 }
 640                 else
 641                         return 0; /* no scaling needed */
 642         }
 643
 644
 645         /* get scaling factor */
 646         scale=moldyn->t_ref/moldyn->t;
 647         if(equi_init&TRUE)
 648                 scale*=2.0;
 649         else
 650                 if(moldyn->pt_scale&T_SCALE_BERENDSEN)
 651                         scale=1.0+(scale-1.0)/moldyn->t_tc;
 652         scale=sqrt(scale);
 653
 654         /* velocity scaling */
 655         for(i=0;i<moldyn->count;i++) {
 656                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
 657                         v3_scale(&(atom[i].v),&(atom[i].v),scale);
 658         }
 659
 660         return 0;
 661 }
 662
 663 double ideal_gas_law_pressure(t_moldyn *moldyn) {
 664
 665         double p;
 666
 667         p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume;
 668
 669         return p;
 670 }
 671
 672 double pressure_calc(t_moldyn *moldyn) {
 673
 674         int i;
 675         double v;
 676         t_virial *virial;
 677
 678         /*
 679          * P = 1/(3V) sum_i ( p_i^2 / 2m + f_i r_i )
 680          *
 681          * virial = f_i r_i
 682          */
 683
 684         v=0.0;
 685         for(i=0;i<moldyn->count;i++) {
 686                 virial=&(moldyn->atom[i].virial);
 687                 v+=(virial->xx+virial->yy+virial->zz);
 688         }
 689
 690         /* assume up to date kinetic energy */
 691         moldyn->p=2.0*moldyn->ekin+v;
 692         moldyn->p/=(3.0*moldyn->volume);
 693
 694         return moldyn->p;
 695 }
 696
 697 double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 698
 699         t_3dvec dim,*tp;
 700         double u,p;
 701         double scale;
 702         t_atom *store;
 703
 704         tp=&(moldyn->tp);
 705         store=malloc(moldyn->count*sizeof(t_atom));
 706         if(store==NULL) {
 707                 printf("[moldyn] allocating store mem failed\n");
 708                 return -1;
 709         }
 710
 711         /* save unscaled potential energy + atom/dim configuration */
 712         u=moldyn->energy;
 713         memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
 714         dim=moldyn->dim;
 715
 716         /* derivative with respect to x direction */
 717         scale=1.0+moldyn->dv/(moldyn->dim.y*moldyn->dim.z);
 718         scale_dim(moldyn,scale,TRUE,0,0);
 719         scale_atoms(moldyn,scale,TRUE,0,0);
 720         link_cell_shutdown(moldyn);
 721         link_cell_init(moldyn,QUIET);
 722         potential_force_calc(moldyn);
 723         tp->x=(moldyn->energy-u)/moldyn->dv;
 724         p=tp->x*tp->x;
 725
 726         /* restore atomic configuration + dim */
 727         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 728         moldyn->dim=dim;
 729
 730         /* derivative with respect to y direction */
 731         scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.z);
 732         scale_dim(moldyn,scale,0,TRUE,0);
 733         scale_atoms(moldyn,scale,0,TRUE,0);
 734         link_cell_shutdown(moldyn);
 735         link_cell_init(moldyn,QUIET);
 736         potential_force_calc(moldyn);
 737         tp->y=(moldyn->energy-u)/moldyn->dv;
 738         p+=tp->y*tp->y;
 739
 740         /* restore atomic configuration + dim */
 741         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 742         moldyn->dim=dim;
 743
 744         /* derivative with respect to z direction */
 745         scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.y);
 746         scale_dim(moldyn,scale,0,0,TRUE);
 747         scale_atoms(moldyn,scale,0,0,TRUE);
 748         link_cell_shutdown(moldyn);
 749         link_cell_init(moldyn,QUIET);
 750         potential_force_calc(moldyn);
 751         tp->z=(moldyn->energy-u)/moldyn->dv;
 752         p+=tp->z*tp->z;
 753
 754         /* restore atomic configuration + dim */
 755         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 756         moldyn->dim=dim;
 757
 758         printf("dU/dV komp addiert = %f %f %f\n",tp->x,tp->y,tp->z);
 759
 760         scale=1.0+pow(moldyn->dv/moldyn->volume,ONE_THIRD);
 761
 762 printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x);
 763         scale_dim(moldyn,scale,1,1,1);
 764         scale_atoms(moldyn,scale,1,1,1);
 765         link_cell_shutdown(moldyn);
 766         link_cell_init(moldyn,QUIET);
 767         potential_force_calc(moldyn);
 768 printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x);
 769
 770         printf("dU/dV einfach = %f\n",((moldyn->energy-u)/moldyn->dv)/ATM);
 771
 772         /* restore atomic configuration + dim */
 773         memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 774         moldyn->dim=dim;
 775
 776         /* restore energy */
 777         moldyn->energy=u;
 778
 779         link_cell_shutdown(moldyn);
 780         link_cell_init(moldyn,QUIET);
 781
 782         return sqrt(p);
 783 }
 784
 785 double get_pressure(t_moldyn *moldyn) {
 786
 787         return moldyn->p;
 788
 789 }
 790
 791 int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) {
 792
 793         t_3dvec *dim;
 794
 795         dim=&(moldyn->dim);
 796
 797         if(x) dim->x*=scale;
 798         if(y) dim->y*=scale;
 799         if(z) dim->z*=scale;
 800
 801         return 0;
 802 }
 803
 804 int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) {
 805
 806         int i;
 807         t_3dvec *r;
 808
 809         for(i=0;i<moldyn->count;i++) {
 810                 r=&(moldyn->atom[i].r);
 811                 if(x) r->x*=scale;
 812                 if(y) r->y*=scale;
 813                 if(z) r->z*=scale;
 814         }
 815
 816         return 0;
 817 }
 818
 819 int scale_volume(t_moldyn *moldyn) {
 820
 821         t_3dvec *dim,*vdim;
 822         double scale;
 823         t_linkcell *lc;
 824
 825         vdim=&(moldyn->vis.dim);
 826         dim=&(moldyn->dim);
 827         lc=&(moldyn->lc);
 828
 829         /* scaling factor */
 830         if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
 831                 scale=1.0-(moldyn->p_ref-moldyn->p)/moldyn->p_tc;
 832                 scale=pow(scale,ONE_THIRD);
 833         }
 834         else {
 835                 scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
 836         }
 837 moldyn->debug=scale;
 838
 839         /* scale the atoms and dimensions */
 840         scale_atoms(moldyn,scale,TRUE,TRUE,TRUE);
 841         scale_dim(moldyn,scale,TRUE,TRUE,TRUE);
 842
 843         /* visualize dimensions */
 844         if(vdim->x!=0) {
 845                 vdim->x=dim->x;
 846                 vdim->y=dim->y;
 847                 vdim->z=dim->z;
 848         }
 849
 850         /* recalculate scaled volume */
 851         moldyn->volume=dim->x*dim->y*dim->z;
 852
 853         /* adjust/reinit linkcell */
 854         if(((int)(dim->x/moldyn->cutoff)!=lc->nx)||
 855            ((int)(dim->y/moldyn->cutoff)!=lc->ny)||
 856            ((int)(dim->z/moldyn->cutoff)!=lc->nx)) {
 857                 link_cell_shutdown(moldyn);
 858                 link_cell_init(moldyn,QUIET);
 859         } else {
 860                 lc->x*=scale;
 861                 lc->y*=scale;
 862                 lc->z*=scale;
 863         }
 864
 865         return 0;
 866
 867 }
 868
 869 double get_e_kin(t_moldyn *moldyn) {
 870
 871         int i;
 872         t_atom *atom;
 873
 874         atom=moldyn->atom;
 875         moldyn->ekin=0.0;
 876
 877         for(i=0;i<moldyn->count;i++)
 878                 moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 879
 880         return moldyn->ekin;
 881 }
 882
 883 double update_e_kin(t_moldyn *moldyn) {
 884
 885         return(get_e_kin(moldyn));
 886 }
 887
 888 double get_total_energy(t_moldyn *moldyn) {
 889
 890         return(moldyn->ekin+moldyn->energy);
 891 }
 892
 893 t_3dvec get_total_p(t_moldyn *moldyn) {
 894
 895         t_3dvec p,p_total;
 896         int i;
 897         t_atom *atom;
 898
 899         atom=moldyn->atom;
 900
 901         v3_zero(&p_total);
 902         for(i=0;i<moldyn->count;i++) {
 903                 v3_scale(&p,&(atom[i].v),atom[i].mass);
 904                 v3_add(&p_total,&p_total,&p);
 905         }
 906
 907         return p_total;
 908 }
 909
 910 double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
 911
 912         double tau;
 913
 914         /* nn_dist is the nearest neighbour distance */
 915
 916         tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
 917
 918         return tau;
 919 }
 920
 921 /*
 922  * numerical tricks
 923  */
 924
 925 /* linked list / cell method */
 926
 927 int link_cell_init(t_moldyn *moldyn,u8 vol) {
 928
 929         t_linkcell *lc;
 930         int i;
 931
 932         lc=&(moldyn->lc);
 933
 934         /* partitioning the md cell */
 935         lc->nx=moldyn->dim.x/moldyn->cutoff;
 936         lc->x=moldyn->dim.x/lc->nx;
 937         lc->ny=moldyn->dim.y/moldyn->cutoff;
 938         lc->y=moldyn->dim.y/lc->ny;
 939         lc->nz=moldyn->dim.z/moldyn->cutoff;
 940         lc->z=moldyn->dim.z/lc->nz;
 941
 942         lc->cells=lc->nx*lc->ny*lc->nz;
 943         lc->subcell=malloc(lc->cells*sizeof(t_list));
 944
 945         if(lc->cells<27)
 946                 printf("[moldyn] FATAL: less then 27 subcells!\n");
 947
 948         if(vol) printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
 949
 950         for(i=0;i<lc->cells;i++)
 951                 list_init_f(&(lc->subcell[i]));
 952
 953         link_cell_update(moldyn);
 954
 955         return 0;
 956 }
 957
 958 int link_cell_update(t_moldyn *moldyn) {
 959
 960         int count,i,j,k;
 961         int nx,ny;
 962         t_atom *atom;
 963         t_linkcell *lc;
 964         double x,y,z;
 965
 966         atom=moldyn->atom;
 967         lc=&(moldyn->lc);
 968
 969         nx=lc->nx;
 970         ny=lc->ny;
 971
 972         x=moldyn->dim.x/2;
 973         y=moldyn->dim.y/2;
 974         z=moldyn->dim.z/2;
 975
 976         for(i=0;i<lc->cells;i++)
 977                 list_destroy_f(&(lc->subcell[i]));
 978
 979         for(count=0;count<moldyn->count;count++) {
 980                 i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
 981                 j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
 982                 k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
 983                 list_add_immediate_f(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]),
 984                                      &(atom[count]));
 985         }
 986
 987         return 0;
 988 }
 989
 990 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
 991
 992         t_linkcell *lc;
 993         int a;
 994         int count1,count2;
 995         int ci,cj,ck;
 996         int nx,ny,nz;
 997         int x,y,z;
 998         u8 bx,by,bz;
 999
1000         lc=&(moldyn->lc);
1001         nx=lc->nx;
1002         ny=lc->ny;
1003         nz=lc->nz;
1004         count1=1;
1005         count2=27;
1006         a=nx*ny;
1007
1008         cell[0]=lc->subcell[i+j*nx+k*a];
1009         for(ci=-1;ci<=1;ci++) {
1010                 bx=0;
1011                 x=i+ci;
1012                 if((x<0)||(x>=nx)) {
1013                         x=(x+nx)%nx;
1014                         bx=1;
1015                 }
1016                 for(cj=-1;cj<=1;cj++) {
1017                         by=0;
1018                         y=j+cj;
1019                         if((y<0)||(y>=ny)) {
1020                                 y=(y+ny)%ny;
1021                                 by=1;
1022                         }
1023                         for(ck=-1;ck<=1;ck++) {
1024                                 bz=0;
1025                                 z=k+ck;
1026                                 if((z<0)||(z>=nz)) {
1027                                         z=(z+nz)%nz;
1028                                         bz=1;
1029                                 }
1030                                 if(!(ci|cj|ck)) continue;
1031                                 if(bx|by|bz) {
1032                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
1033                                 }
1034                                 else {
1035                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
1036                                 }
1037                         }
1038                 }
1039         }
1040
1041         lc->dnlc=count1;
1042
1043         return count1;
1044 }
1045
1046 int link_cell_shutdown(t_moldyn *moldyn) {
1047
1048         int i;
1049         t_linkcell *lc;
1050
1051         lc=&(moldyn->lc);
1052
1053         for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
1054                 list_destroy_f(&(moldyn->lc.subcell[i]));
1055
1056         free(lc->subcell);
1057
1058         return 0;
1059 }
1060
1061 int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
1062
1063         int count;
1064         void *ptr;
1065         t_moldyn_schedule *schedule;
1066
1067         schedule=&(moldyn->schedule);
1068         count=++(schedule->total_sched);
1069
1070         ptr=realloc(schedule->runs,count*sizeof(int));
1071         if(!ptr) {
1072                 perror("[moldyn] realloc (runs)");
1073                 return -1;
1074         }
1075         schedule->runs=ptr;
1076         schedule->runs[count-1]=runs;
1077
1078         ptr=realloc(schedule->tau,count*sizeof(double));
1079         if(!ptr) {
1080                 perror("[moldyn] realloc (tau)");
1081                 return -1;
1082         }
1083         schedule->tau=ptr;
1084         schedule->tau[count-1]=tau;
1085
1086         printf("[moldyn] schedule added:\n");
1087         printf("  number: %d | runs: %d | tau: %f\n",count-1,runs,tau);
1088
1089
1090         return 0;
1091 }
1092
1093 int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
1094
1095         moldyn->schedule.hook=hook;
1096         moldyn->schedule.hook_params=hook_params;
1097
1098         return 0;
1099 }
1100
1101 /*
1102  *
1103  * 'integration of newtons equation' - algorithms
1104  *
1105  */
1106
1107 /* start the integration */
1108
1109 int moldyn_integrate(t_moldyn *moldyn) {
1110
1111         int i;
1112         unsigned int e,m,s,v;
1113         t_3dvec p;
1114         t_moldyn_schedule *sched;
1115         t_atom *atom;
1116         int fd;
1117         char dir[128];
1118         double ds;
1119         double energy_scale;
1120
1121         sched=&(moldyn->schedule);
1122         atom=moldyn->atom;
1123
1124         /* initialize linked cell method */
1125         link_cell_init(moldyn,VERBOSE);
1126
1127         /* logging & visualization */
1128         e=moldyn->ewrite;
1129         m=moldyn->mwrite;
1130         s=moldyn->swrite;
1131         v=moldyn->vwrite;
1132
1133         /* sqaure of some variables */
1134         moldyn->tau_square=moldyn->tau*moldyn->tau;
1135         moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
1136
1137         /* energy scaling factor */
1138         energy_scale=moldyn->count*EV;
1139
1140         /* calculate initial forces */
1141         potential_force_calc(moldyn);
1142
1143         /* some stupid checks before we actually start calculating bullshit */
1144         if(moldyn->cutoff>0.5*moldyn->dim.x)
1145                 printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
1146         if(moldyn->cutoff>0.5*moldyn->dim.y)
1147                 printf("[moldyn] warning: cutoff > 0.5 x dim.y\n");
1148         if(moldyn->cutoff>0.5*moldyn->dim.z)
1149                 printf("[moldyn] warning: cutoff > 0.5 x dim.z\n");
1150         ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
1151         if(ds>0.05*moldyn->nnd)
1152                 printf("[moldyn] warning: forces too high / tau too small!\n");
1153
1154         /* zero absolute time */
1155         moldyn->time=0.0;
1156
1157         /* debugging, ignore */
1158         moldyn->debug=0;
1159
1160         /* tell the world */
1161         printf("[moldyn] integration start, go get a coffee ...\n");
1162
1163         /* executing the schedule */
1164         for(sched->count=0;sched->count<sched->total_sched;sched->count++) {
1165
1166                 /* setting amount of runs and finite time step size */
1167                 moldyn->tau=sched->tau[sched->count];
1168                 moldyn->tau_square=moldyn->tau*moldyn->tau;
1169                 moldyn->time_steps=sched->runs[sched->count];
1170
1171         /* integration according to schedule */
1172
1173         for(i=0;i<moldyn->time_steps;i++) {
1174
1175                 /* integration step */
1176                 moldyn->integrate(moldyn);
1177
1178                 /* calculate kinetic energy, temperature and pressure */
1179                 update_e_kin(moldyn);
1180                 temperature_calc(moldyn);
1181                 pressure_calc(moldyn);
1182                 //thermodynamic_pressure_calc(moldyn);
1183
1184                 /* p/t scaling */
1185                 if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
1186                         scale_velocity(moldyn,FALSE);
1187                 if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
1188                         scale_volume(moldyn);
1189
1190                 /* check for log & visualization */
1191                 if(e) {
1192                         if(!(i%e))
1193                                 dprintf(moldyn->efd,
1194                                         "%f %f %f %f\n",
1195                                         moldyn->time,moldyn->ekin/energy_scale,
1196                                         moldyn->energy/energy_scale,
1197                                         get_total_energy(moldyn)/energy_scale);
1198                 }
1199                 if(m) {
1200                         if(!(i%m)) {
1201                                 p=get_total_p(moldyn);
1202                                 dprintf(moldyn->mfd,
1203                                         "%f %f\n",moldyn->time,v3_norm(&p));
1204                         }
1205                 }
1206                 if(s) {
1207                         if(!(i%s)) {
1208                                 snprintf(dir,128,"%s/s-%07.f.save",
1209                                          moldyn->vlsdir,moldyn->time);
1210                                 fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT);
1211                                 if(fd<0) perror("[moldyn] save fd open");
1212                                 else {
1213                                         write(fd,moldyn,sizeof(t_moldyn));
1214                                         write(fd,moldyn->atom,
1215                                               moldyn->count*sizeof(t_atom));
1216                                 }
1217                                 close(fd);
1218                         }
1219                 }
1220                 if(v) {
1221                         if(!(i%v)) {
1222                                 visual_atoms(&(moldyn->vis),moldyn->time,
1223                                              moldyn->atom,moldyn->count);
1224                                 printf("\rsched: %d, steps: %d, T: %f, P: %f V: %f",
1225                                        sched->count,i,
1226                                        moldyn->t,moldyn->p/ATM,moldyn->volume);
1227                                 fflush(stdout);
1228                         }
1229                 }
1230
1231                 /* increase absolute time */
1232                 moldyn->time+=moldyn->tau;
1233
1234         }
1235
1236                 /* check for hooks */
1237                 if(sched->hook)
1238                         sched->hook(moldyn,sched->hook_params);
1239
1240                 /* get a new info line */
1241                 printf("\n");
1242
1243         }
1244
1245         return 0;
1246 }
1247
1248 /* velocity verlet */
1249
1250 int velocity_verlet(t_moldyn *moldyn) {
1251
1252         int i,count;
1253         double tau,tau_square,h;
1254         t_3dvec delta;
1255         t_atom *atom;
1256
1257         atom=moldyn->atom;
1258         count=moldyn->count;
1259         tau=moldyn->tau;
1260         tau_square=moldyn->tau_square;
1261
1262         for(i=0;i<count;i++) {
1263                 /* new positions */
1264                 h=0.5/atom[i].mass;
1265                 v3_scale(&delta,&(atom[i].v),tau);
1266                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1267                 v3_scale(&delta,&(atom[i].f),h*tau_square);
1268                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
1269                 check_per_bound(moldyn,&(atom[i].r));
1270
1271                 /* velocities [actually v(t+tau/2)] */
1272                 v3_scale(&delta,&(atom[i].f),h*tau);
1273                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1274         }
1275
1276         /* neighbour list update */
1277         link_cell_update(moldyn);
1278
1279         /* forces depending on chosen potential */
1280         potential_force_calc(moldyn);
1281
1282         for(i=0;i<count;i++) {
1283                 /* again velocities [actually v(t+tau)] */
1284                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
1285                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
1286         }
1287
1288         return 0;
1289 }
1290
1291
1292 /*
1293  *
1294  * potentials & corresponding forces & virial routine
1295  *
1296  */
1297
1298 /* generic potential and force calculation */
1299
1300 int potential_force_calc(t_moldyn *moldyn) {
1301
1302         int i,j,k,count;
1303         t_atom *itom,*jtom,*ktom;
1304         t_virial *virial;
1305         t_linkcell *lc;
1306         t_list neighbour_i[27];
1307         t_list neighbour_i2[27];
1308         t_list *this,*that;
1309         u8 bc_ij,bc_ik;
1310         int dnlc;
1311
1312         count=moldyn->count;
1313         itom=moldyn->atom;
1314         lc=&(moldyn->lc);
1315
1316         /* reset energy */
1317         moldyn->energy=0.0;
1318
1319         /* reset force, site energy and virial of every atom */
1320         for(i=0;i<count;i++) {
1321
1322                 /* reset force */
1323                 v3_zero(&(itom[i].f));
1324
1325                 /* reset virial */
1326                 virial=(&(itom[i].virial));
1327                 virial->xx=0.0;
1328                 virial->yy=0.0;
1329                 virial->zz=0.0;
1330                 virial->xy=0.0;
1331                 virial->xz=0.0;
1332                 virial->yz=0.0;
1333
1334                 /* reset site energy */
1335                 itom[i].e=0.0;
1336
1337         }
1338
1339         /* get energy,force and virial of every atom */
1340         for(i=0;i<count;i++) {
1341
1342                 /* single particle potential/force */
1343                 if(itom[i].attr&ATOM_ATTR_1BP)
1344                         moldyn->func1b(moldyn,&(itom[i]));
1345
1346                 if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
1347                         continue;
1348
1349                 /* 2 body pair potential/force */
1350
1351                 link_cell_neighbour_index(moldyn,
1352                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
1353                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
1354                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
1355                                           neighbour_i);
1356
1357                 dnlc=lc->dnlc;
1358
1359                 for(j=0;j<27;j++) {
1360
1361                         this=&(neighbour_i[j]);
1362                         list_reset_f(this);
1363
1364                         if(this->start==NULL)
1365                                 continue;
1366
1367                         bc_ij=(j<dnlc)?0:1;
1368
1369                         do {
1370                                 jtom=this->current->data;
1371
1372                                 if(jtom==&(itom[i]))
1373                                         continue;
1374
1375                                 if((jtom->attr&ATOM_ATTR_2BP)&
1376                                    (itom[i].attr&ATOM_ATTR_2BP)) {
1377                                         moldyn->func2b(moldyn,
1378                                                        &(itom[i]),
1379                                                        jtom,
1380                                                        bc_ij);
1381                                 }
1382
1383                                 /* 3 body potential/force */
1384
1385                                 if(!(itom[i].attr&ATOM_ATTR_3BP)||
1386                                    !(jtom->attr&ATOM_ATTR_3BP))
1387                                         continue;
1388
1389                                 /* copy the neighbour lists */
1390                                 memcpy(neighbour_i2,neighbour_i,
1391                                        27*sizeof(t_list));
1392
1393                                 /* get neighbours of i */
1394                                 for(k=0;k<27;k++) {
1395
1396                                         that=&(neighbour_i2[k]);
1397                                         list_reset_f(that);
1398
1399                                         if(that->start==NULL)
1400                                                 continue;
1401
1402                                         bc_ik=(k<dnlc)?0:1;
1403
1404                                         do {
1405
1406                                                 ktom=that->current->data;
1407
1408                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
1409                                                         continue;
1410
1411                                                 if(ktom==jtom)
1412                                                         continue;
1413
1414                                                 if(ktom==&(itom[i]))
1415                                                         continue;
1416
1417                                                 moldyn->func3b(moldyn,
1418                                                                &(itom[i]),
1419                                                                jtom,
1420                                                                ktom,
1421                                                                bc_ik|bc_ij);
1422
1423                                         } while(list_next_f(that)!=\
1424                                                 L_NO_NEXT_ELEMENT);
1425
1426                                 }
1427
1428                                 /* 2bp post function */
1429                                 if(moldyn->func2b_post) {
1430                                         moldyn->func2b_post(moldyn,
1431                                                             &(itom[i]),
1432                                                             jtom,bc_ij);
1433                                 }
1434
1435                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1436
1437                 }
1438
1439         }
1440
1441 #ifdef DEBUG
1442 printf("\n\n");
1443 #endif
1444 #ifdef VDEBUG
1445 printf("\n\n");
1446 #endif
1447
1448         return 0;
1449 }
1450
1451 /*
1452  * virial calculation
1453  */
1454
1455 //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
1456 int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
1457
1458         a->virial.xx+=f->x*d->x;
1459         a->virial.yy+=f->y*d->y;
1460         a->virial.zz+=f->z*d->z;
1461         a->virial.xy+=f->x*d->y;
1462         a->virial.xz+=f->x*d->z;
1463         a->virial.yz+=f->y*d->z;
1464
1465         return 0;
1466 }
1467
1468 /*
1469  * periodic boundayr checking
1470  */
1471
1472 //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
1473 int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
1474
1475         double x,y,z;
1476         t_3dvec *dim;
1477
1478         dim=&(moldyn->dim);
1479
1480         x=dim->x/2;
1481         y=dim->y/2;
1482         z=dim->z/2;
1483
1484         if(moldyn->status&MOLDYN_STAT_PBX) {
1485                 if(a->x>=x) a->x-=dim->x;
1486                 else if(-a->x>x) a->x+=dim->x;
1487         }
1488         if(moldyn->status&MOLDYN_STAT_PBY) {
1489                 if(a->y>=y) a->y-=dim->y;
1490                 else if(-a->y>y) a->y+=dim->y;
1491         }
1492         if(moldyn->status&MOLDYN_STAT_PBZ) {
1493                 if(a->z>=z) a->z-=dim->z;
1494                 else if(-a->z>z) a->z+=dim->z;
1495         }
1496
1497         return 0;
1498 }
1499
1500 /*
1501  * debugging / critical check functions
1502  */
1503
1504 int moldyn_bc_check(t_moldyn *moldyn) {
1505
1506         t_atom *atom;
1507         t_3dvec *dim;
1508         int i;
1509         double x;
1510         u8 byte;
1511         int j,k;
1512
1513         atom=moldyn->atom;
1514         dim=&(moldyn->dim);
1515         x=dim->x/2;
1516
1517         for(i=0;i<moldyn->count;i++) {
1518                 if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {
1519                         printf("FATAL: atom %d: x: %.20f (%.20f)\n",
1520                                i,atom[i].r.x,dim->x/2);
1521                         printf("diagnostic:\n");
1522                         printf("-----------\natom.r.x:\n");
1523                         for(j=0;j<8;j++) {
1524                                 memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
1525                                 for(k=0;k<8;k++)
1526                                         printf("%d%c",
1527                                         ((byte)&(1<<k))?1:0,
1528                                         (k==7)?'\n':'|');
1529                         }
1530                         printf("---------------\nx=dim.x/2:\n");
1531                         for(j=0;j<8;j++) {
1532                                 memcpy(&byte,(u8 *)(&x)+j,1);
1533                                 for(k=0;k<8;k++)
1534                                         printf("%d%c",
1535                                         ((byte)&(1<<k))?1:0,
1536                                         (k==7)?'\n':'|');
1537                         }
1538                         if(atom[i].r.x==x) printf("the same!\n");
1539                         else printf("different!\n");
1540                 }
1541                 if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
1542                         printf("FATAL: atom %d: y: %.20f (%.20f)\n",
1543                                i,atom[i].r.y,dim->y/2);
1544                 if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
1545                         printf("FATAL: atom %d: z: %.20f (%.20f)\n",
1546                                i,atom[i].r.z,dim->z/2);
1547         }
1548
1549         return 0;
1550 }