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