4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
12 #include <sys/types.h>
29 #include "../moldyn.h"
30 #include "../math/math.h"
34 extern pthread_mutex_t *amutex;
35 extern pthread_mutex_t emutex;
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
51 int albe_potential_force_calc(t_moldyn *moldyn) {
54 t_atom *itom,*jtom,*ktom;
64 t_list neighbour_i[27];
65 t_list neighbour_i2[27];
73 pthread_t kthread[27];
80 t_albe_mult_params *params;
81 t_albe_exchange *exchange;
95 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
96 double f_c_ik,df_c_ik;
99 double f_a,df_a,b,db,f_c,df_c;
108 double dijdik_inv,fcdg,dfcg;
109 t_3dvec dcosdrj,dcosdrk;
126 params=moldyn->pot_params;
127 exchange=&(params->exchange);
133 /* reset global virial */
134 memset(&(moldyn->gvir),0,sizeof(t_virial));
136 /* reset force, site energy and virial of every atom */
138 #pragma omp parallel for private(virial)
140 for(i=0;i<count;i++) {
143 v3_zero(&(itom[i].f));
146 virial=(&(itom[i].virial));
154 /* reset site energy */
159 /* get energy, force and virial of every atom */
161 /* first (and only) loop over atoms i */
162 for(i=0;i<count;i++) {
164 if(!(itom[i].attr&ATOM_ATTR_3BP))
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,
175 /* copy the neighbour lists */
179 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
185 /* loop over atoms j */
192 while(neighbour_i[j][p]!=-1) {
194 jtom=&(itom[neighbour_i[j][p]]);
202 p=lc->subcell->list[p];
204 this=&(neighbour_i[j]);
207 if(this->start==NULL)
212 jtom=this->current->data;
218 if(!(jtom->attr&ATOM_ATTR_3BP))
225 /* j1 func here ... */
226 /* albe 3 body potential function (first ij loop) */
232 * set ij depending values
235 if(brand_i==jtom->brand) {
236 S2=params->S2[brand_i];
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);
247 /* if d_ij2 > S2 => no force & potential energy contribution */
254 /* reset k counter for first k loop */
257 /* first loop over atoms k */
264 while(neighbour_i[k][q]!=-1) {
266 ktom=&(itom[neighbour_i[k][q]]);
274 q=lc->subcell->list[q];
276 that=&(neighbour_i2[k]);
279 if(that->start==NULL)
283 ktom=that->current->data;
286 if(!(ktom->attr&ATOM_ATTR_3BP))
295 /* k1 func here ... */
296 /* albe 3 body potential function (first k loop) */
298 if(kcount>ALBE_MAXN) {
299 printf("FATAL: neighbours = %d\n",kcount);
300 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
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];
321 /* albe needs i,k depending c,d,h and gamma values */
322 gamma_i=params->gamma_m;
326 ci2=params->c2_mixed;
327 di2=params->d2_mixed;
328 ci2di2=params->c2d2_m;
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);
336 /* store data for second k loop */
337 exchange->dist_ik[kcount]=dist_ik;
338 exchange->d_ik2[kcount]=d_ik2;
340 /* return if not within cutoff */
350 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
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..
360 h_cos=h_i+cos_theta; // + in albe formalism
361 d2_h_cos2=di2+(h_cos*h_cos);
363 g=gamma_i*(1.0+ci2di2-frac);
364 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
366 /* zeta sum += f_c_ik * g_ijk */
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));
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;
388 /* increase k counter */
396 } while(list_next_f(that)!=\
402 /* j2 func here ... */
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];
421 lambda=params->lambda_m;
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));
437 f_a=-B*exp(-mu*(d_ij-r0));
441 f_r=A*exp(-lambda*(d_ij-r0));
442 df_r=lambda*f_r/d_ij;
450 b=1.0/sqrt(1.0+zeta_ij);
451 db=-0.5*b/(1.0+zeta_ij);
454 /* force contribution for atom i */
455 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
456 v3_scale(&force,&(dist_ij),scale);
457 v3_add(&(ai->f),&(ai->f),&force);
459 /* force contribution for atom j */
460 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
461 v3_add(&(jtom->f),&(jtom->f),&force);
464 virial_calc(ai,&force,&(dist_ij));
467 if(moldyn->time>DSTART&&moldyn->time<DEND) {
468 if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
469 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
470 printf(" adding %f %f %f\n",force.x,force.y,force.z);
471 if(ai==&(moldyn->atom[0]))
472 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
473 if(jtom==&(moldyn->atom[0]))
474 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
475 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
477 printf(" %f %f %f\n",zeta_ij,.0,.0);
482 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
483 pre_dzeta=0.5*f_a*f_c*db;
485 /* energy contribution */
486 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
487 moldyn->energy+=energy;
490 /* reset k counter for second k loop */
494 /* second loop over atoms k */
501 while(neighbour_i[k][q]!=-1) {
503 ktom=&(itom[neighbour_i[k][q]]);
511 q=lc->subcell->list[q];
513 that=&(neighbour_i2[k]);
516 if(that->start==NULL)
520 ktom=that->current->data;
523 if(!(ktom->attr&ATOM_ATTR_3BP))
533 /* k2 func here ... */
534 /* albe 3 body potential function (second k loop) */
537 printf("FATAL: neighbours!\n");
540 d_ik2=exchange->d_ik2[kcount];
542 if(brand_i==ktom->brand)
543 Sk2=params->S2[brand_i];
547 /* return if d_ik > S */
554 dist_ik=exchange->dist_ik[kcount];
555 d_ik=exchange->d_ik[kcount];
557 /* f_c_ik, df_c_ik */
558 f_c_ik=exchange->f_c_ik[kcount];
559 df_c_ik=exchange->df_c_ik[kcount];
561 /* g, dg, cos_theta */
562 g=exchange->g[kcount];
563 dg=exchange->dg[kcount];
564 cos_theta=exchange->cos_theta[kcount];
566 /* cos_theta derivatives wrt j,k */
567 dijdik_inv=1.0/(d_ij*d_ik);
568 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
569 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
570 v3_add(&dcosdrj,&dcosdrj,&tmp);
571 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
572 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
573 v3_add(&dcosdrk,&dcosdrk,&tmp);
575 /* f_c_ik * dg, df_c_ik * g */
579 /* derivative wrt j */
580 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
582 /* force contribution */
583 v3_add(&(jtom->f),&(jtom->f),&force);
586 if(moldyn->time>DSTART&&moldyn->time<DEND) {
587 if(jtom==&(moldyn->atom[DATOM])) {
588 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
589 printf(" adding %f %f %f\n",force.x,force.y,force.z);
590 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
591 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
592 printf(" d ij ik = %f %f\n",d_ij,d_ik);
598 virial_calc(ai,&force,&dist_ij);
600 /* force contribution to atom i */
601 v3_scale(&force,&force,-1.0);
602 v3_add(&(ai->f),&(ai->f),&force);
604 /* derivative wrt k */
605 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
606 v3_scale(&tmp,&dcosdrk,fcdg);
607 v3_add(&force,&force,&tmp);
608 v3_scale(&force,&force,pre_dzeta);
610 /* force contribution */
611 v3_add(&(ktom->f),&(ktom->f),&force);
614 if(moldyn->time>DSTART&&moldyn->time<DEND) {
615 if(ktom==&(moldyn->atom[DATOM])) {
616 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
617 printf(" adding %f %f %f\n",force.x,force.y,force.z);
618 printf(" total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
619 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
620 printf(" d ij ik = %f %f\n",d_ij,d_ik);
626 virial_calc(ai,&force,&dist_ik);
628 /* force contribution to atom i */
629 v3_scale(&force,&force,-1.0);
630 v3_add(&(ai->f),&(ai->f),&force);
632 /* increase k counter */
642 } while(list_next_f(that)!=\
653 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
668 //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
669 if(moldyn->time>DSTART&&moldyn->time<DEND) {
671 printf(" x: %0.40f\n",moldyn->atom[DATOM].f.x);
672 printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y);
673 printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z);
677 /* some postprocessing */
679 #pragma omp parallel for
681 for(i=0;i<count;i++) {
682 /* calculate global virial */
683 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
684 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
685 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
686 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
687 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
688 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
690 /* check forces regarding the given timestep */
691 if(v3_norm(&(itom[i].f))>\
692 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
693 printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
704 typedef struct s_pft_data {
709 void *potential_force_thread(void *ptr) {
711 t_pft_data *pft_data;
716 t_atom *itom,*jtom,*ktom;
719 int *neighbour_i[27];
725 t_list neighbour_i[27];
726 t_list neighbour_i2[27];
736 t_albe_mult_params *params;
737 t_albe_exchange *exchange;
751 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
752 double f_c_ik,df_c_ik;
755 double f_a,df_a,b,db,f_c,df_c;
764 double dijdik_inv,fcdg,dfcg;
765 t_3dvec dcosdrj,dcosdrk;
778 moldyn=pft_data->moldyn;
786 params=moldyn->pot_params;
788 /* get energy, force and virial for atoms */
790 for(i=pft_data->start;i<pft_data->end;i++) {
792 if(!(itom[i].attr&ATOM_ATTR_3BP))
795 link_cell_neighbour_index(moldyn,
796 (itom[i].r.x+moldyn->dim.x/2)/lc->x,
797 (itom[i].r.y+moldyn->dim.y/2)/lc->y,
798 (itom[i].r.z+moldyn->dim.z/2)/lc->z,
803 /* copy the neighbour lists */
807 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
813 /* loop over atoms j */
820 while(neighbour_i[j][p]!=-1) {
822 jtom=&(itom[neighbour_i[j][p]]);
830 p=lc->subcell->list[p];
832 this=&(neighbour_i[j]);
835 if(this->start==NULL)
840 jtom=this->current->data;
846 if(!(jtom->attr&ATOM_ATTR_3BP))
853 /* j1 func here ... */
854 /* albe 3 body potential function (first ij loop) */
860 * set ij depending values
863 if(brand_i==jtom->brand) {
864 S2=params->S2[brand_i];
871 v3_sub(&dist_ij,&(jtom->r),&(ai->r));
872 if(bc_ij) check_per_bound(moldyn,&dist_ij);
873 d_ij2=v3_absolute_square(&dist_ij);
875 /* if d_ij2 > S2 => no force & potential energy contribution */
882 /* reset k counter for first k loop */
885 /* first loop over atoms k */
892 while(neighbour_i[k][q]!=-1) {
894 ktom=&(itom[neighbour_i[k][q]]);
902 q=lc->subcell->list[q];
904 that=&(neighbour_i2[k]);
907 if(that->start==NULL)
911 ktom=that->current->data;
914 if(!(ktom->attr&ATOM_ATTR_3BP))
923 /* k1 func here ... */
924 /* albe 3 body potential function (first k loop) */
926 if(kcount>ALBE_MAXN) {
927 printf("FATAL: neighbours = %d\n",kcount);
928 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
932 if(brand_i==ktom->brand) {
933 Rk=params->R[brand_i];
934 Sk=params->S[brand_i];
935 Sk2=params->S2[brand_i];
936 /* albe needs i,k depending c,d,h and gamma values */
937 gamma_i=params->gamma[brand_i];
938 c_i=params->c[brand_i];
939 d_i=params->d[brand_i];
940 h_i=params->h[brand_i];
941 ci2=params->c2[brand_i];
942 di2=params->d2[brand_i];
943 ci2di2=params->c2d2[brand_i];
949 /* albe needs i,k depending c,d,h and gamma values */
950 gamma_i=params->gamma_m;
954 ci2=params->c2_mixed;
955 di2=params->d2_mixed;
956 ci2di2=params->c2d2_m;
960 v3_sub(&dist_ik,&(ktom->r),&(ai->r));
961 if(bc_ik) check_per_bound(moldyn,&dist_ik);
962 d_ik2=v3_absolute_square(&dist_ik);
964 /* store data for second k loop */
965 exchange->dist_ik[kcount]=dist_ik;
966 exchange->d_ik2[kcount]=d_ik2;
968 /* return if not within cutoff */
978 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
981 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
982 d2_h_cos2=exchange->di2+(h_cos*h_cos);
983 frac=exchange->ci2/d2_h_cos2;
984 g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
985 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
988 h_cos=h_i+cos_theta; // + in albe formalism
989 d2_h_cos2=di2+(h_cos*h_cos);
991 g=gamma_i*(1.0+ci2di2-frac);
992 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
994 /* zeta sum += f_c_ik * g_ijk */
1002 arg=M_PI*(d_ik-Rk)/s_r;
1003 f_c_ik=0.5+0.5*cos(arg);
1004 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
1008 /* store even more data for second k loop */
1009 exchange->g[kcount]=g;
1010 exchange->dg[kcount]=dg;
1011 exchange->d_ik[kcount]=d_ik;
1012 exchange->cos_theta[kcount]=cos_theta;
1013 exchange->f_c_ik[kcount]=f_c_ik;
1014 exchange->df_c_ik[kcount]=df_c_ik;
1016 /* increase k counter */
1024 } while(list_next_f(that)!=\
1030 /* j2 func here ... */
1033 if(brand_i==jtom->brand) {
1034 S=params->S[brand_i];
1035 R=params->R[brand_i];
1036 B=params->B[brand_i];
1037 A=params->A[brand_i];
1038 r0=params->r0[brand_i];
1039 mu=params->mu[brand_i];
1040 lambda=params->lambda[brand_i];
1047 r0=params->r0_mixed;
1049 lambda=params->lambda_m;
1059 arg=M_PI*(d_ij-R)/s_r;
1060 f_c=0.5+0.5*cos(arg);
1061 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
1065 f_a=-B*exp(-mu*(d_ij-r0));
1069 f_r=A*exp(-lambda*(d_ij-r0));
1070 df_r=lambda*f_r/d_ij;
1078 b=1.0/sqrt(1.0+zeta_ij);
1079 db=-0.5*b/(1.0+zeta_ij);
1082 /* force contribution for atom i */
1083 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
1084 v3_scale(&force,&(dist_ij),scale);
1085 pthread_mutex_lock(&(amutex[ai->tag]));
1086 v3_add(&(ai->f),&(ai->f),&force);
1087 pthread_mutex_unlock(&(amutex[ai->tag]));
1089 /* force contribution for atom j */
1090 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
1091 pthread_mutex_lock(&(amutex[jtom->tag]));
1092 v3_add(&(jtom->f),&(jtom->f),&force);
1093 pthread_mutex_unlock(&(amutex[jtom->tag]));
1096 pthread_mutex_lock(&(amutex[ai->tag]));
1097 virial_calc(ai,&force,&(dist_ij));
1098 pthread_mutex_unlock(&(amutex[ai->tag]));
1101 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1102 if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
1103 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
1104 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1105 if(ai==&(moldyn->atom[0]))
1106 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
1107 if(jtom==&(moldyn->atom[0]))
1108 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1109 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
1111 printf(" %f %f %f\n",zeta_ij,.0,.0);
1116 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
1117 pre_dzeta=0.5*f_a*f_c*db;
1119 /* energy contribution */
1120 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
1121 pthread_mutex_lock(&emutex);
1122 moldyn->energy+=energy;
1123 pthread_mutex_unlock(&emutex);
1124 pthread_mutex_lock(&(amutex[ai->tag]));
1126 pthread_mutex_unlock(&(amutex[ai->tag]));
1128 /* reset k counter for second k loop */
1132 /* second loop over atoms k */
1139 while(neighbour_i[k][q]!=-1) {
1141 ktom=&(itom[neighbour_i[k][q]]);
1149 q=lc->subcell->list[q];
1151 that=&(neighbour_i2[k]);
1154 if(that->start==NULL)
1158 ktom=that->current->data;
1161 if(!(ktom->attr&ATOM_ATTR_3BP))
1167 if(ktom==&(itom[i]))
1171 /* k2 func here ... */
1172 /* albe 3 body potential function (second k loop) */
1174 if(kcount>ALBE_MAXN)
1175 printf("FATAL: neighbours!\n");
1178 d_ik2=exchange->d_ik2[kcount];
1180 if(brand_i==ktom->brand)
1181 Sk2=params->S2[brand_i];
1183 Sk2=params->S2mixed;
1185 /* return if d_ik > S */
1192 dist_ik=exchange->dist_ik[kcount];
1193 d_ik=exchange->d_ik[kcount];
1195 /* f_c_ik, df_c_ik */
1196 f_c_ik=exchange->f_c_ik[kcount];
1197 df_c_ik=exchange->df_c_ik[kcount];
1199 /* g, dg, cos_theta */
1200 g=exchange->g[kcount];
1201 dg=exchange->dg[kcount];
1202 cos_theta=exchange->cos_theta[kcount];
1204 /* cos_theta derivatives wrt j,k */
1205 dijdik_inv=1.0/(d_ij*d_ik);
1206 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
1207 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
1208 v3_add(&dcosdrj,&dcosdrj,&tmp);
1209 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
1210 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
1211 v3_add(&dcosdrk,&dcosdrk,&tmp);
1213 /* f_c_ik * dg, df_c_ik * g */
1217 /* derivative wrt j */
1218 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
1220 /* force contribution */
1221 pthread_mutex_lock(&(amutex[jtom->tag]));
1222 v3_add(&(jtom->f),&(jtom->f),&force);
1223 pthread_mutex_unlock(&(amutex[jtom->tag]));
1226 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1227 if(jtom==&(moldyn->atom[DATOM])) {
1228 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1229 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1230 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1231 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1232 printf(" d ij ik = %f %f\n",d_ij,d_ik);
1238 pthread_mutex_lock(&(amutex[ai->tag]));
1239 virial_calc(ai,&force,&dist_ij);
1241 /* force contribution to atom i */
1242 v3_scale(&force,&force,-1.0);
1243 v3_add(&(ai->f),&(ai->f),&force);
1244 pthread_mutex_unlock(&(amutex[ai->tag]));
1246 /* derivative wrt k */
1247 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
1248 v3_scale(&tmp,&dcosdrk,fcdg);
1249 v3_add(&force,&force,&tmp);
1250 v3_scale(&force,&force,pre_dzeta);
1252 /* force contribution */
1253 pthread_mutex_lock(&(amutex[ktom->tag]));
1254 v3_add(&(ktom->f),&(ktom->f),&force);
1255 pthread_mutex_unlock(&(amutex[ktom->tag]));
1258 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1259 if(ktom==&(moldyn->atom[DATOM])) {
1260 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1261 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1262 printf(" total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
1263 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1264 printf(" d ij ik = %f %f\n",d_ij,d_ik);
1270 pthread_mutex_lock(&(amutex[ai->tag]));
1271 virial_calc(ai,&force,&dist_ik);
1273 /* force contribution to atom i */
1274 v3_scale(&force,&force,-1.0);
1275 v3_add(&(ai->f),&(ai->f),&force);
1276 pthread_mutex_unlock(&(amutex[ai->tag]));
1278 /* increase k counter */
1288 } while(list_next_f(that)!=\
1299 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1314 //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
1315 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1317 printf(" x: %0.40f\n",moldyn->atom[DATOM].f.x);
1318 printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y);
1319 printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z);
1328 int albe_potential_force_calc(t_moldyn *moldyn) {
1331 t_pft_data pft_data[MAX_THREADS];
1333 pthread_t pft_thread[MAX_THREADS];
1337 count=moldyn->count;
1343 /* reset global virial */
1344 memset(&(moldyn->gvir),0,sizeof(t_virial));
1346 /* reset force, site energy and virial of every atom */
1347 for(i=0;i<count;i++) {
1350 v3_zero(&(itom[i].f));
1353 virial=&(itom[i].virial);
1361 /* reset site energy */
1367 for(j=0;j<MAX_THREADS;j++) {
1369 /* prepare thread data */
1370 pft_data[j].moldyn=moldyn;
1371 pft_data[j].start=j*count/MAX_THREADS;
1372 if(j==MAX_THREADS-1) {
1373 pft_data[j].end=count;
1376 pft_data[j].end=pft_data[j].start;
1377 pft_data[j].end+=count/MAX_THREADS;
1380 ret=pthread_create(&(pft_thread[j]),NULL,
1381 potential_force_thread,
1384 perror("[albe fast] pf thread create");
1390 for(j=0;j<MAX_THREADS;j++) {
1392 ret=pthread_join(pft_thread[j],NULL);
1394 perror("[albe fast] pf thread join");
1399 /* some postprocessing */
1400 for(i=0;i<count;i++) {
1401 /* calculate global virial */
1402 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
1403 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
1404 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
1405 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
1406 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
1407 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
1409 /* check forces regarding the given timestep */
1410 if(v3_norm(&(itom[i].f))>\
1411 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
1412 printf("[moldyn] WARNING: pfc (high force: atom %d)\n",