potentials/albe_fast.c

   1 /*
   2  * test: albe_new.c
   3  *
   4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
   5  *
   6  */
   7
   8 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <string.h>
  12 #include <sys/types.h>
  13 #include <sys/stat.h>
  14 #include <fcntl.h>
  15 #include <unistd.h>
  16 #include <sys/time.h>
  17 #include <time.h>
  18 #include <math.h>
  19
  20 #ifdef PARALLEL
  21 #include <omp.h>
  22 #endif
  23
  24 #include "../moldyn.h"
  25 #include "../math/math.h"
  26 #include "albe.h"
  27
  28 /*
  29  * virial calculation
  30  */
  31
  32 #define albe_v_calc(a,f,d)      a->virial.xx+=f->x*d->x; \
  33                                 a->virial.yy+=f->y*d->y; \
  34                                 a->virial.zz+=f->z*d->z; \
  35                                 a->virial.xy+=f->x*d->y; \
  36                                 a->virial.xz+=f->x*d->z; \
  37                                 a->virial.yz+=f->y*d->z
  38
  39 int albe_potential_force_calc(t_moldyn *moldyn) {
  40
  41         int i,j,k,count;
  42         t_atom *itom,*jtom,*ktom;
  43         t_virial *virial;
  44         t_linkcell *lc;
  45 #ifdef STATIC_LISTS
  46         int *neighbour_i[27];
  47         int p,q;
  48         t_atom *atom;
  49 #else
  50         t_list neighbour_i[27];
  51         t_list neighbour_i2[27];
  52         t_list *this,*that;
  53 #endif
  54         u8 bc_ij,bc_ik;
  55         int dnlc;
  56
  57         // needed to work
  58         t_atom *ai;
  59
  60         // optimized
  61         t_albe_mult_params *params;
  62         t_albe_exchange *exchange;
  63         t_3dvec dist_ij;
  64         double d_ij2;
  65         double d_ij;
  66         u8 brand_i;
  67         double S2;
  68         int kcount;
  69         double zeta_ij;
  70         double pre_dzeta;
  71
  72         // more ...
  73         double Rk,Sk,Sk2;
  74         t_3dvec dist_ik;
  75         double d_ik2,d_ik;
  76         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
  77         double f_c_ik,df_c_ik;
  78
  79         t_3dvec force;
  80         double f_a,df_a,b,db,f_c,df_c;
  81         double f_r,df_r;
  82         double scale;
  83         double mu,B;
  84         double lambda,A;
  85         double r0;
  86         double S,R;
  87         double energy;
  88
  89         double dijdik_inv,fcdg,dfcg;
  90         t_3dvec dcosdrj,dcosdrk;
  91         t_3dvec tmp;
  92
  93         // even more ...
  94         double gamma_i;
  95         double c_i;
  96         double d_i;
  97         double h_i;
  98         double ci2;
  99         double di2;
 100         double ci2di2;
 101
 102         count=moldyn->count;
 103         itom=moldyn->atom;
 104         lc=&(moldyn->lc);
 105 #ifdef STATIC_LISTS
 106         atom=moldyn->atom;
 107 #endif
 108
 109         // optimized
 110         params=moldyn->pot_params;
 111         exchange=&(params->exchange);
 112
 113
 114         /* reset energy */
 115         moldyn->energy=0.0;
 116
 117         /* reset global virial */
 118         memset(&(moldyn->gvir),0,sizeof(t_virial));
 119
 120         /* reset force, site energy and virial of every atom */
 121 #ifdef PARALLEL
 122         #pragma omp parallel for private(virial)
 123 #endif
 124         for(i=0;i<count;i++) {
 125
 126                 /* reset force */
 127                 v3_zero(&(itom[i].f));
 128
 129                 /* reset virial */
 130                 virial=(&(itom[i].virial));
 131                 virial->xx=0.0;
 132                 virial->yy=0.0;
 133                 virial->zz=0.0;
 134                 virial->xy=0.0;
 135                 virial->xz=0.0;
 136                 virial->yz=0.0;
 137
 138                 /* reset site energy */
 139                 itom[i].e=0.0;
 140
 141         }
 142
 143         /* get energy, force and virial of every atom */
 144
 145         /* first (and only) loop over atoms i */
 146         for(i=0;i<count;i++) {
 147
 148                 if(!(itom[i].attr&ATOM_ATTR_3BP))
 149                         continue;
 150
 151                 link_cell_neighbour_index(moldyn,
 152                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
 153                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
 154                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
 155                                           neighbour_i);
 156
 157                 dnlc=lc->dnlc;
 158
 159                 /* copy the neighbour lists */
 160 #ifndef STATIC_LISTS
 161                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
 162 #endif
 163
 164                 ai=&(itom[i]);
 165                 brand_i=ai->brand;
 166
 167                 /* loop over atoms j */
 168                 for(j=0;j<27;j++) {
 169
 170                         bc_ij=(j<dnlc)?0:1;
 171 #ifdef STATIC_LISTS
 172                         p=0;
 173
 174                         while(neighbour_i[j][p]!=0) {
 175
 176                                 jtom=&(atom[neighbour_i[j][p]]);
 177                                 p++;
 178 #else
 179                         this=&(neighbour_i[j]);
 180                         list_reset_f(this);
 181
 182                         if(this->start==NULL)
 183                                 continue;
 184
 185                         do {
 186
 187                                 jtom=this->current->data;
 188 #endif
 189
 190                                 if(jtom==&(itom[i]))
 191                                         continue;
 192
 193                                 if(!(jtom->attr&ATOM_ATTR_3BP))
 194                                         continue;
 195
 196                                 /* reset 3bp run */
 197                                 moldyn->run3bp=1;
 198
 199
 200 /* j1 func here ... */
 201 /* albe 3 body potential function (first ij loop) */
 202
 203         /* reset zeta sum */
 204         zeta_ij=0.0;
 205
 206         /*
 207          * set ij depending values
 208          */
 209
 210         if(brand_i==jtom->brand) {
 211                 S2=params->S2[brand_i];
 212         }
 213         else {
 214                 S2=params->S2mixed;
 215         }
 216
 217         /* dist_ij, d_ij2 */
 218         v3_sub(&dist_ij,&(jtom->r),&(ai->r));
 219         if(bc_ij) check_per_bound(moldyn,&dist_ij);
 220         d_ij2=v3_absolute_square(&dist_ij);
 221
 222         /* if d_ij2 > S2 => no force & potential energy contribution */
 223         if(d_ij2>S2)
 224                 continue;
 225
 226         /* d_ij */
 227         d_ij=sqrt(d_ij2);
 228
 229         /* reset k counter for first k loop */
 230         kcount=0;
 231
 232                                 /* first loop over atoms k */
 233                                 for(k=0;k<27;k++) {
 234
 235                                         bc_ik=(k<dnlc)?0:1;
 236 #ifdef STATIC_LISTS
 237                                         q=0;
 238
 239                                         while(neighbour_i[j][q]!=0) {
 240
 241                                                 ktom=&(atom[neighbour_i[k][q]]);
 242                                                 q++;
 243 #else
 244                                         that=&(neighbour_i2[k]);
 245                                         list_reset_f(that);
 246
 247                                         if(that->start==NULL)
 248                                                 continue;
 249
 250                                         do {
 251                                                 ktom=that->current->data;
 252 #endif
 253
 254                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
 255                                                         continue;
 256
 257                                                 if(ktom==jtom)
 258                                                         continue;
 259
 260                                                 if(ktom==&(itom[i]))
 261                                                         continue;
 262
 263
 264 /* k1 func here ... */
 265 /* albe 3 body potential function (first k loop) */
 266
 267         if(kcount>ALBE_MAXN) {
 268                 printf("FATAL: neighbours = %d\n",kcount);
 269                 printf("  -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
 270         }
 271
 272         /* ik constants */
 273         if(brand_i==ktom->brand) {
 274                 Rk=params->R[brand_i];
 275                 Sk=params->S[brand_i];
 276                 Sk2=params->S2[brand_i];
 277                 /* albe needs i,k depending c,d,h and gamma values */
 278                 gamma_i=params->gamma[brand_i];
 279                 c_i=params->c[brand_i];
 280                 d_i=params->d[brand_i];
 281                 h_i=params->h[brand_i];
 282                 ci2=params->c2[brand_i];
 283                 di2=params->d2[brand_i];
 284                 ci2di2=params->c2d2[brand_i];
 285         }
 286         else {
 287                 Rk=params->Rmixed;
 288                 Sk=params->Smixed;
 289                 Sk2=params->S2mixed;
 290                 /* albe needs i,k depending c,d,h and gamma values */
 291                 gamma_i=params->gamma_m;
 292                 c_i=params->c_mixed;
 293                 d_i=params->d_mixed;
 294                 h_i=params->h_mixed;
 295                 ci2=params->c2_mixed;
 296                 di2=params->d2_mixed;
 297                 ci2di2=params->c2d2_m;
 298         }
 299
 300         /* dist_ik, d_ik2 */
 301         v3_sub(&dist_ik,&(ktom->r),&(ai->r));
 302         if(bc_ik) check_per_bound(moldyn,&dist_ik);
 303         d_ik2=v3_absolute_square(&dist_ik);
 304
 305         /* store data for second k loop */
 306         exchange->dist_ik[kcount]=dist_ik;
 307         exchange->d_ik2[kcount]=d_ik2;
 308
 309         /* return if not within cutoff */
 310         if(d_ik2>Sk2) {
 311                 kcount++;
 312                 continue;
 313         }
 314
 315         /* d_ik */
 316         d_ik=sqrt(d_ik2);
 317
 318         /* cos theta */
 319         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 320
 321         /* g_ijk
 322         h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
 323         d2_h_cos2=exchange->di2+(h_cos*h_cos);
 324         frac=exchange->ci2/d2_h_cos2;
 325         g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
 326         dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
 327         */
 328
 329         h_cos=h_i+cos_theta; // + in albe formalism
 330         d2_h_cos2=di2+(h_cos*h_cos);
 331         frac=ci2/d2_h_cos2;
 332         g=gamma_i*(1.0+ci2di2-frac);
 333         dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
 334
 335         /* zeta sum += f_c_ik * g_ijk */
 336         if(d_ik<=Rk) {
 337                 zeta_ij+=g;
 338                 f_c_ik=1.0;
 339                 df_c_ik=0.0;
 340         }
 341         else {
 342                 s_r=Sk-Rk;
 343                 arg=M_PI*(d_ik-Rk)/s_r;
 344                 f_c_ik=0.5+0.5*cos(arg);
 345                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 346                 zeta_ij+=f_c_ik*g;
 347         }
 348
 349         /* store even more data for second k loop */
 350         exchange->g[kcount]=g;
 351         exchange->dg[kcount]=dg;
 352         exchange->d_ik[kcount]=d_ik;
 353         exchange->cos_theta[kcount]=cos_theta;
 354         exchange->f_c_ik[kcount]=f_c_ik;
 355         exchange->df_c_ik[kcount]=df_c_ik;
 356
 357         /* increase k counter */
 358         kcount++;
 359
 360 #ifdef STATIC_LISTS
 361                                         }
 362 #else
 363                                         } while(list_next_f(that)!=\
 364                                                 L_NO_NEXT_ELEMENT);
 365 #endif
 366
 367                                 }
 368
 369
 370 /* j2 func here ... */
 371
 372
 373         if(brand_i==jtom->brand) {
 374                 S=params->S[brand_i];
 375                 R=params->R[brand_i];
 376                 B=params->B[brand_i];
 377                 A=params->A[brand_i];
 378                 r0=params->r0[brand_i];
 379                 mu=params->mu[brand_i];
 380                 lambda=params->lambda[brand_i];
 381         }
 382         else {
 383                 S=params->Smixed;
 384                 R=params->Rmixed;
 385                 B=params->Bmixed;
 386                 A=params->Amixed;
 387                 r0=params->r0_mixed;
 388                 mu=params->mu_m;
 389                 lambda=params->lambda_m;
 390         }
 391
 392         /* f_c, df_c */
 393         if(d_ij<R) {
 394                 f_c=1.0;
 395                 df_c=0.0;
 396         }
 397         else {
 398                 s_r=S-R;
 399                 arg=M_PI*(d_ij-R)/s_r;
 400                 f_c=0.5+0.5*cos(arg);
 401                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 402         }
 403
 404         /* f_a, df_a */
 405         f_a=-B*exp(-mu*(d_ij-r0));
 406         df_a=mu*f_a/d_ij;
 407
 408         /* f_r, df_r */
 409         f_r=A*exp(-lambda*(d_ij-r0));
 410         df_r=lambda*f_r/d_ij;
 411
 412         /* b, db */
 413         if(zeta_ij==0.0) {
 414                 b=1.0;
 415                 db=0.0;
 416         }
 417         else {
 418                 b=1.0/sqrt(1.0+zeta_ij);
 419                 db=-0.5*b/(1.0+zeta_ij);
 420         }
 421
 422         /* force contribution for atom i */
 423         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
 424         v3_scale(&force,&(dist_ij),scale);
 425         v3_add(&(ai->f),&(ai->f),&force);
 426
 427         /* force contribution for atom j */
 428         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 429         v3_add(&(jtom->f),&(jtom->f),&force);
 430
 431         /* virial */
 432         virial_calc(ai,&force,&(dist_ij));
 433
 434 #ifdef DEBUG
 435 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 436         if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
 437                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
 438                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 439                 if(ai==&(moldyn->atom[0]))
 440                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 441                 if(jtom==&(moldyn->atom[0]))
 442                         printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
 443                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 444                                                     f_c,b,f_a,f_r);
 445                 printf("          %f %f %f\n",zeta_ij,.0,.0);
 446         }
 447 }
 448 #endif
 449
 450         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
 451         pre_dzeta=0.5*f_a*f_c*db;
 452
 453         /* energy contribution */
 454         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
 455         moldyn->energy+=energy;
 456         ai->e+=energy;
 457
 458         /* reset k counter for second k loop */
 459         kcount=0;
 460
 461
 462                                 /* second loop over atoms k */
 463                                 for(k=0;k<27;k++) {
 464
 465                                         bc_ik=(k<dnlc)?0:1;
 466 #ifdef STATIC_LISTS
 467                                         q=0;
 468
 469                                         while(neighbour_i[j][q]!=0) {
 470
 471                                                 ktom=&(atom[neighbour_i[k][q]]);
 472                                                 q++;
 473 #else
 474                                         that=&(neighbour_i2[k]);
 475                                         list_reset_f(that);
 476
 477                                         if(that->start==NULL)
 478                                                 continue;
 479
 480                                         do {
 481                                                 ktom=that->current->data;
 482 #endif
 483
 484                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
 485                                                         continue;
 486
 487                                                 if(ktom==jtom)
 488                                                         continue;
 489
 490                                                 if(ktom==&(itom[i]))
 491                                                         continue;
 492
 493
 494 /* k2 func here ... */
 495 /* albe 3 body potential function (second k loop) */
 496
 497         if(kcount>ALBE_MAXN)
 498                 printf("FATAL: neighbours!\n");
 499
 500         /* d_ik2 */
 501         d_ik2=exchange->d_ik2[kcount];
 502
 503         if(brand_i==ktom->brand)
 504                 Sk2=params->S2[brand_i];
 505         else
 506                 Sk2=params->S2mixed;
 507
 508         /* return if d_ik > S */
 509         if(d_ik2>Sk2) {
 510                 kcount++;
 511                 continue;
 512         }
 513
 514         /* dist_ik, d_ik */
 515         dist_ik=exchange->dist_ik[kcount];
 516         d_ik=exchange->d_ik[kcount];
 517
 518         /* f_c_ik, df_c_ik */
 519         f_c_ik=exchange->f_c_ik[kcount];
 520         df_c_ik=exchange->df_c_ik[kcount];
 521
 522         /* g, dg, cos_theta */
 523         g=exchange->g[kcount];
 524         dg=exchange->dg[kcount];
 525         cos_theta=exchange->cos_theta[kcount];
 526
 527         /* cos_theta derivatives wrt j,k */
 528         dijdik_inv=1.0/(d_ij*d_ik);
 529         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
 530         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
 531         v3_add(&dcosdrj,&dcosdrj,&tmp);
 532         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
 533         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 534         v3_add(&dcosdrk,&dcosdrk,&tmp);
 535
 536         /* f_c_ik * dg, df_c_ik * g */
 537         fcdg=f_c_ik*dg;
 538         dfcg=df_c_ik*g;
 539
 540         /* derivative wrt j */
 541         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 542
 543         /* force contribution */
 544         v3_add(&(jtom->f),&(jtom->f),&force);
 545
 546 #ifdef DEBUG
 547 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 548         if(jtom==&(moldyn->atom[DATOM])) {
 549                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
 550                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 551                 printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
 552                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 553                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 554         }
 555 }
 556 #endif
 557
 558         /* virial */
 559         virial_calc(ai,&force,&dist_ij);
 560
 561         /* force contribution to atom i */
 562         v3_scale(&force,&force,-1.0);
 563         v3_add(&(ai->f),&(ai->f),&force);
 564
 565         /* derivative wrt k */
 566         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
 567         v3_scale(&tmp,&dcosdrk,fcdg);
 568         v3_add(&force,&force,&tmp);
 569         v3_scale(&force,&force,pre_dzeta);
 570
 571         /* force contribution */
 572         v3_add(&(ktom->f),&(ktom->f),&force);
 573
 574 #ifdef DEBUG
 575 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 576         if(ktom==&(moldyn->atom[DATOM])) {
 577                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
 578                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 579                 printf("  total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
 580                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 581                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 582         }
 583 }
 584 #endif
 585
 586         /* virial */
 587         virial_calc(ai,&force,&dist_ik);
 588
 589         /* force contribution to atom i */
 590         v3_scale(&force,&force,-1.0);
 591         v3_add(&(ai->f),&(ai->f),&force);
 592
 593         /* increase k counter */
 594         kcount++;
 595
 596
 597
 598 #ifdef STATIC_LISTS
 599                                         }
 600 #else
 601                                         } while(list_next_f(that)!=\
 602                                                 L_NO_NEXT_ELEMENT);
 603 #endif
 604
 605                                 }
 606
 607 #ifdef STATIC_LISTS
 608                         }
 609 #else
 610                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 611 #endif
 612
 613                 }
 614
 615 #ifdef DEBUG
 616         //printf("\n\n");
 617 #endif
 618 #ifdef VDEBUG
 619         printf("\n\n");
 620 #endif
 621
 622         }
 623
 624 #ifdef DEBUG
 625         //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
 626         if(moldyn->time>DSTART&&moldyn->time<DEND) {
 627                 printf("force:\n");
 628                 printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
 629                 printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
 630                 printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
 631         }
 632 #endif
 633
 634         /* some postprocessing */
 635 #ifdef PARALLEL
 636         #pragma omp parallel for
 637 #endif
 638         for(i=0;i<count;i++) {
 639                 /* calculate global virial */
 640                 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
 641                 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
 642                 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
 643                 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
 644                 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
 645                 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
 646
 647                 /* check forces regarding the given timestep */
 648                 if(v3_norm(&(itom[i].f))>\
 649                    0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
 650                         printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
 651                                i);
 652         }
 653
 654         return 0;
 655 }
 656