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 PTHREADS
  25 #include <pthread.h>
  26 #define MAX_THREADS 2
  27 #endif
  28
  29 #include "../moldyn.h"
  30 #include "../math/math.h"
  31 #include "albe.h"
  32
  33 #ifdef PTHREADS
  34 extern pthread_mutex_t *amutex;
  35 extern pthread_mutex_t emutex;
  36 #endif
  37
  38 /*
  39  * virial calculation
  40  */
  41
  42 #define albe_v_calc(a,f,d)      (a)->virial.xx+=(f)->x*(d)->x; \
  43                                 (a)->virial.yy+=(f)->y*(d)->y; \
  44                                 (a)->virial.zz+=(f)->z*(d)->z; \
  45                                 (a)->virial.xy+=(f)->x*(d)->y; \
  46                                 (a)->virial.xz+=(f)->x*(d)->z; \
  47                                 (a)->virial.yz+=(f)->y*(d)->z
  48
  49 #ifndef PTHREADS
  50
  51 int albe_potential_force_calc(t_moldyn *moldyn) {
  52
  53         int i,j,k,count;
  54         t_atom *itom,*jtom,*ktom;
  55         t_virial *virial;
  56         t_linkcell *lc;
  57 #ifdef STATIC_LISTS
  58         int *neighbour_i[27];
  59         int p,q;
  60 #elif LOWMEM_LISTS
  61         int neighbour_i[27];
  62         int p,q;
  63 #else
  64         t_list neighbour_i[27];
  65         t_list neighbour_i2[27];
  66         t_list *this,*that;
  67 #endif
  68         u8 bc_ij,bc_ik;
  69         int dnlc;
  70 #ifdef PTHREADS
  71         int ret;
  72         t_kdata kdata[27];
  73         pthread_t kthread[27];
  74 #endif
  75
  76         // needed to work
  77         t_atom *ai;
  78
  79         // optimized
  80         t_albe_mult_params *params;
  81         t_albe_exchange *exchange;
  82         t_3dvec dist_ij;
  83         double d_ij2;
  84         double d_ij;
  85         u8 brand_i;
  86         double S2;
  87         int kcount;
  88         double zeta_ij;
  89         double pre_dzeta;
  90
  91         // more ...
  92         double Rk,Sk,Sk2;
  93         t_3dvec dist_ik;
  94         double d_ik2,d_ik;
  95         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
  96         double f_c_ik,df_c_ik;
  97
  98         t_3dvec force;
  99         double f_a,df_a,b,db,f_c,df_c;
 100         double f_r,df_r;
 101         double scale;
 102         double mu,B;
 103         double lambda,A;
 104         double r0;
 105         double S,R;
 106         double energy;
 107
 108         double dijdik_inv,fcdg,dfcg;
 109         t_3dvec dcosdrj,dcosdrk;
 110         t_3dvec tmp;
 111
 112         // even more ...
 113         double gamma_i;
 114         double c_i;
 115         double d_i;
 116         double h_i;
 117         double ci2;
 118         double di2;
 119         double ci2di2;
 120
 121         count=moldyn->count;
 122         itom=moldyn->atom;
 123         lc=&(moldyn->lc);
 124
 125         // optimized
 126         params=moldyn->pot_params;
 127         exchange=&(params->exchange);
 128
 129
 130         /* reset energy */
 131         moldyn->energy=0.0;
 132
 133         /* reset global virial */
 134         memset(&(moldyn->gvir),0,sizeof(t_virial));
 135
 136         /* reset force, site energy and virial of every atom */
 137 #ifdef PARALLEL
 138         #pragma omp parallel for private(virial)
 139 #endif
 140         for(i=0;i<count;i++) {
 141
 142                 /* reset force */
 143                 v3_zero(&(itom[i].f));
 144
 145                 /* reset virial */
 146                 virial=(&(itom[i].virial));
 147                 virial->xx=0.0;
 148                 virial->yy=0.0;
 149                 virial->zz=0.0;
 150                 virial->xy=0.0;
 151                 virial->xz=0.0;
 152                 virial->yz=0.0;
 153
 154                 /* reset site energy */
 155                 itom[i].e=0.0;
 156
 157         }
 158
 159         /* get energy, force and virial of every atom */
 160
 161         /* first (and only) loop over atoms i */
 162         for(i=0;i<count;i++) {
 163
 164                 if(!(itom[i].attr&ATOM_ATTR_3BP))
 165                         continue;
 166
 167                 link_cell_neighbour_index(moldyn,
 168                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
 169                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
 170                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
 171                                           neighbour_i);
 172
 173                 dnlc=lc->dnlc;
 174
 175                 /* copy the neighbour lists */
 176 #ifdef STATIC_LISTS
 177 #elif LOWMEM_LISTS
 178 #else
 179                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
 180 #endif
 181
 182                 ai=&(itom[i]);
 183                 brand_i=ai->brand;
 184
 185                 /* loop over atoms j */
 186                 for(j=0;j<27;j++) {
 187
 188                         bc_ij=(j<dnlc)?0:1;
 189 #ifdef STATIC_LISTS
 190                         p=0;
 191
 192                         while(neighbour_i[j][p]!=-1) {
 193
 194                                 jtom=&(itom[neighbour_i[j][p]]);
 195                                 p++;
 196 #elif LOWMEM_LISTS
 197                         p=neighbour_i[j];
 198
 199                         while(p!=-1) {
 200
 201                                 jtom=&(itom[p]);
 202                                 p=lc->subcell->list[p];
 203 #else
 204                         this=&(neighbour_i[j]);
 205                         list_reset_f(this);
 206
 207                         if(this->start==NULL)
 208                                 continue;
 209
 210                         do {
 211
 212                                 jtom=this->current->data;
 213 #endif
 214
 215                                 if(jtom==&(itom[i]))
 216                                         continue;
 217
 218                                 if(!(jtom->attr&ATOM_ATTR_3BP))
 219                                         continue;
 220
 221                                 /* reset 3bp run */
 222                                 moldyn->run3bp=1;
 223
 224
 225 /* j1 func here ... */
 226 /* albe 3 body potential function (first ij loop) */
 227
 228         /* reset zeta sum */
 229         zeta_ij=0.0;
 230
 231         /*
 232          * set ij depending values
 233          */
 234
 235         if(brand_i==jtom->brand) {
 236                 S2=params->S2[brand_i];
 237         }
 238         else {
 239                 S2=params->S2mixed;
 240         }
 241
 242         /* dist_ij, d_ij2 */
 243         v3_sub(&dist_ij,&(jtom->r),&(ai->r));
 244         if(bc_ij) check_per_bound(moldyn,&dist_ij);
 245         d_ij2=v3_absolute_square(&dist_ij);
 246
 247         /* if d_ij2 > S2 => no force & potential energy contribution */
 248         if(d_ij2>S2)
 249                 continue;
 250
 251         /* d_ij */
 252         d_ij=sqrt(d_ij2);
 253
 254         /* reset k counter for first k loop */
 255         kcount=0;
 256
 257                                 /* first loop over atoms k */
 258                                 for(k=0;k<27;k++) {
 259
 260                                         bc_ik=(k<dnlc)?0:1;
 261 #ifdef STATIC_LISTS
 262                                         q=0;
 263
 264                                         while(neighbour_i[k][q]!=-1) {
 265
 266                                                 ktom=&(itom[neighbour_i[k][q]]);
 267                                                 q++;
 268 #elif LOWMEM_LISTS
 269                                         q=neighbour_i[k];
 270
 271                                         while(q!=-1) {
 272
 273                                                 ktom=&(itom[q]);
 274                                                 q=lc->subcell->list[q];
 275 #else
 276                                         that=&(neighbour_i2[k]);
 277                                         list_reset_f(that);
 278
 279                                         if(that->start==NULL)
 280                                                 continue;
 281
 282                                         do {
 283                                                 ktom=that->current->data;
 284 #endif
 285
 286                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
 287                                                         continue;
 288
 289                                                 if(ktom==jtom)
 290                                                         continue;
 291
 292                                                 if(ktom==&(itom[i]))
 293                                                         continue;
 294
 295 /* k1 func here ... */
 296 /* albe 3 body potential function (first k loop) */
 297
 298         if(kcount>ALBE_MAXN) {
 299                 printf("FATAL: neighbours = %d\n",kcount);
 300                 printf("  -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
 301         }
 302
 303         /* ik constants */
 304         if(brand_i==ktom->brand) {
 305                 Rk=params->R[brand_i];
 306                 Sk=params->S[brand_i];
 307                 Sk2=params->S2[brand_i];
 308                 /* albe needs i,k depending c,d,h and gamma values */
 309                 gamma_i=params->gamma[brand_i];
 310                 c_i=params->c[brand_i];
 311                 d_i=params->d[brand_i];
 312                 h_i=params->h[brand_i];
 313                 ci2=params->c2[brand_i];
 314                 di2=params->d2[brand_i];
 315                 ci2di2=params->c2d2[brand_i];
 316         }
 317         else {
 318                 Rk=params->Rmixed;
 319                 Sk=params->Smixed;
 320                 Sk2=params->S2mixed;
 321                 /* albe needs i,k depending c,d,h and gamma values */
 322                 gamma_i=params->gamma_m;
 323                 c_i=params->c_mixed;
 324                 d_i=params->d_mixed;
 325                 h_i=params->h_mixed;
 326                 ci2=params->c2_mixed;
 327                 di2=params->d2_mixed;
 328                 ci2di2=params->c2d2_m;
 329         }
 330
 331         /* dist_ik, d_ik2 */
 332         v3_sub(&dist_ik,&(ktom->r),&(ai->r));
 333         if(bc_ik) check_per_bound(moldyn,&dist_ik);
 334         d_ik2=v3_absolute_square(&dist_ik);
 335
 336         /* store data for second k loop */
 337         exchange->dist_ik[kcount]=dist_ik;
 338         exchange->d_ik2[kcount]=d_ik2;
 339
 340         /* return if not within cutoff */
 341         if(d_ik2>Sk2) {
 342                 kcount++;
 343                 continue;
 344         }
 345
 346         /* d_ik */
 347         d_ik=sqrt(d_ik2);
 348
 349         /* cos theta */
 350         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 351
 352         /* g_ijk
 353         h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
 354         d2_h_cos2=exchange->di2+(h_cos*h_cos);
 355         frac=exchange->ci2/d2_h_cos2;
 356         g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
 357         dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
 358         */
 359
 360         h_cos=h_i+cos_theta; // + in albe formalism
 361         d2_h_cos2=di2+(h_cos*h_cos);
 362         frac=ci2/d2_h_cos2;
 363         g=gamma_i*(1.0+ci2di2-frac);
 364         dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
 365
 366         /* zeta sum += f_c_ik * g_ijk */
 367         if(d_ik<=Rk) {
 368                 zeta_ij+=g;
 369                 f_c_ik=1.0;
 370                 df_c_ik=0.0;
 371         }
 372         else {
 373                 s_r=Sk-Rk;
 374                 arg=M_PI*(d_ik-Rk)/s_r;
 375                 f_c_ik=0.5+0.5*cos(arg);
 376                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
 377                 zeta_ij+=f_c_ik*g;
 378         }
 379
 380         /* store even more data for second k loop */
 381         exchange->g[kcount]=g;
 382         exchange->dg[kcount]=dg;
 383         exchange->d_ik[kcount]=d_ik;
 384         exchange->cos_theta[kcount]=cos_theta;
 385         exchange->f_c_ik[kcount]=f_c_ik;
 386         exchange->df_c_ik[kcount]=df_c_ik;
 387
 388         /* increase k counter */
 389         kcount++;
 390
 391 #ifdef STATIC_LISTS
 392                                         }
 393 #elif LOWMEM_LISTS
 394                                         }
 395 #else
 396                                         } while(list_next_f(that)!=\
 397                                                 L_NO_NEXT_ELEMENT);
 398 #endif
 399
 400                                 }
 401
 402 /* j2 func here ... */
 403
 404
 405         if(brand_i==jtom->brand) {
 406                 S=params->S[brand_i];
 407                 R=params->R[brand_i];
 408                 B=params->B[brand_i];
 409                 A=params->A[brand_i];
 410                 r0=params->r0[brand_i];
 411                 mu=params->mu[brand_i];
 412                 lambda=params->lambda[brand_i];
 413         }
 414         else {
 415                 S=params->Smixed;
 416                 R=params->Rmixed;
 417                 B=params->Bmixed;
 418                 A=params->Amixed;
 419                 r0=params->r0_mixed;
 420                 mu=params->mu_m;
 421                 lambda=params->lambda_m;
 422         }
 423
 424         /* f_c, df_c */
 425         if(d_ij<R) {
 426                 f_c=1.0;
 427                 df_c=0.0;
 428         }
 429         else {
 430                 s_r=S-R;
 431                 arg=M_PI*(d_ij-R)/s_r;
 432                 f_c=0.5+0.5*cos(arg);
 433                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
 434         }
 435
 436         /* f_a, df_a */
 437         f_a=-B*exp(-mu*(d_ij-r0));
 438         df_a=mu*f_a/d_ij;
 439
 440         /* f_r, df_r */
 441         f_r=A*exp(-lambda*(d_ij-r0));
 442         df_r=lambda*f_r/d_ij;
 443
 444         /* b, db */
 445         if(zeta_ij==0.0) {
 446                 b=1.0;
 447                 db=0.0;
 448         }
 449         else {
 450                 b=1.0/sqrt(1.0+zeta_ij);
 451                 db=-0.5*b/(1.0+zeta_ij);
 452         }
 453
 454         /* force contribution for atom i */
 455         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
 456         v3_scale(&force,&(dist_ij),scale);
 457         v3_add(&(ai->f),&(ai->f),&force);
 458
 459         /* force contribution for atom j */
 460         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 461         v3_add(&(jtom->f),&(jtom->f),&force);
 462
 463         /* virial */
 464         albe_v_calc(ai,&force,&(dist_ij));
 465         //virial_calc(ai,&force,&(dist_ij));
 466
 467 #ifdef DEBUG
 468 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 469         if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
 470                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
 471                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 472                 if(ai==&(moldyn->atom[0]))
 473                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 474                 if(jtom==&(moldyn->atom[0]))
 475                         printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
 476                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 477                                                     f_c,b,f_a,f_r);
 478                 printf("          %f %f %f\n",zeta_ij,.0,.0);
 479         }
 480 }
 481 #endif
 482
 483         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
 484         pre_dzeta=0.5*f_a*f_c*db;
 485
 486         /* energy contribution */
 487         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
 488         moldyn->energy+=energy;
 489         ai->e+=energy;
 490
 491         /* reset k counter for second k loop */
 492         kcount=0;
 493
 494
 495                                 /* second loop over atoms k */
 496                                 for(k=0;k<27;k++) {
 497
 498                                         bc_ik=(k<dnlc)?0:1;
 499 #ifdef STATIC_LISTS
 500                                         q=0;
 501
 502                                         while(neighbour_i[k][q]!=-1) {
 503
 504                                                 ktom=&(itom[neighbour_i[k][q]]);
 505                                                 q++;
 506 #elif LOWMEM_LISTS
 507                                         q=neighbour_i[k];
 508
 509                                         while(q!=-1) {
 510
 511                                                 ktom=&(itom[q]);
 512                                                 q=lc->subcell->list[q];
 513 #else
 514                                         that=&(neighbour_i2[k]);
 515                                         list_reset_f(that);
 516
 517                                         if(that->start==NULL)
 518                                                 continue;
 519
 520                                         do {
 521                                                 ktom=that->current->data;
 522 #endif
 523
 524                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
 525                                                         continue;
 526
 527                                                 if(ktom==jtom)
 528                                                         continue;
 529
 530                                                 if(ktom==&(itom[i]))
 531                                                         continue;
 532
 533
 534 /* k2 func here ... */
 535 /* albe 3 body potential function (second k loop) */
 536
 537         if(kcount>ALBE_MAXN)
 538                 printf("FATAL: neighbours!\n");
 539
 540         /* d_ik2 */
 541         d_ik2=exchange->d_ik2[kcount];
 542
 543         if(brand_i==ktom->brand)
 544                 Sk2=params->S2[brand_i];
 545         else
 546                 Sk2=params->S2mixed;
 547
 548         /* return if d_ik > S */
 549         if(d_ik2>Sk2) {
 550                 kcount++;
 551                 continue;
 552         }
 553
 554         /* dist_ik, d_ik */
 555         dist_ik=exchange->dist_ik[kcount];
 556         d_ik=exchange->d_ik[kcount];
 557
 558         /* f_c_ik, df_c_ik */
 559         f_c_ik=exchange->f_c_ik[kcount];
 560         df_c_ik=exchange->df_c_ik[kcount];
 561
 562         /* g, dg, cos_theta */
 563         g=exchange->g[kcount];
 564         dg=exchange->dg[kcount];
 565         cos_theta=exchange->cos_theta[kcount];
 566
 567         /* cos_theta derivatives wrt j,k */
 568         dijdik_inv=1.0/(d_ij*d_ik);
 569         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
 570         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
 571         v3_add(&dcosdrj,&dcosdrj,&tmp);
 572         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
 573         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 574         v3_add(&dcosdrk,&dcosdrk,&tmp);
 575
 576         /* f_c_ik * dg, df_c_ik * g */
 577         fcdg=f_c_ik*dg;
 578         dfcg=df_c_ik*g;
 579
 580         /* derivative wrt j */
 581         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 582
 583         /* force contribution */
 584         v3_add(&(jtom->f),&(jtom->f),&force);
 585
 586 #ifdef DEBUG
 587 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 588         if(jtom==&(moldyn->atom[DATOM])) {
 589                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
 590                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 591                 printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
 592                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 593                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 594         }
 595 }
 596 #endif
 597
 598         /* virial */
 599         albe_v_calc(ai,&force,&dist_ij);
 600         //virial_calc(ai,&force,&dist_ij);
 601
 602         /* force contribution to atom i */
 603         v3_scale(&force,&force,-1.0);
 604         v3_add(&(ai->f),&(ai->f),&force);
 605
 606         /* derivative wrt k */
 607         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
 608         v3_scale(&tmp,&dcosdrk,fcdg);
 609         v3_add(&force,&force,&tmp);
 610         v3_scale(&force,&force,pre_dzeta);
 611
 612         /* force contribution */
 613         v3_add(&(ktom->f),&(ktom->f),&force);
 614
 615 #ifdef DEBUG
 616 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 617         if(ktom==&(moldyn->atom[DATOM])) {
 618                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
 619                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 620                 printf("  total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
 621                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 622                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 623         }
 624 }
 625 #endif
 626
 627         /* virial */
 628         albe_v_calc(ai,&force,&dist_ik);
 629         //virial_calc(ai,&force,&dist_ik);
 630
 631         /* force contribution to atom i */
 632         v3_scale(&force,&force,-1.0);
 633         v3_add(&(ai->f),&(ai->f),&force);
 634
 635         /* increase k counter */
 636         kcount++;
 637
 638
 639
 640 #ifdef STATIC_LISTS
 641                                         }
 642 #elif LOWMEM_LISTS
 643                                         }
 644 #else
 645                                         } while(list_next_f(that)!=\
 646                                                 L_NO_NEXT_ELEMENT);
 647 #endif
 648
 649                                 }
 650
 651 #ifdef STATIC_LISTS
 652                         }
 653 #elif LOWMEM_LISTS
 654                         }
 655 #else
 656                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 657 #endif
 658
 659                 }
 660
 661 #ifdef DEBUG
 662         //printf("\n\n");
 663 #endif
 664 #ifdef VDEBUG
 665         printf("\n\n");
 666 #endif
 667
 668         }
 669
 670 #ifdef DEBUG
 671         //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
 672         if(moldyn->time>DSTART&&moldyn->time<DEND) {
 673                 printf("force:\n");
 674                 printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
 675                 printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
 676                 printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
 677         }
 678 #endif
 679
 680         /* some postprocessing */
 681 #ifdef PARALLEL
 682         #pragma omp parallel for
 683 #endif
 684         for(i=0;i<count;i++) {
 685                 /* calculate global virial */
 686                 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
 687                 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
 688                 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
 689                 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
 690                 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
 691                 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
 692
 693                 /* check forces regarding the given timestep */
 694                 if(v3_norm(&(itom[i].f))>\
 695                    0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
 696                         printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
 697                                i);
 698         }
 699
 700         return 0;
 701 }
 702
 703
 704 #else // PTHREADS
 705
 706
 707 typedef struct s_pft_data {
 708         t_moldyn *moldyn;
 709         int start,end;
 710 } t_pft_data;
 711
 712 void *potential_force_thread(void *ptr) {
 713
 714         t_pft_data *pft_data;
 715         t_moldyn *moldyn;
 716         t_albe_exchange ec;
 717
 718         int i,j,k,count;
 719         t_atom *itom,*jtom,*ktom;
 720         t_linkcell *lc;
 721 #ifdef STATIC_LISTS
 722         int *neighbour_i[27];
 723         int p,q;
 724 #elif LOWMEM_LISTS
 725         int neighbour_i[27];
 726         int p,q;
 727 #else
 728         t_list neighbour_i[27];
 729         t_list neighbour_i2[27];
 730         t_list *this,*that;
 731 #endif
 732         u8 bc_ij,bc_ik;
 733         int dnlc;
 734
 735         // needed to work
 736         t_atom *ai;
 737
 738         // optimized
 739         t_albe_mult_params *params;
 740         t_albe_exchange *exchange;
 741         t_3dvec dist_ij;
 742         double d_ij2;
 743         double d_ij;
 744         u8 brand_i;
 745         double S2;
 746         int kcount;
 747         double zeta_ij;
 748         double pre_dzeta;
 749
 750         // more ...
 751         double Rk,Sk,Sk2;
 752         t_3dvec dist_ik;
 753         double d_ik2,d_ik;
 754         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
 755         double f_c_ik,df_c_ik;
 756
 757         t_3dvec force;
 758         double f_a,df_a,b,db,f_c,df_c;
 759         double f_r,df_r;
 760         double scale;
 761         double mu,B;
 762         double lambda,A;
 763         double r0;
 764         double S,R;
 765         double energy;
 766
 767         double dijdik_inv,fcdg,dfcg;
 768         t_3dvec dcosdrj,dcosdrk;
 769         t_3dvec tmp;
 770
 771         // even more ...
 772         double gamma_i;
 773         double c_i;
 774         double d_i;
 775         double h_i;
 776         double ci2;
 777         double di2;
 778         double ci2di2;
 779
 780         pft_data=ptr;
 781         moldyn=pft_data->moldyn;
 782         exchange=&ec;
 783
 784         count=moldyn->count;
 785         itom=moldyn->atom;
 786         lc=&(moldyn->lc);
 787
 788         // optimized
 789         params=moldyn->pot_params;
 790
 791         /* get energy, force and virial for atoms */
 792
 793         for(i=pft_data->start;i<pft_data->end;i++) {
 794
 795                 if(!(itom[i].attr&ATOM_ATTR_3BP))
 796                         return 0;
 797
 798                 link_cell_neighbour_index(moldyn,
 799                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
 800                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
 801                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
 802                                           neighbour_i);
 803
 804                 dnlc=lc->dnlc;
 805
 806                 /* copy the neighbour lists */
 807 #ifdef STATIC_LISTS
 808 #elif LOWMEM_LISTS
 809 #else
 810                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
 811 #endif
 812
 813                 ai=&(itom[i]);
 814                 brand_i=ai->brand;
 815
 816                 /* loop over atoms j */
 817                 for(j=0;j<27;j++) {
 818
 819                         bc_ij=(j<dnlc)?0:1;
 820 #ifdef STATIC_LISTS
 821                         p=0;
 822
 823                         while(neighbour_i[j][p]!=-1) {
 824
 825                                 jtom=&(itom[neighbour_i[j][p]]);
 826                                 p++;
 827 #elif LOWMEM_LISTS
 828                         p=neighbour_i[j];
 829
 830                         while(p!=-1) {
 831
 832                                 jtom=&(itom[p]);
 833                                 p=lc->subcell->list[p];
 834 #else
 835                         this=&(neighbour_i[j]);
 836                         list_reset_f(this);
 837
 838                         if(this->start==NULL)
 839                                 continue;
 840
 841                         do {
 842
 843                                 jtom=this->current->data;
 844 #endif
 845
 846                                 if(jtom==&(itom[i]))
 847                                         continue;
 848
 849                                 if(!(jtom->attr&ATOM_ATTR_3BP))
 850                                         continue;
 851
 852                                 /* reset 3bp run */
 853                                 moldyn->run3bp=1;
 854
 855
 856 /* j1 func here ... */
 857 /* albe 3 body potential function (first ij loop) */
 858
 859         /* reset zeta sum */
 860         zeta_ij=0.0;
 861
 862         /*
 863          * set ij depending values
 864          */
 865
 866         if(brand_i==jtom->brand) {
 867                 S2=params->S2[brand_i];
 868         }
 869         else {
 870                 S2=params->S2mixed;
 871         }
 872
 873         /* dist_ij, d_ij2 */
 874         v3_sub(&dist_ij,&(jtom->r),&(ai->r));
 875         if(bc_ij) check_per_bound(moldyn,&dist_ij);
 876         d_ij2=v3_absolute_square(&dist_ij);
 877
 878         /* if d_ij2 > S2 => no force & potential energy contribution */
 879         if(d_ij2>S2)
 880                 continue;
 881
 882         /* d_ij */
 883         d_ij=sqrt(d_ij2);
 884
 885         /* reset k counter for first k loop */
 886         kcount=0;
 887
 888                                 /* first loop over atoms k */
 889                                 for(k=0;k<27;k++) {
 890
 891                                         bc_ik=(k<dnlc)?0:1;
 892 #ifdef STATIC_LISTS
 893                                         q=0;
 894
 895                                         while(neighbour_i[k][q]!=-1) {
 896
 897                                                 ktom=&(itom[neighbour_i[k][q]]);
 898                                                 q++;
 899 #elif LOWMEM_LISTS
 900                                         q=neighbour_i[k];
 901
 902                                         while(q!=-1) {
 903
 904                                                 ktom=&(itom[q]);
 905                                                 q=lc->subcell->list[q];
 906 #else
 907                                         that=&(neighbour_i2[k]);
 908                                         list_reset_f(that);
 909
 910                                         if(that->start==NULL)
 911                                                 continue;
 912
 913                                         do {
 914                                                 ktom=that->current->data;
 915 #endif
 916
 917                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
 918                                                         continue;
 919
 920                                                 if(ktom==jtom)
 921                                                         continue;
 922
 923                                                 if(ktom==&(itom[i]))
 924                                                         continue;
 925
 926 /* k1 func here ... */
 927 /* albe 3 body potential function (first k loop) */
 928
 929         if(kcount>ALBE_MAXN) {
 930                 printf("FATAL: neighbours = %d\n",kcount);
 931                 printf("  -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
 932         }
 933
 934         /* ik constants */
 935         if(brand_i==ktom->brand) {
 936                 Rk=params->R[brand_i];
 937                 Sk=params->S[brand_i];
 938                 Sk2=params->S2[brand_i];
 939                 /* albe needs i,k depending c,d,h and gamma values */
 940                 gamma_i=params->gamma[brand_i];
 941                 c_i=params->c[brand_i];
 942                 d_i=params->d[brand_i];
 943                 h_i=params->h[brand_i];
 944                 ci2=params->c2[brand_i];
 945                 di2=params->d2[brand_i];
 946                 ci2di2=params->c2d2[brand_i];
 947         }
 948         else {
 949                 Rk=params->Rmixed;
 950                 Sk=params->Smixed;
 951                 Sk2=params->S2mixed;
 952                 /* albe needs i,k depending c,d,h and gamma values */
 953                 gamma_i=params->gamma_m;
 954                 c_i=params->c_mixed;
 955                 d_i=params->d_mixed;
 956                 h_i=params->h_mixed;
 957                 ci2=params->c2_mixed;
 958                 di2=params->d2_mixed;
 959                 ci2di2=params->c2d2_m;
 960         }
 961
 962         /* dist_ik, d_ik2 */
 963         v3_sub(&dist_ik,&(ktom->r),&(ai->r));
 964         if(bc_ik) check_per_bound(moldyn,&dist_ik);
 965         d_ik2=v3_absolute_square(&dist_ik);
 966
 967         /* store data for second k loop */
 968         exchange->dist_ik[kcount]=dist_ik;
 969         exchange->d_ik2[kcount]=d_ik2;
 970
 971         /* return if not within cutoff */
 972         if(d_ik2>Sk2) {
 973                 kcount++;
 974                 continue;
 975         }
 976
 977         /* d_ik */
 978         d_ik=sqrt(d_ik2);
 979
 980         /* cos theta */
 981         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 982
 983         /* g_ijk
 984         h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
 985         d2_h_cos2=exchange->di2+(h_cos*h_cos);
 986         frac=exchange->ci2/d2_h_cos2;
 987         g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
 988         dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
 989         */
 990
 991         h_cos=h_i+cos_theta; // + in albe formalism
 992         d2_h_cos2=di2+(h_cos*h_cos);
 993         frac=ci2/d2_h_cos2;
 994         g=gamma_i*(1.0+ci2di2-frac);
 995         dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
 996
 997         /* zeta sum += f_c_ik * g_ijk */
 998         if(d_ik<=Rk) {
 999                 zeta_ij+=g;
1000                 f_c_ik=1.0;
1001                 df_c_ik=0.0;
1002         }
1003         else {
1004                 s_r=Sk-Rk;
1005                 arg=M_PI*(d_ik-Rk)/s_r;
1006                 f_c_ik=0.5+0.5*cos(arg);
1007                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
1008                 zeta_ij+=f_c_ik*g;
1009         }
1010
1011         /* store even more data for second k loop */
1012         exchange->g[kcount]=g;
1013         exchange->dg[kcount]=dg;
1014         exchange->d_ik[kcount]=d_ik;
1015         exchange->cos_theta[kcount]=cos_theta;
1016         exchange->f_c_ik[kcount]=f_c_ik;
1017         exchange->df_c_ik[kcount]=df_c_ik;
1018
1019         /* increase k counter */
1020         kcount++;
1021
1022 #ifdef STATIC_LISTS
1023                                         }
1024 #elif LOWMEM_LISTS
1025                                         }
1026 #else
1027                                         } while(list_next_f(that)!=\
1028                                                 L_NO_NEXT_ELEMENT);
1029 #endif
1030
1031                                 }
1032
1033 /* j2 func here ... */
1034
1035
1036         if(brand_i==jtom->brand) {
1037                 S=params->S[brand_i];
1038                 R=params->R[brand_i];
1039                 B=params->B[brand_i];
1040                 A=params->A[brand_i];
1041                 r0=params->r0[brand_i];
1042                 mu=params->mu[brand_i];
1043                 lambda=params->lambda[brand_i];
1044         }
1045         else {
1046                 S=params->Smixed;
1047                 R=params->Rmixed;
1048                 B=params->Bmixed;
1049                 A=params->Amixed;
1050                 r0=params->r0_mixed;
1051                 mu=params->mu_m;
1052                 lambda=params->lambda_m;
1053         }
1054
1055         /* f_c, df_c */
1056         if(d_ij<R) {
1057                 f_c=1.0;
1058                 df_c=0.0;
1059         }
1060         else {
1061                 s_r=S-R;
1062                 arg=M_PI*(d_ij-R)/s_r;
1063                 f_c=0.5+0.5*cos(arg);
1064                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
1065         }
1066
1067         /* f_a, df_a */
1068         f_a=-B*exp(-mu*(d_ij-r0));
1069         df_a=mu*f_a/d_ij;
1070
1071         /* f_r, df_r */
1072         f_r=A*exp(-lambda*(d_ij-r0));
1073         df_r=lambda*f_r/d_ij;
1074
1075         /* b, db */
1076         if(zeta_ij==0.0) {
1077                 b=1.0;
1078                 db=0.0;
1079         }
1080         else {
1081                 b=1.0/sqrt(1.0+zeta_ij);
1082                 db=-0.5*b/(1.0+zeta_ij);
1083         }
1084
1085         /* force contribution for atom i */
1086         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
1087         v3_scale(&force,&(dist_ij),scale);
1088         pthread_mutex_lock(&(amutex[ai->tag]));
1089         v3_add(&(ai->f),&(ai->f),&force);
1090         pthread_mutex_unlock(&(amutex[ai->tag]));
1091
1092         /* force contribution for atom j */
1093         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
1094         pthread_mutex_lock(&(amutex[jtom->tag]));
1095         v3_add(&(jtom->f),&(jtom->f),&force);
1096         pthread_mutex_unlock(&(amutex[jtom->tag]));
1097
1098         /* virial */
1099         pthread_mutex_lock(&(amutex[ai->tag]));
1100         albe_v_calc(ai,&force,&(dist_ij));
1101         //virial_calc(ai,&force,&(dist_ij));
1102         pthread_mutex_unlock(&(amutex[ai->tag]));
1103
1104 #ifdef DEBUG
1105 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1106         if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
1107                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
1108                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
1109                 if(ai==&(moldyn->atom[0]))
1110                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
1111                 if(jtom==&(moldyn->atom[0]))
1112                         printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1113                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
1114                                                     f_c,b,f_a,f_r);
1115                 printf("          %f %f %f\n",zeta_ij,.0,.0);
1116         }
1117 }
1118 #endif
1119
1120         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
1121         pre_dzeta=0.5*f_a*f_c*db;
1122
1123         /* energy contribution */
1124         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
1125         pthread_mutex_lock(&emutex);
1126         moldyn->energy+=energy;
1127         pthread_mutex_unlock(&emutex);
1128         pthread_mutex_lock(&(amutex[ai->tag]));
1129         ai->e+=energy;
1130         pthread_mutex_unlock(&(amutex[ai->tag]));
1131
1132         /* reset k counter for second k loop */
1133         kcount=0;
1134
1135
1136                                 /* second loop over atoms k */
1137                                 for(k=0;k<27;k++) {
1138
1139                                         bc_ik=(k<dnlc)?0:1;
1140 #ifdef STATIC_LISTS
1141                                         q=0;
1142
1143                                         while(neighbour_i[k][q]!=-1) {
1144
1145                                                 ktom=&(itom[neighbour_i[k][q]]);
1146                                                 q++;
1147 #elif LOWMEM_LISTS
1148                                         q=neighbour_i[k];
1149
1150                                         while(q!=-1) {
1151
1152                                                 ktom=&(itom[q]);
1153                                                 q=lc->subcell->list[q];
1154 #else
1155                                         that=&(neighbour_i2[k]);
1156                                         list_reset_f(that);
1157
1158                                         if(that->start==NULL)
1159                                                 continue;
1160
1161                                         do {
1162                                                 ktom=that->current->data;
1163 #endif
1164
1165                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
1166                                                         continue;
1167
1168                                                 if(ktom==jtom)
1169                                                         continue;
1170
1171                                                 if(ktom==&(itom[i]))
1172                                                         continue;
1173
1174
1175 /* k2 func here ... */
1176 /* albe 3 body potential function (second k loop) */
1177
1178         if(kcount>ALBE_MAXN)
1179                 printf("FATAL: neighbours!\n");
1180
1181         /* d_ik2 */
1182         d_ik2=exchange->d_ik2[kcount];
1183
1184         if(brand_i==ktom->brand)
1185                 Sk2=params->S2[brand_i];
1186         else
1187                 Sk2=params->S2mixed;
1188
1189         /* return if d_ik > S */
1190         if(d_ik2>Sk2) {
1191                 kcount++;
1192                 continue;
1193         }
1194
1195         /* dist_ik, d_ik */
1196         dist_ik=exchange->dist_ik[kcount];
1197         d_ik=exchange->d_ik[kcount];
1198
1199         /* f_c_ik, df_c_ik */
1200         f_c_ik=exchange->f_c_ik[kcount];
1201         df_c_ik=exchange->df_c_ik[kcount];
1202
1203         /* g, dg, cos_theta */
1204         g=exchange->g[kcount];
1205         dg=exchange->dg[kcount];
1206         cos_theta=exchange->cos_theta[kcount];
1207
1208         /* cos_theta derivatives wrt j,k */
1209         dijdik_inv=1.0/(d_ij*d_ik);
1210         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
1211         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
1212         v3_add(&dcosdrj,&dcosdrj,&tmp);
1213         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
1214         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
1215         v3_add(&dcosdrk,&dcosdrk,&tmp);
1216
1217         /* f_c_ik * dg, df_c_ik * g */
1218         fcdg=f_c_ik*dg;
1219         dfcg=df_c_ik*g;
1220
1221         /* derivative wrt j */
1222         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
1223
1224         /* force contribution */
1225         pthread_mutex_lock(&(amutex[jtom->tag]));
1226         v3_add(&(jtom->f),&(jtom->f),&force);
1227         pthread_mutex_unlock(&(amutex[jtom->tag]));
1228
1229 #ifdef DEBUG
1230 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1231         if(jtom==&(moldyn->atom[DATOM])) {
1232                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1233                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
1234                 printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1235                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1236                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
1237         }
1238 }
1239 #endif
1240
1241         /* virial */
1242         pthread_mutex_lock(&(amutex[ai->tag]));
1243         albe_v_calc(ai,&force,&dist_ij);
1244         //virial_calc(ai,&force,&dist_ij);
1245
1246         /* force contribution to atom i */
1247         v3_scale(&force,&force,-1.0);
1248         v3_add(&(ai->f),&(ai->f),&force);
1249         pthread_mutex_unlock(&(amutex[ai->tag]));
1250
1251         /* derivative wrt k */
1252         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
1253         v3_scale(&tmp,&dcosdrk,fcdg);
1254         v3_add(&force,&force,&tmp);
1255         v3_scale(&force,&force,pre_dzeta);
1256
1257         /* force contribution */
1258         pthread_mutex_lock(&(amutex[ktom->tag]));
1259         v3_add(&(ktom->f),&(ktom->f),&force);
1260         pthread_mutex_unlock(&(amutex[ktom->tag]));
1261
1262 #ifdef DEBUG
1263 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1264         if(ktom==&(moldyn->atom[DATOM])) {
1265                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1266                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
1267                 printf("  total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
1268                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1269                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
1270         }
1271 }
1272 #endif
1273
1274         /* virial */
1275         pthread_mutex_lock(&(amutex[ai->tag]));
1276         albe_v_calc(ai,&force,&dist_ik);
1277         //virial_calc(ai,&force,&dist_ik);
1278
1279         /* force contribution to atom i */
1280         v3_scale(&force,&force,-1.0);
1281         v3_add(&(ai->f),&(ai->f),&force);
1282         pthread_mutex_unlock(&(amutex[ai->tag]));
1283
1284         /* increase k counter */
1285         kcount++;
1286
1287
1288
1289 #ifdef STATIC_LISTS
1290                                         }
1291 #elif LOWMEM_LISTS
1292                                         }
1293 #else
1294                                         } while(list_next_f(that)!=\
1295                                                 L_NO_NEXT_ELEMENT);
1296 #endif
1297
1298                                 }
1299
1300 #ifdef STATIC_LISTS
1301                         }
1302 #elif LOWMEM_LISTS
1303                         }
1304 #else
1305                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1306 #endif
1307
1308                 }
1309
1310         } // i loop
1311
1312 #ifdef DEBUG
1313         //printf("\n\n");
1314 #endif
1315 #ifdef VDEBUG
1316         printf("\n\n");
1317 #endif
1318
1319 #ifdef DEBUG
1320         //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
1321         if(moldyn->time>DSTART&&moldyn->time<DEND) {
1322                 printf("force:\n");
1323                 printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
1324                 printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
1325                 printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
1326         }
1327 #endif
1328
1329         pthread_exit(NULL);
1330
1331         return 0;
1332 }
1333
1334 int albe_potential_force_calc(t_moldyn *moldyn) {
1335
1336         int i,j,ret;
1337         t_pft_data pft_data[MAX_THREADS];
1338         int count;
1339         pthread_t pft_thread[MAX_THREADS];
1340         t_atom *itom;
1341         t_virial *virial;
1342
1343         count=moldyn->count;
1344         itom=moldyn->atom;
1345
1346         /* reset energy */
1347         moldyn->energy=0.0;
1348
1349         /* reset global virial */
1350         memset(&(moldyn->gvir),0,sizeof(t_virial));
1351
1352         /* reset force, site energy and virial of every atom */
1353         for(i=0;i<count;i++) {
1354
1355                 /* reset force */
1356                 v3_zero(&(itom[i].f));
1357
1358                 /* reset virial */
1359                 virial=&(itom[i].virial);
1360                 virial->xx=0.0;
1361                 virial->yy=0.0;
1362                 virial->zz=0.0;
1363                 virial->xy=0.0;
1364                 virial->xz=0.0;
1365                 virial->yz=0.0;
1366
1367                 /* reset site energy */
1368                 itom[i].e=0.0;
1369
1370         }
1371
1372         /* start threads */
1373         for(j=0;j<MAX_THREADS;j++) {
1374
1375                 /* prepare thread data */
1376                 pft_data[j].moldyn=moldyn;
1377                 pft_data[j].start=j*(count/MAX_THREADS);
1378                 if(j==MAX_THREADS-1) {
1379                         pft_data[j].end=count;
1380                 }
1381                 else {
1382                         pft_data[j].end=pft_data[j].start;
1383                         pft_data[j].end+=count/MAX_THREADS;
1384                 }
1385
1386                 ret=pthread_create(&(pft_thread[j]),NULL,
1387                                    potential_force_thread,
1388                                    &(pft_data[j]));
1389                 if(ret)  {
1390                         perror("[albe fast] pf thread create");
1391                         return ret;
1392                 }
1393         }
1394
1395         /* join threads */
1396         for(j=0;j<MAX_THREADS;j++) {
1397
1398                 ret=pthread_join(pft_thread[j],NULL);
1399                 if(ret) {
1400                         perror("[albe fast] pf thread join");
1401                         return ret;
1402                 }
1403         }
1404
1405         /* some postprocessing */
1406         for(i=0;i<count;i++) {
1407                 /* calculate global virial */
1408                 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
1409                 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
1410                 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
1411                 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
1412                 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
1413                 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
1414
1415                 /* check forces regarding the given timestep */
1416                 if(v3_norm(&(itom[i].f))>\
1417                    0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
1418                         printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
1419                                i);
1420         }
1421
1422         return 0;
1423 }
1424
1425 #endif // PTHREADS