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
  20 #include "math/math.h"
  21 #include "init/init.h"
  22 #include "random/random.h"
  23 #include "visual/visual.h"
  24 #include "list/list.h"
  25
  26
  27 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
  28
  29         //int ret;
  30
  31         //ret=moldyn_parse_argv(moldyn,argc,argv);
  32         //if(ret<0) return ret;
  33
  34         memset(moldyn,0,sizeof(t_moldyn));
  35
  36         rand_init(&(moldyn->random),NULL,1);
  37         moldyn->random.status|=RAND_STAT_VERBOSE;
  38
  39         return 0;
  40 }
  41
  42 int moldyn_shutdown(t_moldyn *moldyn) {
  43
  44         printf("[moldyn] shutdown\n");
  45         moldyn_log_shutdown(moldyn);
  46         link_cell_shutdown(moldyn);
  47         rand_close(&(moldyn->random));
  48         free(moldyn->atom);
  49
  50         return 0;
  51 }
  52
  53 int set_int_alg(t_moldyn *moldyn,u8 algo) {
  54
  55         switch(algo) {
  56                 case MOLDYN_INTEGRATE_VERLET:
  57                         moldyn->integrate=velocity_verlet;
  58                         break;
  59                 default:
  60                         printf("unknown integration algorithm: %02x\n",algo);
  61                         return -1;
  62         }
  63
  64         return 0;
  65 }
  66
  67 int set_cutoff(t_moldyn *moldyn,double cutoff) {
  68
  69         moldyn->cutoff=cutoff;
  70
  71         return 0;
  72 }
  73
  74 int set_temperature(t_moldyn *moldyn,double t_ref) {
  75
  76         moldyn->t_ref=t_ref;
  77
  78         return 0;
  79 }
  80
  81 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
  82
  83         moldyn->pt_scale=(ptype|ttype);
  84         moldyn->t_tc=ttc;
  85         moldyn->p_tc=ptc;
  86
  87         return 0;
  88 }
  89
  90 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
  91
  92         moldyn->dim.x=x;
  93         moldyn->dim.y=y;
  94         moldyn->dim.z=z;
  95
  96         if(visualize) {
  97                 moldyn->vis.dim.x=x;
  98                 moldyn->vis.dim.y=y;
  99                 moldyn->vis.dim.z=z;
 100         }
 101
 102         return 0;
 103 }
 104
 105 int set_nn_dist(t_moldyn *moldyn,double dist) {
 106
 107         moldyn->nnd=dist;
 108
 109         return 0;
 110 }
 111
 112 int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
 113
 114         if(x)
 115                 moldyn->status|=MOLDYN_STAT_PBX;
 116
 117         if(y)
 118                 moldyn->status|=MOLDYN_STAT_PBY;
 119
 120         if(z)
 121                 moldyn->status|=MOLDYN_STAT_PBZ;
 122
 123         return 0;
 124 }
 125
 126 int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params) {
 127
 128         moldyn->func1b=func;
 129         moldyn->pot1b_params=params;
 130
 131         return 0;
 132 }
 133
 134 int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) {
 135
 136         moldyn->func2b=func;
 137         moldyn->pot2b_params=params;
 138
 139         return 0;
 140 }
 141
 142 int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) {
 143
 144         moldyn->func2b_post=func;
 145         moldyn->pot2b_params=params;
 146
 147         return 0;
 148 }
 149
 150 int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) {
 151
 152         moldyn->func3b=func;
 153         moldyn->pot3b_params=params;
 154
 155         return 0;
 156 }
 157
 158 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
 159
 160         strncpy(moldyn->vlsdir,dir,127);
 161
 162         return 0;
 163 }
 164
 165 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 166
 167         char filename[128];
 168         int ret;
 169
 170         switch(type) {
 171                 case LOG_TOTAL_ENERGY:
 172                         moldyn->ewrite=timer;
 173                         snprintf(filename,127,"%s/energy",moldyn->vlsdir);
 174                         moldyn->efd=open(filename,
 175                                          O_WRONLY|O_CREAT|O_EXCL,
 176                                          S_IRUSR|S_IWUSR);
 177                         if(moldyn->efd<0) {
 178                                 perror("[moldyn] energy log fd open");
 179                                 return moldyn->efd;
 180                         }
 181                         dprintf(moldyn->efd,"# total energy log file\n");
 182                         break;
 183                 case LOG_TOTAL_MOMENTUM:
 184                         moldyn->mwrite=timer;
 185                         snprintf(filename,127,"%s/momentum",moldyn->vlsdir);
 186                         moldyn->mfd=open(filename,
 187                                          O_WRONLY|O_CREAT|O_EXCL,
 188                                          S_IRUSR|S_IWUSR);
 189                         if(moldyn->mfd<0) {
 190                                 perror("[moldyn] momentum log fd open");
 191                                 return moldyn->mfd;
 192                         }
 193                         dprintf(moldyn->efd,"# total momentum log file\n");
 194                         break;
 195                 case SAVE_STEP:
 196                         moldyn->swrite=timer;
 197                         break;
 198                 case VISUAL_STEP:
 199                         moldyn->vwrite=timer;
 200                         ret=visual_init(&(moldyn->vis),moldyn->vlsdir);
 201                         if(ret<0) {
 202                                 printf("[moldyn] visual init failure\n");
 203                                 return ret;
 204                         }
 205                         break;
 206                 default:
 207                         printf("[moldyn] unknown log mechanism: %02x\n",type);
 208                         return -1;
 209         }
 210
 211         return 0;
 212 }
 213
 214 int moldyn_log_shutdown(t_moldyn *moldyn) {
 215
 216         printf("[moldyn] log shutdown\n");
 217         if(moldyn->efd) close(moldyn->efd);
 218         if(moldyn->mfd) close(moldyn->mfd);
 219         if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
 220
 221         return 0;
 222 }
 223
 224 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 225                    u8 attr,u8 bnum,int a,int b,int c) {
 226
 227         int count;
 228         int ret;
 229         t_3dvec origin;
 230
 231         count=a*b*c;
 232
 233         if(type==FCC) count*=4;
 234
 235         if(type==DIAMOND) count*=8;
 236
 237         moldyn->atom=malloc(count*sizeof(t_atom));
 238         if(moldyn->atom==NULL) {
 239                 perror("malloc (atoms)");
 240                 return -1;
 241         }
 242
 243         v3_zero(&origin);
 244
 245         switch(type) {
 246                 case FCC:
 247                         ret=fcc_init(a,b,c,lc,moldyn->atom,&origin);
 248                         break;
 249                 case DIAMOND:
 250                         ret=diamond_init(a,b,c,lc,moldyn->atom,&origin);
 251                         break;
 252                 default:
 253                         printf("unknown lattice type (%02x)\n",type);
 254                         return -1;
 255         }
 256
 257         /* debug */
 258         if(ret!=count) {
 259                 printf("ok, there is something wrong ...\n");
 260                 printf("calculated -> %d atoms\n",count);
 261                 printf("created -> %d atoms\n",ret);
 262                 return -1;
 263         }
 264
 265         moldyn->count=count;
 266         printf("[moldyn] created lattice with %d atoms\n",count);
 267
 268         while(count) {
 269                 count-=1;
 270                 moldyn->atom[count].element=element;
 271                 moldyn->atom[count].mass=mass;
 272                 moldyn->atom[count].attr=attr;
 273                 moldyn->atom[count].bnum=bnum;
 274                 check_per_bound(moldyn,&(moldyn->atom[count].r));
 275         }
 276
 277
 278         return ret;
 279 }
 280
 281 int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
 282              t_3dvec *r,t_3dvec *v) {
 283
 284         t_atom *atom;
 285         void *ptr;
 286         int count;
 287
 288         atom=moldyn->atom;
 289         count=++(moldyn->count);
 290
 291         ptr=realloc(atom,count*sizeof(t_atom));
 292         if(!ptr) {
 293                 perror("[moldyn] realloc (add atom)");
 294                 return -1;
 295         }
 296         moldyn->atom=ptr;
 297
 298         atom=moldyn->atom;
 299         atom[count-1].r=*r;
 300         atom[count-1].v=*v;
 301         atom[count-1].element=element;
 302         atom[count-1].mass=mass;
 303         atom[count-1].bnum=bnum;
 304         atom[count-1].attr=attr;
 305
 306         return 0;
 307 }
 308
 309 int destroy_atoms(t_moldyn *moldyn) {
 310
 311         if(moldyn->atom) free(moldyn->atom);
 312
 313         return 0;
 314 }
 315
 316 int thermal_init(t_moldyn *moldyn,u8 equi_init) {
 317
 318         /*
 319          * - gaussian distribution of velocities
 320          * - zero total momentum
 321          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 322          */
 323
 324         int i;
 325         double v,sigma;
 326         t_3dvec p_total,delta;
 327         t_atom *atom;
 328         t_random *random;
 329
 330         atom=moldyn->atom;
 331         random=&(moldyn->random);
 332
 333         /* gaussian distribution of velocities */
 334         v3_zero(&p_total);
 335         for(i=0;i<moldyn->count;i++) {
 336                 sigma=sqrt(2.0*K_BOLTZMANN*moldyn->t_ref/atom[i].mass);
 337                 /* x direction */
 338                 v=sigma*rand_get_gauss(random);
 339                 atom[i].v.x=v;
 340                 p_total.x+=atom[i].mass*v;
 341                 /* y direction */
 342                 v=sigma*rand_get_gauss(random);
 343                 atom[i].v.y=v;
 344                 p_total.y+=atom[i].mass*v;
 345                 /* z direction */
 346                 v=sigma*rand_get_gauss(random);
 347                 atom[i].v.z=v;
 348                 p_total.z+=atom[i].mass*v;
 349         }
 350
 351         /* zero total momentum */
 352         v3_scale(&p_total,&p_total,1.0/moldyn->count);
 353         for(i=0;i<moldyn->count;i++) {
 354                 v3_scale(&delta,&p_total,1.0/atom[i].mass);
 355                 v3_sub(&(atom[i].v),&(atom[i].v),&delta);
 356         }
 357
 358         /* velocity scaling */
 359         scale_velocity(moldyn,equi_init);
 360
 361         return 0;
 362 }
 363
 364 int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
 365
 366         int i;
 367         double e,scale;
 368         t_atom *atom;
 369         int count;
 370
 371         atom=moldyn->atom;
 372
 373         /*
 374          * - velocity scaling (E = 3/2 N k T), E: kinetic energy
 375          */
 376
 377         /* get kinetic energy / temperature & count involved atoms */
 378         e=0.0;
 379         count=0;
 380         for(i=0;i<moldyn->count;i++) {
 381                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) {
 382                         e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 383                         count+=1;
 384                 }
 385         }
 386         if(count!=0) moldyn->t=(2.0*e)/(3.0*count*K_BOLTZMANN);
 387         else return 0;  /* no atoms involved in scaling! */
 388
 389         /* (temporary) hack for e,t = 0 */
 390         if(e==0.0) {
 391         moldyn->t=0.0;
 392                 if(moldyn->t_ref!=0.0) {
 393                         thermal_init(moldyn,equi_init);
 394                         return 0;
 395                 }
 396                 else
 397                         return 0; /* no scaling needed */
 398         }
 399
 400
 401         /* get scaling factor */
 402         scale=moldyn->t_ref/moldyn->t;
 403         if(equi_init&TRUE)
 404                 scale*=2.0;
 405         else
 406                 if(moldyn->pt_scale&T_SCALE_BERENDSEN)
 407                         scale=1.0+(scale-1.0)/moldyn->t_tc;
 408         scale=sqrt(scale);
 409
 410         /* velocity scaling */
 411         for(i=0;i<moldyn->count;i++) {
 412                 if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
 413                         v3_scale(&(atom[i].v),&(atom[i].v),scale);
 414         }
 415
 416         return 0;
 417 }
 418
 419 double get_e_kin(t_moldyn *moldyn) {
 420
 421         int i;
 422         t_atom *atom;
 423
 424         atom=moldyn->atom;
 425         moldyn->ekin=0.0;
 426
 427         for(i=0;i<moldyn->count;i++)
 428                 moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
 429
 430         return moldyn->ekin;
 431 }
 432
 433 double get_e_pot(t_moldyn *moldyn) {
 434
 435         return moldyn->energy;
 436 }
 437
 438 double update_e_kin(t_moldyn *moldyn) {
 439
 440         return(get_e_kin(moldyn));
 441 }
 442
 443 double get_total_energy(t_moldyn *moldyn) {
 444
 445         return(moldyn->ekin+moldyn->energy);
 446 }
 447
 448 t_3dvec get_total_p(t_moldyn *moldyn) {
 449
 450         t_3dvec p,p_total;
 451         int i;
 452         t_atom *atom;
 453
 454         atom=moldyn->atom;
 455
 456         v3_zero(&p_total);
 457         for(i=0;i<moldyn->count;i++) {
 458                 v3_scale(&p,&(atom[i].v),atom[i].mass);
 459                 v3_add(&p_total,&p_total,&p);
 460         }
 461
 462         return p_total;
 463 }
 464
 465 double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
 466
 467         double tau;
 468
 469         /* nn_dist is the nearest neighbour distance */
 470
 471         tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
 472
 473         return tau;
 474 }
 475
 476 /*
 477  * numerical tricks
 478  */
 479
 480 /* linked list / cell method */
 481
 482 int link_cell_init(t_moldyn *moldyn) {
 483
 484         t_linkcell *lc;
 485         int i;
 486         int fd;
 487
 488         fd=open("/dev/null",O_WRONLY);
 489
 490         lc=&(moldyn->lc);
 491
 492         /* partitioning the md cell */
 493         lc->nx=moldyn->dim.x/moldyn->cutoff;
 494         lc->x=moldyn->dim.x/lc->nx;
 495         lc->ny=moldyn->dim.y/moldyn->cutoff;
 496         lc->y=moldyn->dim.y/lc->ny;
 497         lc->nz=moldyn->dim.z/moldyn->cutoff;
 498         lc->z=moldyn->dim.z/lc->nz;
 499
 500         lc->cells=lc->nx*lc->ny*lc->nz;
 501         lc->subcell=malloc(lc->cells*sizeof(t_list));
 502
 503         printf("[moldyn] initializing linked cells (%d)\n",lc->cells);
 504
 505         for(i=0;i<lc->cells;i++)
 506                 //list_init(&(lc->subcell[i]),1);
 507                 list_init(&(lc->subcell[i]),fd);
 508
 509         link_cell_update(moldyn);
 510
 511         return 0;
 512 }
 513
 514 int link_cell_update(t_moldyn *moldyn) {
 515
 516         int count,i,j,k;
 517         int nx,ny,nz;
 518         t_atom *atom;
 519         t_linkcell *lc;
 520
 521         atom=moldyn->atom;
 522         lc=&(moldyn->lc);
 523
 524         nx=lc->nx;
 525         ny=lc->ny;
 526         nz=lc->nz;
 527
 528         for(i=0;i<lc->cells;i++)
 529                 list_destroy(&(moldyn->lc.subcell[i]));
 530
 531         for(count=0;count<moldyn->count;count++) {
 532                 i=(atom[count].r.x+(moldyn->dim.x/2))/lc->x;
 533                 j=(atom[count].r.y+(moldyn->dim.y/2))/lc->y;
 534                 k=(atom[count].r.z+(moldyn->dim.z/2))/lc->z;
 535                 list_add_immediate_ptr(&(moldyn->lc.subcell[i+j*nx+k*nx*ny]),
 536                                        &(atom[count]));
 537         }
 538
 539         return 0;
 540 }
 541
 542 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
 543
 544         t_linkcell *lc;
 545         int a;
 546         int count1,count2;
 547         int ci,cj,ck;
 548         int nx,ny,nz;
 549         int x,y,z;
 550         u8 bx,by,bz;
 551
 552         lc=&(moldyn->lc);
 553         nx=lc->nx;
 554         ny=lc->ny;
 555         nz=lc->nz;
 556         count1=1;
 557         count2=27;
 558         a=nx*ny;
 559
 560         cell[0]=lc->subcell[i+j*nx+k*a];
 561         for(ci=-1;ci<=1;ci++) {
 562                 bx=0;
 563                 x=i+ci;
 564                 if((x<0)||(x>=nx)) {
 565                         x=(x+nx)%nx;
 566                         bx=1;
 567                 }
 568                 for(cj=-1;cj<=1;cj++) {
 569                         by=0;
 570                         y=j+cj;
 571                         if((y<0)||(y>=ny)) {
 572                                 y=(y+ny)%ny;
 573                                 by=1;
 574                         }
 575                         for(ck=-1;ck<=1;ck++) {
 576                                 bz=0;
 577                                 z=k+ck;
 578                                 if((z<0)||(z>=nz)) {
 579                                         z=(z+nz)%nz;
 580                                         bz=1;
 581                                 }
 582                                 if(!(ci|cj|ck)) continue;
 583                                 if(bx|by|bz) {
 584                                         cell[--count2]=lc->subcell[x+y*nx+z*a];
 585                                 }
 586                                 else {
 587                                         cell[count1++]=lc->subcell[x+y*nx+z*a];
 588                                 }
 589                         }
 590                 }
 591         }
 592
 593         lc->dnlc=count1;
 594
 595         return count1;
 596 }
 597
 598 int link_cell_shutdown(t_moldyn *moldyn) {
 599
 600         int i;
 601         t_linkcell *lc;
 602
 603         lc=&(moldyn->lc);
 604
 605         for(i=0;i<lc->nx*lc->ny*lc->nz;i++)
 606                 list_shutdown(&(moldyn->lc.subcell[i]));
 607
 608         return 0;
 609 }
 610
 611 int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
 612
 613         int count;
 614         void *ptr;
 615         t_moldyn_schedule *schedule;
 616
 617         schedule=&(moldyn->schedule);
 618         count=++(schedule->content_count);
 619
 620         ptr=realloc(moldyn->schedule.runs,count*sizeof(int));
 621         if(!ptr) {
 622                 perror("[moldyn] realloc (runs)");
 623                 return -1;
 624         }
 625         moldyn->schedule.runs=ptr;
 626         moldyn->schedule.runs[count-1]=runs;
 627
 628         ptr=realloc(schedule->tau,count*sizeof(double));
 629         if(!ptr) {
 630                 perror("[moldyn] realloc (tau)");
 631                 return -1;
 632         }
 633         moldyn->schedule.tau=ptr;
 634         moldyn->schedule.tau[count-1]=tau;
 635
 636         return 0;
 637 }
 638
 639 int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params) {
 640
 641         moldyn->schedule.hook=hook;
 642         moldyn->schedule.hook_params=hook_params;
 643
 644         return 0;
 645 }
 646
 647 /*
 648  *
 649  * 'integration of newtons equation' - algorithms
 650  *
 651  */
 652
 653 /* start the integration */
 654
 655 int moldyn_integrate(t_moldyn *moldyn) {
 656
 657         int i,sched;
 658         unsigned int e,m,s,v;
 659         t_3dvec p;
 660         t_moldyn_schedule *schedule;
 661         t_atom *atom;
 662         int fd;
 663         char dir[128];
 664         double ds;
 665
 666         schedule=&(moldyn->schedule);
 667         atom=moldyn->atom;
 668
 669         /* initialize linked cell method */
 670         link_cell_init(moldyn);
 671
 672         /* logging & visualization */
 673         e=moldyn->ewrite;
 674         m=moldyn->mwrite;
 675         s=moldyn->swrite;
 676         v=moldyn->vwrite;
 677
 678         /* sqaure of some variables */
 679         moldyn->tau_square=moldyn->tau*moldyn->tau;
 680         moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
 681
 682         /* calculate initial forces */
 683         potential_force_calc(moldyn);
 684
 685         /* some stupid checks before we actually start calculating bullshit */
 686         if(moldyn->cutoff>0.5*moldyn->dim.x)
 687                 printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
 688         if(moldyn->cutoff>0.5*moldyn->dim.y)
 689                 printf("[moldyn] warning: cutoff > 0.5 x dim.y\n");
 690         if(moldyn->cutoff>0.5*moldyn->dim.z)
 691                 printf("[moldyn] warning: cutoff > 0.5 x dim.z\n");
 692         ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
 693         if(ds>0.05*moldyn->nnd)
 694                 printf("[moldyn] warning: forces too high / tau too small!\n");
 695
 696         /* zero absolute time */
 697         moldyn->time=0.0;
 698
 699         /* debugging, ignore */
 700         moldyn->debug=0;
 701
 702         /* executing the schedule */
 703         for(sched=0;sched<moldyn->schedule.content_count;sched++) {
 704
 705                 /* setting amount of runs and finite time step size */
 706                 moldyn->tau=schedule->tau[sched];
 707                 moldyn->tau_square=moldyn->tau*moldyn->tau;
 708                 moldyn->time_steps=schedule->runs[sched];
 709
 710         /* integration according to schedule */
 711
 712         for(i=0;i<moldyn->time_steps;i++) {
 713
 714                 /* integration step */
 715                 moldyn->integrate(moldyn);
 716
 717                 /* p/t scaling */
 718                 if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
 719                         scale_velocity(moldyn,FALSE);
 720
 721                 /* check for log & visualization */
 722                 if(e) {
 723                         if(!(i%e))
 724                                 dprintf(moldyn->efd,
 725                                         "%f %f %f %f\n",
 726                                         moldyn->time,update_e_kin(moldyn),
 727                                         moldyn->energy,
 728                                         get_total_energy(moldyn));
 729                 }
 730                 if(m) {
 731                         if(!(i%m)) {
 732                                 p=get_total_p(moldyn);
 733                                 dprintf(moldyn->mfd,
 734                                         "%f %f\n",moldyn->time,v3_norm(&p));
 735                         }
 736                 }
 737                 if(s) {
 738                         if(!(i%s)) {
 739                                 snprintf(dir,128,"%s/s-%07.f.save",
 740                                          moldyn->vlsdir,moldyn->time);
 741                                 fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT);
 742                                 if(fd<0) perror("[moldyn] save fd open");
 743                                 else {
 744                                         write(fd,moldyn,sizeof(t_moldyn));
 745                                         write(fd,moldyn->atom,
 746                                               moldyn->count*sizeof(t_atom));
 747                                 }
 748                                 close(fd);
 749                         }
 750                 }
 751                 if(v) {
 752                         if(!(i%v)) {
 753                                 visual_atoms(&(moldyn->vis),moldyn->time,
 754                                              moldyn->atom,moldyn->count);
 755                                 printf("\rsched: %d, steps: %d, debug: %d",
 756                                        sched,i,moldyn->debug);
 757                                 fflush(stdout);
 758                         }
 759                 }
 760
 761                 /* increase absolute time */
 762                 moldyn->time+=moldyn->tau;
 763
 764         }
 765
 766                 /* check for hooks */
 767                 if(schedule->hook)
 768                         schedule->hook(moldyn,schedule->hook_params);
 769
 770                 /* get a new info line */
 771                 printf("\n");
 772
 773         }
 774
 775         return 0;
 776 }
 777
 778 /* velocity verlet */
 779
 780 int velocity_verlet(t_moldyn *moldyn) {
 781
 782         int i,count;
 783         double tau,tau_square;
 784         t_3dvec delta;
 785         t_atom *atom;
 786
 787         atom=moldyn->atom;
 788         count=moldyn->count;
 789         tau=moldyn->tau;
 790         tau_square=moldyn->tau_square;
 791
 792         for(i=0;i<count;i++) {
 793                 /* new positions */
 794                 v3_scale(&delta,&(atom[i].v),tau);
 795                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 796                 v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
 797                 v3_add(&(atom[i].r),&(atom[i].r),&delta);
 798                 check_per_bound(moldyn,&(atom[i].r));
 799
 800                 /* velocities */
 801                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 802                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 803         }
 804
 805         /* neighbour list update */
 806         link_cell_update(moldyn);
 807
 808         /* forces depending on chosen potential */
 809         potential_force_calc(moldyn);
 810
 811         for(i=0;i<count;i++) {
 812                 /* again velocities */
 813                 v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
 814                 v3_add(&(atom[i].v),&(atom[i].v),&delta);
 815         }
 816
 817         return 0;
 818 }
 819
 820
 821 /*
 822  *
 823  * potentials & corresponding forces
 824  *
 825  */
 826
 827 /* generic potential and force calculation */
 828
 829 int potential_force_calc(t_moldyn *moldyn) {
 830
 831         int i,j,k,count;
 832         t_atom *itom,*jtom,*ktom;
 833         t_linkcell *lc;
 834         t_list neighbour_i[27];
 835         t_list neighbour_i2[27];
 836         //t_list neighbour_j[27];
 837         t_list *this,*that;
 838         u8 bc_ij,bc_ik;
 839         int dnlc;
 840
 841         count=moldyn->count;
 842         itom=moldyn->atom;
 843         lc=&(moldyn->lc);
 844
 845         /* reset energy */
 846         moldyn->energy=0.0;
 847
 848         /* get energy and force of every atom */
 849         for(i=0;i<count;i++) {
 850
 851                 /* reset force */
 852                 v3_zero(&(itom[i].f));
 853
 854                 /* single particle potential/force */
 855                 if(itom[i].attr&ATOM_ATTR_1BP)
 856                         moldyn->func1b(moldyn,&(itom[i]));
 857
 858                 if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
 859                         continue;
 860
 861                 /* 2 body pair potential/force */
 862
 863                 link_cell_neighbour_index(moldyn,
 864                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
 865                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
 866                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
 867                                           neighbour_i);
 868
 869                 dnlc=lc->dnlc;
 870
 871                 for(j=0;j<27;j++) {
 872
 873                         this=&(neighbour_i[j]);
 874                         list_reset(this);
 875
 876                         if(this->start==NULL)
 877                                 continue;
 878
 879                         bc_ij=(j<dnlc)?0:1;
 880
 881                         do {
 882                                 jtom=this->current->data;
 883
 884                                 if(jtom==&(itom[i]))
 885                                         continue;
 886
 887                                 if((jtom->attr&ATOM_ATTR_2BP)&
 888                                    (itom[i].attr&ATOM_ATTR_2BP))
 889                                         moldyn->func2b(moldyn,
 890                                                        &(itom[i]),
 891                                                        jtom,
 892                                                        bc_ij);
 893
 894                                 /* 3 body potential/force */
 895
 896                                 if(!(itom[i].attr&ATOM_ATTR_3BP)||
 897                                    !(jtom->attr&ATOM_ATTR_3BP))
 898                                         continue;
 899
 900                                 /* copy the neighbour lists */
 901                                 memcpy(neighbour_i2,neighbour_i,
 902                                        27*sizeof(t_list));
 903
 904                                 /* get neighbours of i */
 905                                 for(k=0;k<27;k++) {
 906
 907                                         that=&(neighbour_i2[k]);
 908                                         list_reset(that);
 909
 910                                         if(that->start==NULL)
 911                                                 continue;
 912
 913                                         bc_ik=(k<dnlc)?0:1;
 914
 915                                         do {
 916
 917                                                 ktom=that->current->data;
 918
 919                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
 920                                                         continue;
 921
 922                                                 if(ktom==jtom)
 923                                                         continue;
 924
 925                                                 if(ktom==&(itom[i]))
 926                                                         continue;
 927
 928                                                 moldyn->func3b(moldyn,
 929                                                                &(itom[i]),
 930                                                                jtom,
 931                                                                ktom,
 932                                                                bc_ik|bc_ij);
 933
 934                                         } while(list_next(that)!=\
 935                                                 L_NO_NEXT_ELEMENT);
 936
 937                                 }
 938
 939                                 /* 2bp post function */
 940                                 if(moldyn->func2b_post) {
 941                                         moldyn->func2b_post(moldyn,
 942                                                             &(itom[i]),
 943                                                             jtom,bc_ij);
 944                                 }
 945
 946                         } while(list_next(this)!=L_NO_NEXT_ELEMENT);
 947
 948                 }
 949
 950         }
 951
 952         return 0;
 953 }
 954
 955 /*
 956  * periodic boundayr checking
 957  */
 958
 959 int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
 960
 961         double x,y,z;
 962         t_3dvec *dim;
 963
 964         dim=&(moldyn->dim);
 965
 966         x=0.5*dim->x;
 967         y=0.5*dim->y;
 968         z=0.5*dim->z;
 969
 970         if(moldyn->status&MOLDYN_STAT_PBX) {
 971                 if(a->x>=x) a->x-=dim->x;
 972                 else if(-a->x>x) a->x+=dim->x;
 973         }
 974         if(moldyn->status&MOLDYN_STAT_PBY) {
 975                 if(a->y>=y) a->y-=dim->y;
 976                 else if(-a->y>y) a->y+=dim->y;
 977         }
 978         if(moldyn->status&MOLDYN_STAT_PBZ) {
 979                 if(a->z>=z) a->z-=dim->z;
 980                 else if(-a->z>z) a->z+=dim->z;
 981         }
 982
 983         return 0;
 984 }
 985
 986
 987 /*
 988  * example potentials
 989  */
 990
 991 /* harmonic oscillator potential and force */
 992
 993 int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 994
 995         t_ho_params *params;
 996         t_3dvec force,distance;
 997         double d;
 998         double sc,equi_dist;
 999
1000         params=moldyn->pot2b_params;
1001         sc=params->spring_constant;
1002         equi_dist=params->equilibrium_distance;
1003
1004         v3_sub(&distance,&(aj->r),&(ai->r));
1005
1006         if(bc) check_per_bound(moldyn,&distance);
1007         d=v3_norm(&distance);
1008         if(d<=moldyn->cutoff) {
1009                 /* energy is 1/2 (d-d0)^2, but we will add this twice ... */
1010                 moldyn->energy+=(0.25*sc*(d-equi_dist)*(d-equi_dist));
1011                 /* f = -grad E; grad r_ij = -1 1/r_ij distance */
1012                 v3_scale(&force,&distance,sc*(1.0-(equi_dist/d)));
1013                 v3_add(&(ai->f),&(ai->f),&force);
1014         }
1015
1016         return 0;
1017 }
1018
1019 /* lennard jones potential & force for one sort of atoms */
1020
1021 int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
1022
1023         t_lj_params *params;
1024         t_3dvec force,distance;
1025         double d,h1,h2;
1026         double eps,sig6,sig12;
1027
1028         params=moldyn->pot2b_params;
1029         eps=params->epsilon4;
1030         sig6=params->sigma6;
1031         sig12=params->sigma12;
1032
1033         v3_sub(&distance,&(aj->r),&(ai->r));
1034         if(bc) check_per_bound(moldyn,&distance);
1035         d=v3_absolute_square(&distance);        /* 1/r^2 */
1036         if(d<=moldyn->cutoff_square) {
1037                 d=1.0/d;                        /* 1/r^2 */
1038                 h2=d*d;                         /* 1/r^4 */
1039                 h2*=d;                          /* 1/r^6 */
1040                 h1=h2*h2;                       /* 1/r^12 */
1041                 /* energy is eps*..., but we will add this twice ... */
1042                 moldyn->energy+=0.5*eps*(sig12*h1-sig6*h2);
1043                 h2*=d;                          /* 1/r^8 */
1044                 h1*=d;                          /* 1/r^14 */
1045                 h2*=6*sig6;
1046                 h1*=12*sig12;
1047                 d=+h1-h2;
1048                 d*=eps;
1049                 v3_scale(&force,&distance,-1.0*d); /* f = - grad E */
1050                 v3_add(&(ai->f),&(ai->f),&force);
1051         }
1052
1053         return 0;
1054 }
1055
1056 /*
1057  * tersoff potential & force for 2 sorts of atoms
1058  */
1059
1060 /* create mixed terms from parameters and set them */
1061 int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
1062
1063         printf("[moldyn] tersoff parameter completion\n");
1064         p->Smixed=sqrt(p->S[0]*p->S[1]);
1065         p->Rmixed=sqrt(p->R[0]*p->R[1]);
1066         p->Amixed=sqrt(p->A[0]*p->A[1]);
1067         p->Bmixed=sqrt(p->B[0]*p->B[1]);
1068         p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
1069         p->mu_m=0.5*(p->mu[0]+p->mu[1]);
1070
1071         printf("[moldyn] tersoff mult parameter info:\n");
1072         printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
1073         printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
1074         printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
1075         printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
1076         printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],
1077                                           p->lambda_m);
1078         printf("  mu     | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m);
1079         printf("  beta   | %.10f | %.10f\n",p->beta[0],p->beta[1]);
1080         printf("  n      | %f | %f\n",p->n[0],p->n[1]);
1081         printf("  c      | %f | %f\n",p->c[0],p->c[1]);
1082         printf("  d      | %f | %f\n",p->d[0],p->d[1]);
1083         printf("  h      | %f | %f\n",p->h[0],p->h[1]);
1084         printf("  chi    | %f \n",p->chi);
1085
1086         return 0;
1087 }
1088
1089 /* tersoff 1 body part */
1090 int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
1091
1092         int num;
1093         t_tersoff_mult_params *params;
1094         t_tersoff_exchange *exchange;
1095
1096         num=ai->bnum;
1097         params=moldyn->pot1b_params;
1098         exchange=&(params->exchange);
1099
1100         /*
1101          * simple: point constant parameters only depending on atom i to
1102          *         their right values
1103          */
1104
1105         exchange->beta_i=&(params->beta[num]);
1106         exchange->n_i=&(params->n[num]);
1107         exchange->c_i=&(params->c[num]);
1108         exchange->d_i=&(params->d[num]);
1109         exchange->h_i=&(params->h[num]);
1110
1111         exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
1112         exchange->ci2=params->c[num]*params->c[num];
1113         exchange->di2=params->d[num]*params->d[num];
1114         exchange->ci2di2=exchange->ci2/exchange->di2;
1115
1116         return 0;
1117 }
1118
1119 /* tersoff 2 body part */
1120 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
1121
1122         t_tersoff_mult_params *params;
1123         t_tersoff_exchange *exchange;
1124         t_3dvec dist_ij,force;
1125         double d_ij;
1126         double A,B,R,S,lambda,mu;
1127         double f_r,df_r;
1128         double f_c,df_c;
1129         int num;
1130         double s_r;
1131         double arg;
1132
1133         params=moldyn->pot2b_params;
1134         num=aj->bnum;
1135         exchange=&(params->exchange);
1136
1137         /* clear 3bp and 2bp post run */
1138         exchange->run3bp=0;
1139         exchange->run2bp_post=0;
1140
1141         /* reset S > r > R mark */
1142         exchange->d_ij_between_rs=0;
1143
1144         /*
1145          * calc of 2bp contribution of V_ij and dV_ij/ji
1146          *
1147          * for Vij and dV_ij we need:
1148          * - f_c_ij, df_c_ij
1149          * - f_r_ij, df_r_ij
1150          *
1151          * for dV_ji we need:
1152          * - f_c_ji = f_c_ij, df_c_ji = df_c_ij
1153          * - f_r_ji = f_r_ij; df_r_ji = df_r_ij
1154          *
1155          */
1156
1157         /* dist_ij, d_ij */
1158         v3_sub(&dist_ij,&(aj->r),&(ai->r));
1159         if(bc) check_per_bound(moldyn,&dist_ij);
1160         d_ij=v3_norm(&dist_ij);
1161
1162         /* save for use in 3bp */
1163         exchange->d_ij=d_ij;
1164         exchange->dist_ij=dist_ij;
1165
1166         /* constants */
1167         if(num==ai->bnum) {
1168                 S=params->S[num];
1169                 R=params->R[num];
1170                 A=params->A[num];
1171                 B=params->B[num];
1172                 lambda=params->lambda[num];
1173                 mu=params->mu[num];
1174                 exchange->chi=1.0;
1175         }
1176         else {
1177                 S=params->Smixed;
1178                 R=params->Rmixed;
1179                 A=params->Amixed;
1180                 B=params->Bmixed;
1181                 lambda=params->lambda_m;
1182                 mu=params->mu_m;
1183                 params->exchange.chi=params->chi;
1184         }
1185
1186         /* if d_ij > S => no force & potential energy contribution */
1187         if(d_ij>S)
1188                 return 0;
1189
1190         /* more constants */
1191         exchange->beta_j=&(params->beta[num]);
1192         exchange->n_j=&(params->n[num]);
1193         exchange->c_j=&(params->c[num]);
1194         exchange->d_j=&(params->d[num]);
1195         exchange->h_j=&(params->h[num]);
1196         if(num==ai->bnum) {
1197                 exchange->betajnj=exchange->betaini;
1198                 exchange->cj2=exchange->ci2;
1199                 exchange->dj2=exchange->di2;
1200                 exchange->cj2dj2=exchange->ci2di2;
1201         }
1202         else {
1203                 exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
1204                 exchange->cj2=params->c[num]*params->c[num];
1205                 exchange->dj2=params->d[num]*params->d[num];
1206                 exchange->cj2dj2=exchange->cj2/exchange->dj2;
1207         }
1208
1209         /* f_r_ij = f_r_ji, df_r_ij = df_r_ji */
1210         f_r=A*exp(-lambda*d_ij);
1211         df_r=lambda*f_r/d_ij;
1212
1213         /* f_a, df_a calc (again, same for ij and ji) | save for later use! */
1214         exchange->f_a=-B*exp(-mu*d_ij);
1215         exchange->df_a=-mu*exchange->f_a/d_ij;
1216
1217         /* f_c, df_c calc (again, same for ij and ji) */
1218         if(d_ij<R) {
1219                 /* f_c = 1, df_c = 0 */
1220                 f_c=1.0;
1221                 df_c=0.0;
1222                 /* two body contribution (ij, ji) */
1223                 v3_scale(&force,&dist_ij,-df_r);
1224         }
1225         else {
1226                 s_r=S-R;
1227                 arg=M_PI*(d_ij-R)/s_r;
1228                 f_c=0.5+0.5*cos(arg);
1229                 df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij));
1230                 /* two body contribution (ij, ji) */
1231                 v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
1232                 /* tell 3bp that S > r > R */
1233                 exchange->d_ij_between_rs=1;
1234         }
1235
1236         /* add forces of 2bp (ij, ji) contribution
1237          * dVij = dVji and we sum up both: no 1/2) */
1238         v3_add(&(ai->f),&(ai->f),&force);
1239
1240         /* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */
1241         moldyn->energy+=(0.5*f_r*f_c);
1242
1243         /* save for use in 3bp */
1244         exchange->f_c=f_c;
1245         exchange->df_c=df_c;
1246
1247         /* enable the run of 3bp function and 2bp post processing */
1248         exchange->run3bp=1;
1249         exchange->run2bp_post=1;
1250
1251         /* reset 3bp sums */
1252         exchange->zeta_ij=0.0;
1253         exchange->zeta_ji=0.0;
1254         v3_zero(&(exchange->dzeta_ij));
1255         v3_zero(&(exchange->dzeta_ji));
1256
1257         return 0;
1258 }
1259
1260 /* tersoff 2 body post part */
1261
1262 int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
1263
1264         /*
1265          * here we have to allow for the 3bp sums
1266          *
1267          * that is:
1268          * - zeta_ij, dzeta_ij
1269          * - zeta_ji, dzeta_ji
1270          *
1271          * to compute the 3bp contribution to:
1272          * - Vij, dVij
1273          * - dVji
1274          *
1275          */
1276
1277         t_tersoff_mult_params *params;
1278         t_tersoff_exchange *exchange;
1279
1280         t_3dvec force,temp;
1281         t_3dvec *dist_ij;
1282         double b,db,tmp;
1283         double f_c,df_c,f_a,df_a;
1284         double chi,ni,betaini,nj,betajnj;
1285         double zeta;
1286
1287         params=moldyn->pot2b_params;
1288         exchange=&(params->exchange);
1289
1290         /* we do not run if f_c_ij was detected to be 0! */
1291         if(!(exchange->run2bp_post))
1292                 return 0;
1293
1294         f_c=exchange->f_c;
1295         df_c=exchange->df_c;
1296         f_a=exchange->f_a;
1297         df_a=exchange->df_a;
1298         betaini=exchange->betaini;
1299         betajnj=exchange->betajnj;
1300         ni=*(exchange->n_i);
1301         nj=*(exchange->n_j);
1302         chi=exchange->chi;
1303         dist_ij=&(exchange->dist_ij);
1304
1305         /* Vij and dVij */
1306         zeta=exchange->zeta_ij;
1307         if(zeta==0.0) {
1308                 moldyn->debug++;                /* just for debugging ... */
1309                 db=0.0;
1310                 b=chi;
1311                 v3_scale(&force,dist_ij,df_a*b*f_c);
1312         }
1313         else {
1314                 tmp=betaini*pow(zeta,ni-1.0);           /* beta^n * zeta^n-1 */
1315                 b=(1+zeta*tmp);                         /* 1 + beta^n zeta^n */
1316                 db=chi*pow(b,-1.0/(2*ni)-1);            /* x(...)^(-1/2n - 1) */
1317                 b=db*b;                                 /* b_ij */
1318                 db*=-0.5*tmp;                           /* db_ij */
1319                 v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
1320                 v3_scale(&temp,dist_ij,df_a*b);
1321                 v3_add(&force,&force,&temp);
1322                 v3_scale(&force,&force,f_c);
1323         }
1324         v3_scale(&temp,dist_ij,df_c*b*f_a);
1325         v3_add(&force,&force,&temp);
1326         v3_scale(&force,&force,-0.5);
1327
1328         /* add force */
1329         v3_add(&(ai->f),&(ai->f),&force);
1330
1331         /* add energy of 3bp sum */
1332         moldyn->energy+=(0.5*f_c*b*f_a);
1333
1334         /* dVji */
1335         zeta=exchange->zeta_ji;
1336         if(zeta==0.0) {
1337                 moldyn->debug++;
1338                 b=chi;
1339                 v3_scale(&force,dist_ij,df_a*b*f_c);
1340         }
1341         else {
1342                 tmp=betajnj*pow(zeta,nj-1.0);           /* beta^n * zeta^n-1 */
1343                 b=(1+zeta*tmp);                         /* 1 + beta^n zeta^n */
1344                 db=chi*pow(b,-1.0/(2*nj)-1);            /* x(...)^(-1/2n - 1) */
1345                 b=db*b;                                 /* b_ij */
1346                 db*=-0.5*tmp;                           /* db_ij */
1347                 v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
1348                 v3_scale(&temp,dist_ij,df_a*b);
1349                 v3_add(&force,&force,&temp);
1350                 v3_scale(&force,&force,f_c);
1351         }
1352         v3_scale(&temp,dist_ij,df_c*b*f_a);
1353         v3_add(&force,&force,&temp);
1354         v3_scale(&force,&force,-0.5);
1355
1356         /* add force */
1357         v3_sub(&(ai->f),&(ai->f),&force);
1358
1359         return 0;
1360 }
1361
1362 /* tersoff 3 body part */
1363
1364 int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
1365
1366         t_tersoff_mult_params *params;
1367         t_tersoff_exchange *exchange;
1368         t_3dvec dist_ij,dist_ik,dist_jk;
1369         t_3dvec temp1,temp2;
1370         t_3dvec *dzeta;
1371         double R,S,s_r;
1372         double B,mu;
1373         double d_ij,d_ik,d_jk;
1374         double rr,dd;
1375         double f_c,df_c;
1376         double f_c_ik,df_c_ik,arg;
1377         double f_c_jk;
1378         double n,c,d,h;
1379         double c2,d2,c2d2;
1380         double cos_theta,d_costheta1,d_costheta2;
1381         double h_cos,d2_h_cos2;
1382         double frac,g,zeta,chi;
1383         double tmp;
1384         int num;
1385
1386         params=moldyn->pot3b_params;
1387         exchange=&(params->exchange);
1388
1389         if(!(exchange->run3bp))
1390                 return 0;
1391
1392         /*
1393          * calc of 3bp contribution of V_ij and dV_ij/ji/jk &
1394          * 2bp contribution of dV_jk
1395          *
1396          * for Vij and dV_ij we still need:
1397          * - b_ij, db_ij (zeta_ij)
1398          *   - f_c_ik, df_c_ik, constants_i, cos_theta_ijk, d_costheta_ijk
1399          *
1400          * for dV_ji we still need:
1401          * - b_ji, db_ji (zeta_ji)
1402          *   - f_c_jk, d_c_jk, constants_j, cos_theta_jik, d_costheta_jik
1403          *
1404          * for dV_jk we need:
1405          * - f_c_jk
1406          * - f_a_jk
1407          * - db_jk (zeta_jk)
1408          *   - f_c_ji, df_c_ji, constants_j, cos_theta_jki, d_costheta_jki
1409          *
1410          */
1411
1412         /*
1413          * get exchange data
1414          */
1415
1416         /* dist_ij, d_ij - this is < S_ij ! */
1417         dist_ij=exchange->dist_ij;
1418         d_ij=exchange->d_ij;
1419
1420         /* f_c_ij, df_c_ij (same for ji) */
1421         f_c=exchange->f_c;
1422         df_c=exchange->df_c;
1423
1424         /*
1425          * calculate unknown values now ...
1426          */
1427
1428         /* V_ij and dV_ij stuff (in b_ij there is f_c_ik) */
1429
1430         /* dist_ik, d_ik */
1431         v3_sub(&dist_ik,&(ak->r),&(ai->r));
1432         if(bc) check_per_bound(moldyn,&dist_ik);
1433         d_ik=v3_norm(&dist_ik);
1434
1435         /* ik constants */
1436         num=ai->bnum;
1437         if(num==ak->bnum) {
1438                 R=params->R[num];
1439                 S=params->S[num];
1440         }
1441         else {
1442                 R=params->Rmixed;
1443                 S=params->Smixed;
1444         }
1445
1446         /* zeta_ij/dzeta_ij contribution only for d_ik < S */
1447         if(d_ik<S) {
1448
1449                 /* get constants_i from exchange data */
1450                 n=*(exchange->n_i);
1451                 c=*(exchange->c_i);
1452                 d=*(exchange->d_i);
1453                 h=*(exchange->h_i);
1454                 c2=exchange->ci2;
1455                 d2=exchange->di2;
1456                 c2d2=exchange->ci2di2;
1457
1458                 /* cosine of theta_ijk by scalaproduct */
1459                 rr=v3_scalar_product(&dist_ij,&dist_ik);
1460                 dd=d_ij*d_ik;
1461                 cos_theta=rr/dd;
1462
1463                 /* d_costheta */
1464                 tmp=1.0/dd;
1465                 d_costheta1=cos_theta/(d_ij*d_ij)-tmp;
1466                 d_costheta2=cos_theta/(d_ik*d_ik)-tmp;
1467
1468                 /* some usefull values */
1469                 h_cos=(h-cos_theta);
1470                 d2_h_cos2=d2+(h_cos*h_cos);
1471                 frac=c2/(d2_h_cos2);
1472
1473                 /* g(cos_theta) */
1474                 g=1.0+c2d2-frac;
1475
1476                 /* d_costheta_ij and dg(cos_theta) - needed in any case! */
1477                 v3_scale(&temp1,&dist_ij,d_costheta1);
1478                 v3_scale(&temp2,&dist_ik,d_costheta2);
1479                 v3_add(&temp1,&temp1,&temp2);
1480                 v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
1481
1482                 /* f_c_ik & df_c_ik + {d,}zeta contribution */
1483                 dzeta=&(exchange->dzeta_ij);
1484                 if(d_ik<R) {
1485                         /* {d,}f_c_ik */
1486                         // => f_c_ik=1.0;
1487                         // => df_c_ik=0.0; of course we do not set this!
1488
1489                         /* zeta_ij */
1490                         exchange->zeta_ij+=g;
1491
1492                         /* dzeta_ij */
1493                         v3_add(dzeta,dzeta,&temp1);
1494                 }
1495                 else {
1496                         /* {d,}f_c_ik */
1497                         s_r=S-R;
1498                         arg=M_PI*(d_ik-R)/s_r;
1499                         f_c_ik=0.5+0.5*cos(arg);
1500                         df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik));
1501
1502                         /* zeta_ij */
1503                         exchange->zeta_ij+=f_c_ik*g;
1504
1505                         /* dzeta_ij */
1506                         v3_scale(&temp1,&temp1,f_c_ik);
1507                         v3_scale(&temp2,&dist_ik,g*df_c_ik);
1508                         v3_add(&temp1,&temp1,&temp2);
1509                         v3_add(dzeta,dzeta,&temp1);
1510                 }
1511         }
1512
1513         /* dV_ji stuff (in b_ji there is f_c_jk) + dV_jk stuff! */
1514
1515         /* dist_jk, d_jk */
1516         v3_sub(&dist_jk,&(ak->r),&(aj->r));
1517         if(bc) check_per_bound(moldyn,&dist_jk);
1518         d_jk=v3_norm(&dist_jk);
1519
1520         /* jk constants */
1521         num=aj->bnum;
1522         if(num==ak->bnum) {
1523                 R=params->R[num];
1524                 S=params->S[num];
1525                 B=params->B[num];
1526                 mu=params->mu[num];
1527                 chi=1.0;
1528         }
1529         else {
1530                 R=params->Rmixed;
1531                 S=params->Smixed;
1532                 B=params->Bmixed;
1533                 mu=params->mu_m;
1534                 chi=params->chi;
1535         }
1536
1537         /* zeta_ji/dzeta_ji contribution only for d_jk < S_jk */
1538         if(d_jk<S) {
1539
1540                 /* constants_j from exchange data */
1541                 n=*(exchange->n_j);
1542                 c=*(exchange->c_j);
1543                 d=*(exchange->d_j);
1544                 h=*(exchange->h_j);
1545                 c2=exchange->cj2;
1546                 d2=exchange->dj2;
1547                 c2d2=exchange->cj2dj2;
1548
1549                 /* cosine of theta_jik by scalaproduct */
1550                 rr=v3_scalar_product(&dist_ij,&dist_jk); /* times -1 */
1551                 dd=d_ij*d_jk;
1552                 cos_theta=rr/dd;
1553
1554                 /* d_costheta */
1555                 d_costheta1=1.0/(d_jk*d_ij);
1556                 d_costheta2=cos_theta/(d_ij*d_ij); /* in fact -cos(), but ^ */
1557
1558                 /* some usefull values */
1559                 h_cos=(h-cos_theta);
1560                 d2_h_cos2=d2+(h_cos*h_cos);
1561                 frac=c2/(d2_h_cos2);
1562
1563                 /* g(cos_theta) */
1564                 g=1.0+c2d2-frac;
1565
1566                 /* d_costheta_ij and dg(cos_theta) - needed in any case! */
1567                 v3_scale(&temp1,&dist_jk,d_costheta1);
1568                 v3_scale(&temp2,&dist_ij,-d_costheta2); /* ji -> ij => -1 */
1569                 v3_add(&temp1,&temp1,&temp2);
1570                 v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */
1571
1572                 /* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */
1573                 dzeta=&(exchange->dzeta_ji);
1574                 if(d_jk<R) {
1575                         /* f_c_jk */
1576                         f_c_jk=1.0;
1577
1578                         /* zeta_ji */
1579                         exchange->zeta_ji+=g;
1580
1581                         /* dzeta_ji */
1582                         v3_add(dzeta,dzeta,&temp1);
1583                 }
1584                 else {
1585                         /* f_c_jk */
1586                         s_r=S-R;
1587                         arg=M_PI*(d_jk-R)/s_r;
1588                         f_c_jk=0.5+0.5*cos(arg);
1589
1590                         /* zeta_ji */
1591                         exchange->zeta_ji+=f_c_jk*g;
1592
1593                         /* dzeta_ij */
1594                         v3_scale(&temp1,&temp1,f_c_jk);
1595                         v3_add(dzeta,dzeta,&temp1);
1596                 }
1597
1598                 /* dV_jk stuff | add force contribution on atom i immediately */
1599                 if(exchange->d_ij_between_rs) {
1600                         zeta=f_c*g;
1601                         v3_scale(&temp1,&temp1,f_c);
1602                         v3_scale(&temp2,&dist_ij,df_c);
1603                         v3_add(&temp1,&temp1,&temp2);
1604                 }
1605                 else {
1606                         zeta=g;
1607                         // dzeta_jk is simply dg, which is temp1
1608                 }
1609                 /* betajnj * zeta_jk ^ nj-1 */
1610                 tmp=exchange->betajnj*pow(zeta,(n-1.0));
1611                 tmp=-chi/2.0*pow(1+tmp*zeta,-1.0/(2.0*n)-1)*tmp;
1612                 v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
1613                 v3_add(&(ai->f),&(ai->f),&temp1); /* -1 skipped in f_a calc ^ */
1614                                                   /* scaled with 0.5 ^ */
1615         }
1616
1617         return 0;
1618 }
1619
1620
1621 /*
1622  * debugging / critical check functions
1623  */
1624
1625 int moldyn_bc_check(t_moldyn *moldyn) {
1626
1627         t_atom *atom;
1628         t_3dvec *dim;
1629         int i;
1630
1631         atom=moldyn->atom;
1632         dim=&(moldyn->dim);
1633
1634         for(i=0;i<moldyn->count;i++) {
1635                 if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2)
1636                         printf("FATAL: atom %d: x: %.20f (%.20f)\n",
1637                                i,atom[i].r.x*1e10,dim->x/2*1e10);
1638                 if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
1639                         printf("FATAL: atom %d: y: %.20f (%.20f)\n",
1640                                i,atom[i].r.y*1e10,dim->y/2*1e10);
1641                 if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
1642                         printf("FATAL: atom %d: z: %.20f (%.20f)\n",
1643                                i,atom[i].r.z*1e10,dim->z/2*1e10);
1644         }
1645
1646         return 0;
1647 }
1648