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 #ifdef MATTONI
 456         scale=-0.5*(f_c*(df_r-b*df_a)); // - in albe formalism
 457 #else
 458         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
 459 #endif
 460         v3_scale(&force,&(dist_ij),scale);
 461         v3_add(&(ai->f),&(ai->f),&force);
 462
 463         /* force contribution for atom j */
 464         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
 465         v3_add(&(jtom->f),&(jtom->f),&force);
 466
 467         /* virial */
 468         albe_v_calc(ai,&force,&(dist_ij));
 469         //virial_calc(ai,&force,&(dist_ij));
 470
 471 #ifdef DEBUG
 472 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 473         if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
 474                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
 475                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 476                 if(ai==&(moldyn->atom[0]))
 477                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
 478                 if(jtom==&(moldyn->atom[0]))
 479                         printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
 480                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
 481                                                     f_c,b,f_a,f_r);
 482                 printf("          %f %f %f\n",zeta_ij,.0,.0);
 483         }
 484 }
 485 #endif
 486
 487         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
 488         pre_dzeta=0.5*f_a*f_c*db;
 489
 490         /* energy contribution */
 491         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
 492         moldyn->energy+=energy;
 493         ai->e+=energy;
 494
 495         /* reset k counter for second k loop */
 496         kcount=0;
 497
 498
 499                                 /* second loop over atoms k */
 500                                 for(k=0;k<27;k++) {
 501
 502                                         bc_ik=(k<dnlc)?0:1;
 503 #ifdef STATIC_LISTS
 504                                         q=0;
 505
 506                                         while(neighbour_i[k][q]!=-1) {
 507
 508                                                 ktom=&(itom[neighbour_i[k][q]]);
 509                                                 q++;
 510 #elif LOWMEM_LISTS
 511                                         q=neighbour_i[k];
 512
 513                                         while(q!=-1) {
 514
 515                                                 ktom=&(itom[q]);
 516                                                 q=lc->subcell->list[q];
 517 #else
 518                                         that=&(neighbour_i2[k]);
 519                                         list_reset_f(that);
 520
 521                                         if(that->start==NULL)
 522                                                 continue;
 523
 524                                         do {
 525                                                 ktom=that->current->data;
 526 #endif
 527
 528                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
 529                                                         continue;
 530
 531                                                 if(ktom==jtom)
 532                                                         continue;
 533
 534                                                 if(ktom==&(itom[i]))
 535                                                         continue;
 536
 537
 538 /* k2 func here ... */
 539 /* albe 3 body potential function (second k loop) */
 540
 541         if(kcount>ALBE_MAXN)
 542                 printf("FATAL: neighbours!\n");
 543
 544         /* d_ik2 */
 545         d_ik2=exchange->d_ik2[kcount];
 546
 547         if(brand_i==ktom->brand)
 548                 Sk2=params->S2[brand_i];
 549         else
 550                 Sk2=params->S2mixed;
 551
 552         /* return if d_ik > S */
 553         if(d_ik2>Sk2) {
 554                 kcount++;
 555                 continue;
 556         }
 557
 558         /* dist_ik, d_ik */
 559         dist_ik=exchange->dist_ik[kcount];
 560         d_ik=exchange->d_ik[kcount];
 561
 562         /* f_c_ik, df_c_ik */
 563         f_c_ik=exchange->f_c_ik[kcount];
 564         df_c_ik=exchange->df_c_ik[kcount];
 565
 566         /* g, dg, cos_theta */
 567         g=exchange->g[kcount];
 568         dg=exchange->dg[kcount];
 569         cos_theta=exchange->cos_theta[kcount];
 570
 571         /* cos_theta derivatives wrt j,k */
 572         dijdik_inv=1.0/(d_ij*d_ik);
 573         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
 574         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
 575         v3_add(&dcosdrj,&dcosdrj,&tmp);
 576         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
 577         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
 578         v3_add(&dcosdrk,&dcosdrk,&tmp);
 579
 580         /* f_c_ik * dg, df_c_ik * g */
 581         fcdg=f_c_ik*dg;
 582         dfcg=df_c_ik*g;
 583
 584         /* derivative wrt j */
 585         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 586
 587         /* force contribution */
 588         v3_add(&(jtom->f),&(jtom->f),&force);
 589
 590 #ifdef DEBUG
 591 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 592         if(jtom==&(moldyn->atom[DATOM])) {
 593                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
 594                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 595                 printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
 596                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 597                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 598         }
 599 }
 600 #endif
 601
 602         /* virial */
 603         albe_v_calc(ai,&force,&dist_ij);
 604         //virial_calc(ai,&force,&dist_ij);
 605
 606         /* force contribution to atom i */
 607         v3_scale(&force,&force,-1.0);
 608         v3_add(&(ai->f),&(ai->f),&force);
 609
 610         /* derivative wrt k */
 611 #ifdef MATTONI
 612         v3_scale(&tmp,&dcosdrk,fcdg);
 613         v3_scale(&force,&tmp,pre_dzeta);
 614 #else
 615         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
 616         v3_scale(&tmp,&dcosdrk,fcdg);
 617         v3_add(&force,&force,&tmp);
 618         v3_scale(&force,&force,pre_dzeta);
 619 #endif
 620
 621         /* force contribution */
 622         v3_add(&(ktom->f),&(ktom->f),&force);
 623
 624 #ifdef DEBUG
 625 if(moldyn->time>DSTART&&moldyn->time<DEND) {
 626         if(ktom==&(moldyn->atom[DATOM])) {
 627                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
 628                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
 629                 printf("  total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
 630                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
 631                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
 632         }
 633 }
 634 #endif
 635
 636         /* virial */
 637         albe_v_calc(ai,&force,&dist_ik);
 638         //virial_calc(ai,&force,&dist_ik);
 639
 640         /* force contribution to atom i */
 641         v3_scale(&force,&force,-1.0);
 642         v3_add(&(ai->f),&(ai->f),&force);
 643
 644         /* increase k counter */
 645         kcount++;
 646
 647
 648
 649 #ifdef STATIC_LISTS
 650                                         }
 651 #elif LOWMEM_LISTS
 652                                         }
 653 #else
 654                                         } while(list_next_f(that)!=\
 655                                                 L_NO_NEXT_ELEMENT);
 656 #endif
 657
 658                                 }
 659
 660 #ifdef STATIC_LISTS
 661                         }
 662 #elif LOWMEM_LISTS
 663                         }
 664 #else
 665                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 666 #endif
 667
 668                 }
 669
 670 #ifdef DEBUG
 671         //printf("\n\n");
 672 #endif
 673 #ifdef VDEBUG
 674         printf("\n\n");
 675 #endif
 676
 677         }
 678
 679 #ifdef DEBUG
 680         //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
 681         if(moldyn->time>DSTART&&moldyn->time<DEND) {
 682                 printf("force:\n");
 683                 printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
 684                 printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
 685                 printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
 686         }
 687 #endif
 688
 689         /* some postprocessing */
 690 #ifdef PARALLEL
 691         #pragma omp parallel for
 692 #endif
 693         for(i=0;i<count;i++) {
 694                 /* calculate global virial */
 695                 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
 696                 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
 697                 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
 698                 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
 699                 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
 700                 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
 701
 702                 /* check forces regarding the given timestep */
 703                 if(v3_norm(&(itom[i].f))>\
 704                    0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
 705                         printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
 706                                i);
 707         }
 708
 709         return 0;
 710 }
 711
 712
 713 #else // PTHREADS
 714
 715
 716 typedef struct s_pft_data {
 717         t_moldyn *moldyn;
 718         int start,end;
 719 } t_pft_data;
 720
 721 void *potential_force_thread(void *ptr) {
 722
 723         t_pft_data *pft_data;
 724         t_moldyn *moldyn;
 725         t_albe_exchange ec;
 726
 727         int i,j,k,count;
 728         t_atom *itom,*jtom,*ktom;
 729         t_linkcell *lc;
 730 #ifdef STATIC_LISTS
 731         int *neighbour_i[27];
 732         int p,q;
 733 #elif LOWMEM_LISTS
 734         int neighbour_i[27];
 735         int p,q;
 736 #else
 737         t_list neighbour_i[27];
 738         t_list neighbour_i2[27];
 739         t_list *this,*that;
 740 #endif
 741         u8 bc_ij,bc_ik;
 742         int dnlc;
 743
 744         // needed to work
 745         t_atom *ai;
 746
 747         // optimized
 748         t_albe_mult_params *params;
 749         t_albe_exchange *exchange;
 750         t_3dvec dist_ij;
 751         double d_ij2;
 752         double d_ij;
 753         u8 brand_i;
 754         double S2;
 755         int kcount;
 756         double zeta_ij;
 757         double pre_dzeta;
 758
 759         // more ...
 760         double Rk,Sk,Sk2;
 761         t_3dvec dist_ik;
 762         double d_ik2,d_ik;
 763         double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
 764         double f_c_ik,df_c_ik;
 765
 766         t_3dvec force;
 767         double f_a,df_a,b,db,f_c,df_c;
 768         double f_r,df_r;
 769         double scale;
 770         double mu,B;
 771         double lambda,A;
 772         double r0;
 773         double S,R;
 774         double energy;
 775
 776         double dijdik_inv,fcdg,dfcg;
 777         t_3dvec dcosdrj,dcosdrk;
 778         t_3dvec tmp;
 779
 780         // even more ...
 781         double gamma_i;
 782         double c_i;
 783         double d_i;
 784         double h_i;
 785         double ci2;
 786         double di2;
 787         double ci2di2;
 788
 789         pft_data=ptr;
 790         moldyn=pft_data->moldyn;
 791         exchange=&ec;
 792
 793         count=moldyn->count;
 794         itom=moldyn->atom;
 795         lc=&(moldyn->lc);
 796
 797         // optimized
 798         params=moldyn->pot_params;
 799
 800         /* get energy, force and virial for atoms */
 801
 802         for(i=pft_data->start;i<pft_data->end;i++) {
 803
 804                 if(!(itom[i].attr&ATOM_ATTR_3BP))
 805                         return 0;
 806
 807                 // thread safe this way!
 808                 dnlc=link_cell_neighbour_index(moldyn,
 809                                           (itom[i].r.x+moldyn->dim.x/2)/lc->x,
 810                                           (itom[i].r.y+moldyn->dim.y/2)/lc->y,
 811                                           (itom[i].r.z+moldyn->dim.z/2)/lc->z,
 812                                           neighbour_i);
 813
 814                 /* copy the neighbour lists */
 815 #ifdef STATIC_LISTS
 816 #elif LOWMEM_LISTS
 817 #else
 818                 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
 819 #endif
 820
 821                 ai=&(itom[i]);
 822                 brand_i=ai->brand;
 823
 824                 /* loop over atoms j */
 825                 for(j=0;j<27;j++) {
 826
 827                         bc_ij=(j<dnlc)?0:1;
 828 #ifdef STATIC_LISTS
 829                         p=0;
 830
 831                         while(neighbour_i[j][p]!=-1) {
 832
 833                                 jtom=&(itom[neighbour_i[j][p]]);
 834                                 p++;
 835 #elif LOWMEM_LISTS
 836                         p=neighbour_i[j];
 837
 838                         while(p!=-1) {
 839
 840                                 jtom=&(itom[p]);
 841                                 p=lc->subcell->list[p]; // thread safe!
 842 #else
 843                         this=&(neighbour_i[j]);
 844                         list_reset_f(this);
 845
 846                         if(this->start==NULL)
 847                                 continue;
 848
 849                         do {
 850
 851                                 jtom=this->current->data;
 852 #endif
 853
 854                                 if(jtom==&(itom[i]))
 855                                         continue;
 856
 857                                 if(!(jtom->attr&ATOM_ATTR_3BP))
 858                                         continue;
 859
 860                                 /* reset 3bp run */
 861                                 moldyn->run3bp=1;
 862
 863
 864 /* j1 func here ... */
 865 /* albe 3 body potential function (first ij loop) */
 866
 867         /* reset zeta sum */
 868         zeta_ij=0.0;
 869
 870         /*
 871          * set ij depending values
 872          */
 873
 874         if(brand_i==jtom->brand) {
 875                 S2=params->S2[brand_i];
 876         }
 877         else {
 878                 S2=params->S2mixed;
 879         }
 880
 881         /* dist_ij, d_ij2 */
 882         v3_sub(&dist_ij,&(jtom->r),&(ai->r));
 883         if(bc_ij) check_per_bound(moldyn,&dist_ij);
 884         d_ij2=v3_absolute_square(&dist_ij);
 885
 886         /* if d_ij2 > S2 => no force & potential energy contribution */
 887         if(d_ij2>S2)
 888                 continue;
 889
 890         /* d_ij */
 891         d_ij=sqrt(d_ij2);
 892
 893         /* reset k counter for first k loop */
 894         kcount=0;
 895
 896                                 /* first loop over atoms k */
 897                                 for(k=0;k<27;k++) {
 898
 899                                         bc_ik=(k<dnlc)?0:1;
 900 #ifdef STATIC_LISTS
 901                                         q=0;
 902
 903                                         while(neighbour_i[k][q]!=-1) {
 904
 905                                                 ktom=&(itom[neighbour_i[k][q]]);
 906                                                 q++;
 907 #elif LOWMEM_LISTS
 908                                         q=neighbour_i[k];
 909
 910                                         while(q!=-1) {
 911
 912                                                 ktom=&(itom[q]);
 913                                                 q=lc->subcell->list[q];
 914 #else
 915                                         that=&(neighbour_i2[k]);
 916                                         list_reset_f(that);
 917
 918                                         if(that->start==NULL)
 919                                                 continue;
 920
 921                                         do {
 922                                                 ktom=that->current->data;
 923 #endif
 924
 925                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
 926                                                         continue;
 927
 928                                                 if(ktom==jtom)
 929                                                         continue;
 930
 931                                                 if(ktom==&(itom[i]))
 932                                                         continue;
 933
 934 /* k1 func here ... */
 935 /* albe 3 body potential function (first k loop) */
 936
 937         if(kcount>ALBE_MAXN) {
 938                 printf("FATAL: neighbours = %d\n",kcount);
 939                 printf("  -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
 940         }
 941
 942         /* ik constants */
 943         if(brand_i==ktom->brand) {
 944                 Rk=params->R[brand_i];
 945                 Sk=params->S[brand_i];
 946                 Sk2=params->S2[brand_i];
 947                 /* albe needs i,k depending c,d,h and gamma values */
 948                 gamma_i=params->gamma[brand_i];
 949                 c_i=params->c[brand_i];
 950                 d_i=params->d[brand_i];
 951                 h_i=params->h[brand_i];
 952                 ci2=params->c2[brand_i];
 953                 di2=params->d2[brand_i];
 954                 ci2di2=params->c2d2[brand_i];
 955         }
 956         else {
 957                 Rk=params->Rmixed;
 958                 Sk=params->Smixed;
 959                 Sk2=params->S2mixed;
 960                 /* albe needs i,k depending c,d,h and gamma values */
 961                 gamma_i=params->gamma_m;
 962                 c_i=params->c_mixed;
 963                 d_i=params->d_mixed;
 964                 h_i=params->h_mixed;
 965                 ci2=params->c2_mixed;
 966                 di2=params->d2_mixed;
 967                 ci2di2=params->c2d2_m;
 968         }
 969
 970         /* dist_ik, d_ik2 */
 971         v3_sub(&dist_ik,&(ktom->r),&(ai->r));
 972         if(bc_ik) check_per_bound(moldyn,&dist_ik);
 973         d_ik2=v3_absolute_square(&dist_ik);
 974
 975         /* store data for second k loop */
 976         exchange->dist_ik[kcount]=dist_ik;
 977         exchange->d_ik2[kcount]=d_ik2;
 978
 979         /* return if not within cutoff */
 980         if(d_ik2>Sk2) {
 981                 kcount++;
 982                 continue;
 983         }
 984
 985         /* d_ik */
 986         d_ik=sqrt(d_ik2);
 987
 988         /* cos theta */
 989         cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 990
 991         /* g_ijk
 992         h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
 993         d2_h_cos2=exchange->di2+(h_cos*h_cos);
 994         frac=exchange->ci2/d2_h_cos2;
 995         g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
 996         dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
 997         */
 998
 999         h_cos=h_i+cos_theta; // + in albe formalism
1000         d2_h_cos2=di2+(h_cos*h_cos);
1001         frac=ci2/d2_h_cos2;
1002         g=gamma_i*(1.0+ci2di2-frac);
1003         dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
1004
1005         /* zeta sum += f_c_ik * g_ijk */
1006         if(d_ik<=Rk) {
1007                 zeta_ij+=g;
1008                 f_c_ik=1.0;
1009                 df_c_ik=0.0;
1010         }
1011         else {
1012                 s_r=Sk-Rk;
1013                 arg=M_PI*(d_ik-Rk)/s_r;
1014                 f_c_ik=0.5+0.5*cos(arg);
1015                 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
1016                 zeta_ij+=f_c_ik*g;
1017         }
1018
1019         /* store even more data for second k loop */
1020         exchange->g[kcount]=g;
1021         exchange->dg[kcount]=dg;
1022         exchange->d_ik[kcount]=d_ik;
1023         exchange->cos_theta[kcount]=cos_theta;
1024         exchange->f_c_ik[kcount]=f_c_ik;
1025         exchange->df_c_ik[kcount]=df_c_ik;
1026
1027         /* increase k counter */
1028         kcount++;
1029
1030 #ifdef STATIC_LISTS
1031                                         }
1032 #elif LOWMEM_LISTS
1033                                         }
1034 #else
1035                                         } while(list_next_f(that)!=\
1036                                                 L_NO_NEXT_ELEMENT);
1037 #endif
1038
1039                                 }
1040
1041 /* j2 func here ... */
1042
1043
1044         if(brand_i==jtom->brand) {
1045                 S=params->S[brand_i];
1046                 R=params->R[brand_i];
1047                 B=params->B[brand_i];
1048                 A=params->A[brand_i];
1049                 r0=params->r0[brand_i];
1050                 mu=params->mu[brand_i];
1051                 lambda=params->lambda[brand_i];
1052         }
1053         else {
1054                 S=params->Smixed;
1055                 R=params->Rmixed;
1056                 B=params->Bmixed;
1057                 A=params->Amixed;
1058                 r0=params->r0_mixed;
1059                 mu=params->mu_m;
1060                 lambda=params->lambda_m;
1061         }
1062
1063         /* f_c, df_c */
1064         if(d_ij<R) {
1065                 f_c=1.0;
1066                 df_c=0.0;
1067         }
1068         else {
1069                 s_r=S-R;
1070                 arg=M_PI*(d_ij-R)/s_r;
1071                 f_c=0.5+0.5*cos(arg);
1072                 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
1073         }
1074
1075         /* f_a, df_a */
1076         f_a=-B*exp(-mu*(d_ij-r0));
1077         df_a=mu*f_a/d_ij;
1078
1079         /* f_r, df_r */
1080         f_r=A*exp(-lambda*(d_ij-r0));
1081         df_r=lambda*f_r/d_ij;
1082
1083         /* b, db */
1084         if(zeta_ij==0.0) {
1085                 b=1.0;
1086                 db=0.0;
1087         }
1088         else {
1089                 b=1.0/sqrt(1.0+zeta_ij);
1090                 db=-0.5*b/(1.0+zeta_ij);
1091         }
1092
1093         /* force contribution for atom i */
1094 #ifdef MATTONI
1095         scale=-0.5*(f_c*(df_r-b*df_a)); // - in albe formalism
1096 #else
1097         scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
1098 #endif
1099         v3_scale(&force,&(dist_ij),scale);
1100         if(pthread_mutex_lock(&(amutex[ai->tag])))
1101                 perror("[albe fast] mutex lock (1)\n");
1102         v3_add(&(ai->f),&(ai->f),&force);
1103         if(pthread_mutex_unlock(&(amutex[ai->tag])))
1104                 perror("[albe fast] mutex unlock (1)\n");
1105
1106         /* force contribution for atom j */
1107         v3_scale(&force,&force,-1.0); // dri rij = - drj rij
1108         if(pthread_mutex_lock(&(amutex[jtom->tag])))
1109                 perror("[albe fast] mutex lock (2)\n");
1110         v3_add(&(jtom->f),&(jtom->f),&force);
1111         if(pthread_mutex_unlock(&(amutex[jtom->tag])))
1112                 perror("[albe fast] mutex unlock (2)\n");
1113
1114         /* virial */
1115         if(pthread_mutex_lock(&(amutex[ai->tag])))
1116                 perror("[albe fast] mutex lock (3)\n");
1117         albe_v_calc(ai,&force,&(dist_ij));
1118         //virial_calc(ai,&force,&(dist_ij));
1119         if(pthread_mutex_unlock(&(amutex[ai->tag])))
1120                 perror("[albe fast] mutex unlock (3)\n");
1121
1122 #ifdef DEBUG
1123 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1124         if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
1125                 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
1126                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
1127                 if(ai==&(moldyn->atom[0]))
1128                         printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
1129                 if(jtom==&(moldyn->atom[0]))
1130                         printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1131                 printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
1132                                                     f_c,b,f_a,f_r);
1133                 printf("          %f %f %f\n",zeta_ij,.0,.0);
1134         }
1135 }
1136 #endif
1137
1138         /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
1139         pre_dzeta=0.5*f_a*f_c*db;
1140
1141         /* energy contribution */
1142         energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
1143         if(pthread_mutex_lock(&emutex))
1144                 perror("[albe fast] mutex lock (energy)\n");
1145         moldyn->energy+=energy;
1146         if(pthread_mutex_unlock(&emutex))
1147                 perror("[albe fast] mutex unlock (energy)\n");
1148         if(pthread_mutex_lock(&(amutex[ai->tag])))
1149                 perror("[albe fast] mutex lock (4)\n");
1150         ai->e+=energy;
1151         if(pthread_mutex_unlock(&(amutex[ai->tag])))
1152                 perror("[albe fast] mutex unlock (4)\n");
1153
1154         /* reset k counter for second k loop */
1155         kcount=0;
1156
1157
1158                                 /* second loop over atoms k */
1159                                 for(k=0;k<27;k++) {
1160
1161                                         bc_ik=(k<dnlc)?0:1;
1162 #ifdef STATIC_LISTS
1163                                         q=0;
1164
1165                                         while(neighbour_i[k][q]!=-1) {
1166
1167                                                 ktom=&(itom[neighbour_i[k][q]]);
1168                                                 q++;
1169 #elif LOWMEM_LISTS
1170                                         q=neighbour_i[k];
1171
1172                                         while(q!=-1) {
1173
1174                                                 ktom=&(itom[q]);
1175                                                 q=lc->subcell->list[q];
1176 #else
1177                                         that=&(neighbour_i2[k]);
1178                                         list_reset_f(that);
1179
1180                                         if(that->start==NULL)
1181                                                 continue;
1182
1183                                         do {
1184                                                 ktom=that->current->data;
1185 #endif
1186
1187                                                 if(!(ktom->attr&ATOM_ATTR_3BP))
1188                                                         continue;
1189
1190                                                 if(ktom==jtom)
1191                                                         continue;
1192
1193                                                 if(ktom==&(itom[i]))
1194                                                         continue;
1195
1196
1197 /* k2 func here ... */
1198 /* albe 3 body potential function (second k loop) */
1199
1200         if(kcount>ALBE_MAXN)
1201                 printf("FATAL: neighbours!\n");
1202
1203         /* d_ik2 */
1204         d_ik2=exchange->d_ik2[kcount];
1205
1206         if(brand_i==ktom->brand)
1207                 Sk2=params->S2[brand_i];
1208         else
1209                 Sk2=params->S2mixed;
1210
1211         /* return if d_ik > S */
1212         if(d_ik2>Sk2) {
1213                 kcount++;
1214                 continue;
1215         }
1216
1217         /* dist_ik, d_ik */
1218         dist_ik=exchange->dist_ik[kcount];
1219         d_ik=exchange->d_ik[kcount];
1220
1221         /* f_c_ik, df_c_ik */
1222         f_c_ik=exchange->f_c_ik[kcount];
1223         df_c_ik=exchange->df_c_ik[kcount];
1224
1225         /* g, dg, cos_theta */
1226         g=exchange->g[kcount];
1227         dg=exchange->dg[kcount];
1228         cos_theta=exchange->cos_theta[kcount];
1229
1230         /* cos_theta derivatives wrt j,k */
1231         dijdik_inv=1.0/(d_ij*d_ik);
1232         v3_scale(&dcosdrj,&dist_ik,dijdik_inv);         // j
1233         v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
1234         v3_add(&dcosdrj,&dcosdrj,&tmp);
1235         v3_scale(&dcosdrk,&dist_ij,dijdik_inv);         // k
1236         v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
1237         v3_add(&dcosdrk,&dcosdrk,&tmp);
1238
1239         /* f_c_ik * dg, df_c_ik * g */
1240         fcdg=f_c_ik*dg;
1241         dfcg=df_c_ik*g;
1242
1243         /* derivative wrt j */
1244         v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
1245
1246         /* force contribution */
1247         if(pthread_mutex_lock(&(amutex[jtom->tag])))
1248                 perror("[albe fast] mutex lock (5)\n");
1249         v3_add(&(jtom->f),&(jtom->f),&force);
1250         if(pthread_mutex_unlock(&(amutex[jtom->tag])))
1251                 perror("[albe fast] mutex unlock (5)\n");
1252
1253 #ifdef DEBUG
1254 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1255         if(jtom==&(moldyn->atom[DATOM])) {
1256                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1257                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
1258                 printf("  total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1259                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1260                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
1261         }
1262 }
1263 #endif
1264
1265         /* virial */
1266         if(pthread_mutex_lock(&(amutex[ai->tag])))
1267                 perror("[albe fast] mutex lock (6)\n");
1268         albe_v_calc(ai,&force,&dist_ij);
1269         //virial_calc(ai,&force,&dist_ij);
1270
1271         /* force contribution to atom i */
1272         v3_scale(&force,&force,-1.0);
1273         v3_add(&(ai->f),&(ai->f),&force);
1274         if(pthread_mutex_unlock(&(amutex[ai->tag])))
1275                 perror("[albe fast] mutex unlock (6)\n");
1276
1277         /* derivative wrt k */
1278 #ifdef MATTONI
1279         v3_scale(&tmp,&dcosdrk,fcdg);
1280         v3_scale(&force,&tmp,pre_dzeta);
1281 #else
1282         v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
1283         v3_scale(&tmp,&dcosdrk,fcdg);
1284         v3_add(&force,&force,&tmp);
1285         v3_scale(&force,&force,pre_dzeta);
1286 #endif
1287
1288         /* force contribution */
1289         if(pthread_mutex_lock(&(amutex[ktom->tag])))
1290                 perror("[albe fast] mutex lock (7)\n");
1291         v3_add(&(ktom->f),&(ktom->f),&force);
1292         if(pthread_mutex_unlock(&(amutex[ktom->tag])))
1293                 perror("[albe fast] mutex unlock (7)\n");
1294
1295 #ifdef DEBUG
1296 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1297         if(ktom==&(moldyn->atom[DATOM])) {
1298                 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1299                 printf("  adding %f %f %f\n",force.x,force.y,force.z);
1300                 printf("  total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
1301                 printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1302                 printf("    d ij ik = %f %f\n",d_ij,d_ik);
1303         }
1304 }
1305 #endif
1306
1307         /* virial */
1308         if(pthread_mutex_lock(&(amutex[ai->tag])))
1309                 perror("[albe fast] mutex lock (8)\n");
1310         albe_v_calc(ai,&force,&dist_ik);
1311         //virial_calc(ai,&force,&dist_ik);
1312
1313         /* force contribution to atom i */
1314         v3_scale(&force,&force,-1.0);
1315         v3_add(&(ai->f),&(ai->f),&force);
1316         if(pthread_mutex_unlock(&(amutex[ai->tag])))
1317                 perror("[albe fast] mutex unlock (8)\n");
1318
1319         /* increase k counter */
1320         kcount++;
1321
1322
1323
1324 #ifdef STATIC_LISTS
1325                                         }
1326 #elif LOWMEM_LISTS
1327                                         }
1328 #else
1329                                         } while(list_next_f(that)!=\
1330                                                 L_NO_NEXT_ELEMENT);
1331 #endif
1332
1333                                 }
1334
1335 #ifdef STATIC_LISTS
1336                         }
1337 #elif LOWMEM_LISTS
1338                         }
1339 #else
1340                         } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1341 #endif
1342
1343                 }
1344
1345         } // i loop
1346
1347 #ifdef DEBUG
1348         //printf("\n\n");
1349 #endif
1350 #ifdef VDEBUG
1351         printf("\n\n");
1352 #endif
1353
1354 #ifdef DEBUG
1355         //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
1356         if(moldyn->time>DSTART&&moldyn->time<DEND) {
1357                 printf("force:\n");
1358                 printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
1359                 printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
1360                 printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
1361         }
1362 #endif
1363
1364         pthread_exit(NULL);
1365
1366         return 0;
1367 }
1368
1369 int albe_potential_force_calc(t_moldyn *moldyn) {
1370
1371         int i,j,ret;
1372         t_pft_data pft_data[MAX_THREADS];
1373         int count;
1374         pthread_t pft_thread[MAX_THREADS];
1375         t_atom *itom;
1376         t_virial *virial;
1377
1378         count=moldyn->count;
1379         itom=moldyn->atom;
1380
1381         /* reset energy */
1382         moldyn->energy=0.0;
1383
1384         /* reset global virial */
1385         memset(&(moldyn->gvir),0,sizeof(t_virial));
1386
1387         /* reset force, site energy and virial of every atom */
1388         for(i=0;i<count;i++) {
1389
1390                 /* reset force */
1391                 v3_zero(&(itom[i].f));
1392
1393                 /* reset virial */
1394                 virial=&(itom[i].virial);
1395                 virial->xx=0.0;
1396                 virial->yy=0.0;
1397                 virial->zz=0.0;
1398                 virial->xy=0.0;
1399                 virial->xz=0.0;
1400                 virial->yz=0.0;
1401
1402                 /* reset site energy */
1403                 itom[i].e=0.0;
1404
1405         }
1406
1407         /* start threads */
1408         for(j=0;j<MAX_THREADS;j++) {
1409
1410                 /* prepare thread data */
1411                 pft_data[j].moldyn=moldyn;
1412                 pft_data[j].start=j*(count/MAX_THREADS);
1413                 if(j==MAX_THREADS-1) {
1414                         pft_data[j].end=count;
1415                 }
1416                 else {
1417                         pft_data[j].end=pft_data[j].start;
1418                         pft_data[j].end+=count/MAX_THREADS;
1419                 }
1420
1421                 ret=pthread_create(&(pft_thread[j]),NULL,
1422                                    potential_force_thread,
1423                                    &(pft_data[j]));
1424                 if(ret)  {
1425                         perror("[albe fast] pf thread create");
1426                         return ret;
1427                 }
1428         }
1429
1430         /* join threads */
1431         for(j=0;j<MAX_THREADS;j++) {
1432
1433                 ret=pthread_join(pft_thread[j],NULL);
1434                 if(ret) {
1435                         perror("[albe fast] pf thread join");
1436                         return ret;
1437                 }
1438         }
1439
1440         /* some postprocessing */
1441         for(i=0;i<count;i++) {
1442                 /* calculate global virial */
1443                 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
1444                 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
1445                 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
1446                 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
1447                 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
1448                 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
1449
1450                 /* check forces regarding the given timestep */
1451                 if(v3_norm(&(itom[i].f))>\
1452                    0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
1453                         printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
1454                                i);
1455         }
1456
1457         return 0;
1458 }
1459
1460 #endif // PTHREADS