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