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