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