4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
12 #include <sys/types.h>
28 #include "../moldyn.h"
29 #include "../math/math.h"
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
45 int albe_potential_force_calc(t_moldyn *moldyn) {
48 t_atom *itom,*jtom,*ktom;
58 t_list neighbour_i[27];
59 t_list neighbour_i2[27];
67 pthread_t kthread[27];
74 t_albe_mult_params *params;
75 t_albe_exchange *exchange;
89 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
90 double f_c_ik,df_c_ik;
93 double f_a,df_a,b,db,f_c,df_c;
102 double dijdik_inv,fcdg,dfcg;
103 t_3dvec dcosdrj,dcosdrk;
120 params=moldyn->pot_params;
121 exchange=&(params->exchange);
127 /* reset global virial */
128 memset(&(moldyn->gvir),0,sizeof(t_virial));
130 /* reset force, site energy and virial of every atom */
132 #pragma omp parallel for private(virial)
134 for(i=0;i<count;i++) {
137 v3_zero(&(itom[i].f));
140 virial=(&(itom[i].virial));
148 /* reset site energy */
153 /* get energy, force and virial of every atom */
155 /* first (and only) loop over atoms i */
156 for(i=0;i<count;i++) {
158 if(!(itom[i].attr&ATOM_ATTR_3BP))
161 link_cell_neighbour_index(moldyn,
162 (itom[i].r.x+moldyn->dim.x/2)/lc->x,
163 (itom[i].r.y+moldyn->dim.y/2)/lc->y,
164 (itom[i].r.z+moldyn->dim.z/2)/lc->z,
169 /* copy the neighbour lists */
173 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
179 /* loop over atoms j */
186 while(neighbour_i[j][p]!=-1) {
188 jtom=&(itom[neighbour_i[j][p]]);
196 p=lc->subcell->list[p];
198 this=&(neighbour_i[j]);
201 if(this->start==NULL)
206 jtom=this->current->data;
212 if(!(jtom->attr&ATOM_ATTR_3BP))
219 /* j1 func here ... */
220 /* albe 3 body potential function (first ij loop) */
226 * set ij depending values
229 if(brand_i==jtom->brand) {
230 S2=params->S2[brand_i];
237 v3_sub(&dist_ij,&(jtom->r),&(ai->r));
238 if(bc_ij) check_per_bound(moldyn,&dist_ij);
239 d_ij2=v3_absolute_square(&dist_ij);
241 /* if d_ij2 > S2 => no force & potential energy contribution */
248 /* reset k counter for first k loop */
251 /* first loop over atoms k */
258 while(neighbour_i[k][q]!=-1) {
260 ktom=&(itom[neighbour_i[k][q]]);
268 q=lc->subcell->list[q];
270 that=&(neighbour_i2[k]);
273 if(that->start==NULL)
277 ktom=that->current->data;
280 if(!(ktom->attr&ATOM_ATTR_3BP))
289 /* k1 func here ... */
290 /* albe 3 body potential function (first k loop) */
292 if(kcount>ALBE_MAXN) {
293 printf("FATAL: neighbours = %d\n",kcount);
294 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
298 if(brand_i==ktom->brand) {
299 Rk=params->R[brand_i];
300 Sk=params->S[brand_i];
301 Sk2=params->S2[brand_i];
302 /* albe needs i,k depending c,d,h and gamma values */
303 gamma_i=params->gamma[brand_i];
304 c_i=params->c[brand_i];
305 d_i=params->d[brand_i];
306 h_i=params->h[brand_i];
307 ci2=params->c2[brand_i];
308 di2=params->d2[brand_i];
309 ci2di2=params->c2d2[brand_i];
315 /* albe needs i,k depending c,d,h and gamma values */
316 gamma_i=params->gamma_m;
320 ci2=params->c2_mixed;
321 di2=params->d2_mixed;
322 ci2di2=params->c2d2_m;
326 v3_sub(&dist_ik,&(ktom->r),&(ai->r));
327 if(bc_ik) check_per_bound(moldyn,&dist_ik);
328 d_ik2=v3_absolute_square(&dist_ik);
330 /* store data for second k loop */
331 exchange->dist_ik[kcount]=dist_ik;
332 exchange->d_ik2[kcount]=d_ik2;
334 /* return if not within cutoff */
344 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
347 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
348 d2_h_cos2=exchange->di2+(h_cos*h_cos);
349 frac=exchange->ci2/d2_h_cos2;
350 g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
351 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
354 h_cos=h_i+cos_theta; // + in albe formalism
355 d2_h_cos2=di2+(h_cos*h_cos);
357 g=gamma_i*(1.0+ci2di2-frac);
358 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
360 /* zeta sum += f_c_ik * g_ijk */
368 arg=M_PI*(d_ik-Rk)/s_r;
369 f_c_ik=0.5+0.5*cos(arg);
370 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
374 /* store even more data for second k loop */
375 exchange->g[kcount]=g;
376 exchange->dg[kcount]=dg;
377 exchange->d_ik[kcount]=d_ik;
378 exchange->cos_theta[kcount]=cos_theta;
379 exchange->f_c_ik[kcount]=f_c_ik;
380 exchange->df_c_ik[kcount]=df_c_ik;
382 /* increase k counter */
390 } while(list_next_f(that)!=\
396 /* j2 func here ... */
399 if(brand_i==jtom->brand) {
400 S=params->S[brand_i];
401 R=params->R[brand_i];
402 B=params->B[brand_i];
403 A=params->A[brand_i];
404 r0=params->r0[brand_i];
405 mu=params->mu[brand_i];
406 lambda=params->lambda[brand_i];
415 lambda=params->lambda_m;
425 arg=M_PI*(d_ij-R)/s_r;
426 f_c=0.5+0.5*cos(arg);
427 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
431 f_a=-B*exp(-mu*(d_ij-r0));
435 f_r=A*exp(-lambda*(d_ij-r0));
436 df_r=lambda*f_r/d_ij;
444 b=1.0/sqrt(1.0+zeta_ij);
445 db=-0.5*b/(1.0+zeta_ij);
448 /* force contribution for atom i */
449 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
450 v3_scale(&force,&(dist_ij),scale);
451 v3_add(&(ai->f),&(ai->f),&force);
453 /* force contribution for atom j */
454 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
455 v3_add(&(jtom->f),&(jtom->f),&force);
458 virial_calc(ai,&force,&(dist_ij));
461 if(moldyn->time>DSTART&&moldyn->time<DEND) {
462 if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
463 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
464 printf(" adding %f %f %f\n",force.x,force.y,force.z);
465 if(ai==&(moldyn->atom[0]))
466 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
467 if(jtom==&(moldyn->atom[0]))
468 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
469 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
471 printf(" %f %f %f\n",zeta_ij,.0,.0);
476 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
477 pre_dzeta=0.5*f_a*f_c*db;
479 /* energy contribution */
480 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
481 moldyn->energy+=energy;
484 /* reset k counter for second k loop */
488 /* second loop over atoms k */
495 while(neighbour_i[k][q]!=-1) {
497 ktom=&(itom[neighbour_i[k][q]]);
505 q=lc->subcell->list[q];
507 that=&(neighbour_i2[k]);
510 if(that->start==NULL)
514 ktom=that->current->data;
517 if(!(ktom->attr&ATOM_ATTR_3BP))
527 /* k2 func here ... */
528 /* albe 3 body potential function (second k loop) */
531 printf("FATAL: neighbours!\n");
534 d_ik2=exchange->d_ik2[kcount];
536 if(brand_i==ktom->brand)
537 Sk2=params->S2[brand_i];
541 /* return if d_ik > S */
548 dist_ik=exchange->dist_ik[kcount];
549 d_ik=exchange->d_ik[kcount];
551 /* f_c_ik, df_c_ik */
552 f_c_ik=exchange->f_c_ik[kcount];
553 df_c_ik=exchange->df_c_ik[kcount];
555 /* g, dg, cos_theta */
556 g=exchange->g[kcount];
557 dg=exchange->dg[kcount];
558 cos_theta=exchange->cos_theta[kcount];
560 /* cos_theta derivatives wrt j,k */
561 dijdik_inv=1.0/(d_ij*d_ik);
562 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
563 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
564 v3_add(&dcosdrj,&dcosdrj,&tmp);
565 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
566 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
567 v3_add(&dcosdrk,&dcosdrk,&tmp);
569 /* f_c_ik * dg, df_c_ik * g */
573 /* derivative wrt j */
574 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
576 /* force contribution */
577 v3_add(&(jtom->f),&(jtom->f),&force);
580 if(moldyn->time>DSTART&&moldyn->time<DEND) {
581 if(jtom==&(moldyn->atom[DATOM])) {
582 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
583 printf(" adding %f %f %f\n",force.x,force.y,force.z);
584 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
585 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
586 printf(" d ij ik = %f %f\n",d_ij,d_ik);
592 virial_calc(ai,&force,&dist_ij);
594 /* force contribution to atom i */
595 v3_scale(&force,&force,-1.0);
596 v3_add(&(ai->f),&(ai->f),&force);
598 /* derivative wrt k */
599 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
600 v3_scale(&tmp,&dcosdrk,fcdg);
601 v3_add(&force,&force,&tmp);
602 v3_scale(&force,&force,pre_dzeta);
604 /* force contribution */
605 v3_add(&(ktom->f),&(ktom->f),&force);
608 if(moldyn->time>DSTART&&moldyn->time<DEND) {
609 if(ktom==&(moldyn->atom[DATOM])) {
610 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
611 printf(" adding %f %f %f\n",force.x,force.y,force.z);
612 printf(" total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
613 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
614 printf(" d ij ik = %f %f\n",d_ij,d_ik);
620 virial_calc(ai,&force,&dist_ik);
622 /* force contribution to atom i */
623 v3_scale(&force,&force,-1.0);
624 v3_add(&(ai->f),&(ai->f),&force);
626 /* increase k counter */
636 } while(list_next_f(that)!=\
647 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
662 //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
663 if(moldyn->time>DSTART&&moldyn->time<DEND) {
665 printf(" x: %0.40f\n",moldyn->atom[DATOM].f.x);
666 printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y);
667 printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z);
671 /* some postprocessing */
673 #pragma omp parallel for
675 for(i=0;i<count;i++) {
676 /* calculate global virial */
677 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
678 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
679 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
680 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
681 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
682 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
684 /* check forces regarding the given timestep */
685 if(v3_norm(&(itom[i].f))>\
686 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
687 printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
698 typedef struct s_pft_data {
703 void *potential_force_thread(void *ptr) {
705 t_pft_data *pft_data;
710 t_atom *itom,*jtom,*ktom;
713 int *neighbour_i[27];
719 t_list neighbour_i[27];
720 t_list neighbour_i2[27];
730 t_albe_mult_params *params;
731 t_albe_exchange *exchange;
745 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
746 double f_c_ik,df_c_ik;
749 double f_a,df_a,b,db,f_c,df_c;
758 double dijdik_inv,fcdg,dfcg;
759 t_3dvec dcosdrj,dcosdrk;
772 moldyn=pft_data->moldyn;
780 params=moldyn->pot_params;
783 /* get energy, force and virial of every atom */
785 /* first (and only) loop over atoms i */
786 for(i=0;i<count;i++) {
788 if(!(itom[i].attr&ATOM_ATTR_3BP))
791 link_cell_neighbour_index(moldyn,
792 (itom[i].r.x+moldyn->dim.x/2)/lc->x,
793 (itom[i].r.y+moldyn->dim.y/2)/lc->y,
794 (itom[i].r.z+moldyn->dim.z/2)/lc->z,
799 /* copy the neighbour lists */
803 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
809 /* loop over atoms j */
816 while(neighbour_i[j][p]!=-1) {
818 jtom=&(itom[neighbour_i[j][p]]);
826 p=lc->subcell->list[p];
828 this=&(neighbour_i[j]);
831 if(this->start==NULL)
836 jtom=this->current->data;
842 if(!(jtom->attr&ATOM_ATTR_3BP))
849 /* j1 func here ... */
850 /* albe 3 body potential function (first ij loop) */
856 * set ij depending values
859 if(brand_i==jtom->brand) {
860 S2=params->S2[brand_i];
867 v3_sub(&dist_ij,&(jtom->r),&(ai->r));
868 if(bc_ij) check_per_bound(moldyn,&dist_ij);
869 d_ij2=v3_absolute_square(&dist_ij);
871 /* if d_ij2 > S2 => no force & potential energy contribution */
878 /* reset k counter for first k loop */
881 /* first loop over atoms k */
888 while(neighbour_i[k][q]!=-1) {
890 ktom=&(itom[neighbour_i[k][q]]);
898 q=lc->subcell->list[q];
900 that=&(neighbour_i2[k]);
903 if(that->start==NULL)
907 ktom=that->current->data;
910 if(!(ktom->attr&ATOM_ATTR_3BP))
919 /* k1 func here ... */
920 /* albe 3 body potential function (first k loop) */
922 if(kcount>ALBE_MAXN) {
923 printf("FATAL: neighbours = %d\n",kcount);
924 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
928 if(brand_i==ktom->brand) {
929 Rk=params->R[brand_i];
930 Sk=params->S[brand_i];
931 Sk2=params->S2[brand_i];
932 /* albe needs i,k depending c,d,h and gamma values */
933 gamma_i=params->gamma[brand_i];
934 c_i=params->c[brand_i];
935 d_i=params->d[brand_i];
936 h_i=params->h[brand_i];
937 ci2=params->c2[brand_i];
938 di2=params->d2[brand_i];
939 ci2di2=params->c2d2[brand_i];
945 /* albe needs i,k depending c,d,h and gamma values */
946 gamma_i=params->gamma_m;
950 ci2=params->c2_mixed;
951 di2=params->d2_mixed;
952 ci2di2=params->c2d2_m;
956 v3_sub(&dist_ik,&(ktom->r),&(ai->r));
957 if(bc_ik) check_per_bound(moldyn,&dist_ik);
958 d_ik2=v3_absolute_square(&dist_ik);
960 /* store data for second k loop */
961 exchange->dist_ik[kcount]=dist_ik;
962 exchange->d_ik2[kcount]=d_ik2;
964 /* return if not within cutoff */
974 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
977 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
978 d2_h_cos2=exchange->di2+(h_cos*h_cos);
979 frac=exchange->ci2/d2_h_cos2;
980 g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
981 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
984 h_cos=h_i+cos_theta; // + in albe formalism
985 d2_h_cos2=di2+(h_cos*h_cos);
987 g=gamma_i*(1.0+ci2di2-frac);
988 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
990 /* zeta sum += f_c_ik * g_ijk */
998 arg=M_PI*(d_ik-Rk)/s_r;
999 f_c_ik=0.5+0.5*cos(arg);
1000 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
1004 /* store even more data for second k loop */
1005 exchange->g[kcount]=g;
1006 exchange->dg[kcount]=dg;
1007 exchange->d_ik[kcount]=d_ik;
1008 exchange->cos_theta[kcount]=cos_theta;
1009 exchange->f_c_ik[kcount]=f_c_ik;
1010 exchange->df_c_ik[kcount]=df_c_ik;
1012 /* increase k counter */
1020 } while(list_next_f(that)!=\
1026 /* j2 func here ... */
1029 if(brand_i==jtom->brand) {
1030 S=params->S[brand_i];
1031 R=params->R[brand_i];
1032 B=params->B[brand_i];
1033 A=params->A[brand_i];
1034 r0=params->r0[brand_i];
1035 mu=params->mu[brand_i];
1036 lambda=params->lambda[brand_i];
1043 r0=params->r0_mixed;
1045 lambda=params->lambda_m;
1055 arg=M_PI*(d_ij-R)/s_r;
1056 f_c=0.5+0.5*cos(arg);
1057 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
1061 f_a=-B*exp(-mu*(d_ij-r0));
1065 f_r=A*exp(-lambda*(d_ij-r0));
1066 df_r=lambda*f_r/d_ij;
1074 b=1.0/sqrt(1.0+zeta_ij);
1075 db=-0.5*b/(1.0+zeta_ij);
1078 /* force contribution for atom i */
1079 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
1080 v3_scale(&force,&(dist_ij),scale);
1081 v3_add(&(ai->f),&(ai->f),&force);
1083 /* force contribution for atom j */
1084 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
1085 v3_add(&(jtom->f),&(jtom->f),&force);
1088 virial_calc(ai,&force,&(dist_ij));
1091 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1092 if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
1093 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
1094 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1095 if(ai==&(moldyn->atom[0]))
1096 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
1097 if(jtom==&(moldyn->atom[0]))
1098 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1099 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
1101 printf(" %f %f %f\n",zeta_ij,.0,.0);
1106 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
1107 pre_dzeta=0.5*f_a*f_c*db;
1109 /* energy contribution */
1110 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
1111 moldyn->energy+=energy;
1114 /* reset k counter for second k loop */
1118 /* second loop over atoms k */
1125 while(neighbour_i[k][q]!=-1) {
1127 ktom=&(itom[neighbour_i[k][q]]);
1135 q=lc->subcell->list[q];
1137 that=&(neighbour_i2[k]);
1140 if(that->start==NULL)
1144 ktom=that->current->data;
1147 if(!(ktom->attr&ATOM_ATTR_3BP))
1153 if(ktom==&(itom[i]))
1157 /* k2 func here ... */
1158 /* albe 3 body potential function (second k loop) */
1160 if(kcount>ALBE_MAXN)
1161 printf("FATAL: neighbours!\n");
1164 d_ik2=exchange->d_ik2[kcount];
1166 if(brand_i==ktom->brand)
1167 Sk2=params->S2[brand_i];
1169 Sk2=params->S2mixed;
1171 /* return if d_ik > S */
1178 dist_ik=exchange->dist_ik[kcount];
1179 d_ik=exchange->d_ik[kcount];
1181 /* f_c_ik, df_c_ik */
1182 f_c_ik=exchange->f_c_ik[kcount];
1183 df_c_ik=exchange->df_c_ik[kcount];
1185 /* g, dg, cos_theta */
1186 g=exchange->g[kcount];
1187 dg=exchange->dg[kcount];
1188 cos_theta=exchange->cos_theta[kcount];
1190 /* cos_theta derivatives wrt j,k */
1191 dijdik_inv=1.0/(d_ij*d_ik);
1192 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
1193 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
1194 v3_add(&dcosdrj,&dcosdrj,&tmp);
1195 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
1196 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
1197 v3_add(&dcosdrk,&dcosdrk,&tmp);
1199 /* f_c_ik * dg, df_c_ik * g */
1203 /* derivative wrt j */
1204 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
1206 /* force contribution */
1207 v3_add(&(jtom->f),&(jtom->f),&force);
1210 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1211 if(jtom==&(moldyn->atom[DATOM])) {
1212 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1213 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1214 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1215 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1216 printf(" d ij ik = %f %f\n",d_ij,d_ik);
1222 virial_calc(ai,&force,&dist_ij);
1224 /* force contribution to atom i */
1225 v3_scale(&force,&force,-1.0);
1226 v3_add(&(ai->f),&(ai->f),&force);
1228 /* derivative wrt k */
1229 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
1230 v3_scale(&tmp,&dcosdrk,fcdg);
1231 v3_add(&force,&force,&tmp);
1232 v3_scale(&force,&force,pre_dzeta);
1234 /* force contribution */
1235 v3_add(&(ktom->f),&(ktom->f),&force);
1238 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1239 if(ktom==&(moldyn->atom[DATOM])) {
1240 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1241 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1242 printf(" total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
1243 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1244 printf(" d ij ik = %f %f\n",d_ij,d_ik);
1250 virial_calc(ai,&force,&dist_ik);
1252 /* force contribution to atom i */
1253 v3_scale(&force,&force,-1.0);
1254 v3_add(&(ai->f),&(ai->f),&force);
1256 /* increase k counter */
1266 } while(list_next_f(that)!=\
1277 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1292 //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
1293 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1295 printf(" x: %0.40f\n",moldyn->atom[DATOM].f.x);
1296 printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y);
1297 printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z);
1301 /* some postprocessing */
1303 #pragma omp parallel for
1305 for(i=0;i<count;i++) {
1306 /* calculate global virial */
1307 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
1308 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
1309 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
1310 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
1311 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
1312 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
1314 /* check forces regarding the given timestep */
1315 if(v3_norm(&(itom[i].f))>\
1316 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
1317 printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
1324 int albe_potential_force_calc(t_moldyn *moldyn) {
1327 t_pft_data *pft_data;
1329 pthread_t *pft_thread;
1333 count=moldyn->count;
1339 /* reset global virial */
1340 memset(&(moldyn->gvir),0,sizeof(t_virial));
1342 /* reset force, site energy and virial of every atom */
1343 for(i=0;i<count;i++) {
1346 v3_zero(&(itom[i].f));
1349 virial=&(itom[i].virial);
1357 /* reset site energy */
1362 /* alloc thread memory */
1363 pft_thread=malloc(count*sizeof(pthread_t));
1364 if(pft_thread==NULL) {
1365 perror("[albe fast] alloc thread mem");
1368 pft_data=malloc(count*sizeof(t_pft_data));
1369 if(pft_data==NULL) {
1370 perror("[albe fast] alloc thread mem");
1375 for(i=0;i<count;i++) {
1377 /* prepare thread data */
1378 pft_data[i].moldyn=moldyn;
1381 ret=pthread_create(&(pft_thread[i]),NULL,
1382 potential_force_thread,&(pft_data[i]));
1384 perror("[albe fast] pf thread create");
1390 for(i=0;i<count;i++) {
1391 ret=pthread_join(pft_thread[i],NULL);
1393 perror("[albe fast] pf thread join");