]> hackdaworld.org Git - physik/posic.git/blob - potentials/albe_fast.c
added partition support in create_lattice function
[physik/posic.git] / 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