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 */
456 scale=-0.5*(f_c*(df_r-b*df_a)); // - in albe formalism
458 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
460 v3_scale(&force,&(dist_ij),scale);
461 v3_add(&(ai->f),&(ai->f),&force);
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);
468 albe_v_calc(ai,&force,&(dist_ij));
469 //virial_calc(ai,&force,&(dist_ij));
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),
482 printf(" %f %f %f\n",zeta_ij,.0,.0);
487 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
488 pre_dzeta=0.5*f_a*f_c*db;
490 /* energy contribution */
491 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
492 moldyn->energy+=energy;
495 /* reset k counter for second k loop */
499 /* second loop over atoms k */
506 while(neighbour_i[k][q]!=-1) {
508 ktom=&(itom[neighbour_i[k][q]]);
516 q=lc->subcell->list[q];
518 that=&(neighbour_i2[k]);
521 if(that->start==NULL)
525 ktom=that->current->data;
528 if(!(ktom->attr&ATOM_ATTR_3BP))
538 /* k2 func here ... */
539 /* albe 3 body potential function (second k loop) */
542 printf("FATAL: neighbours!\n");
545 d_ik2=exchange->d_ik2[kcount];
547 if(brand_i==ktom->brand)
548 Sk2=params->S2[brand_i];
552 /* return if d_ik > S */
559 dist_ik=exchange->dist_ik[kcount];
560 d_ik=exchange->d_ik[kcount];
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];
566 /* g, dg, cos_theta */
567 g=exchange->g[kcount];
568 dg=exchange->dg[kcount];
569 cos_theta=exchange->cos_theta[kcount];
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);
580 /* f_c_ik * dg, df_c_ik * g */
584 /* derivative wrt j */
585 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
587 /* force contribution */
588 v3_add(&(jtom->f),&(jtom->f),&force);
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);
603 albe_v_calc(ai,&force,&dist_ij);
604 //virial_calc(ai,&force,&dist_ij);
606 /* force contribution to atom i */
607 v3_scale(&force,&force,-1.0);
608 v3_add(&(ai->f),&(ai->f),&force);
610 /* derivative wrt k */
612 v3_scale(&tmp,&dcosdrk,fcdg);
613 v3_scale(&force,&tmp,pre_dzeta);
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);
621 /* force contribution */
622 v3_add(&(ktom->f),&(ktom->f),&force);
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);
637 albe_v_calc(ai,&force,&dist_ik);
638 //virial_calc(ai,&force,&dist_ik);
640 /* force contribution to atom i */
641 v3_scale(&force,&force,-1.0);
642 v3_add(&(ai->f),&(ai->f),&force);
644 /* increase k counter */
654 } while(list_next_f(that)!=\
665 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
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) {
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);
689 /* some postprocessing */
691 #pragma omp parallel for
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;
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",
716 typedef struct s_pft_data {
721 void *potential_force_thread(void *ptr) {
723 t_pft_data *pft_data;
728 t_atom *itom,*jtom,*ktom;
731 int *neighbour_i[27];
737 t_list neighbour_i[27];
738 t_list neighbour_i2[27];
748 t_albe_mult_params *params;
749 t_albe_exchange *exchange;
763 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
764 double f_c_ik,df_c_ik;
767 double f_a,df_a,b,db,f_c,df_c;
776 double dijdik_inv,fcdg,dfcg;
777 t_3dvec dcosdrj,dcosdrk;
790 moldyn=pft_data->moldyn;
798 params=moldyn->pot_params;
800 /* get energy, force and virial for atoms */
802 for(i=pft_data->start;i<pft_data->end;i++) {
804 if(!(itom[i].attr&ATOM_ATTR_3BP))
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,
814 /* copy the neighbour lists */
818 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
824 /* loop over atoms j */
831 while(neighbour_i[j][p]!=-1) {
833 jtom=&(itom[neighbour_i[j][p]]);
841 p=lc->subcell->list[p]; // thread safe!
843 this=&(neighbour_i[j]);
846 if(this->start==NULL)
851 jtom=this->current->data;
857 if(!(jtom->attr&ATOM_ATTR_3BP))
864 /* j1 func here ... */
865 /* albe 3 body potential function (first ij loop) */
871 * set ij depending values
874 if(brand_i==jtom->brand) {
875 S2=params->S2[brand_i];
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);
886 /* if d_ij2 > S2 => no force & potential energy contribution */
893 /* reset k counter for first k loop */
896 /* first loop over atoms k */
903 while(neighbour_i[k][q]!=-1) {
905 ktom=&(itom[neighbour_i[k][q]]);
913 q=lc->subcell->list[q];
915 that=&(neighbour_i2[k]);
918 if(that->start==NULL)
922 ktom=that->current->data;
925 if(!(ktom->attr&ATOM_ATTR_3BP))
934 /* k1 func here ... */
935 /* albe 3 body potential function (first k loop) */
937 if(kcount>ALBE_MAXN) {
938 printf("FATAL: neighbours = %d\n",kcount);
939 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
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];
960 /* albe needs i,k depending c,d,h and gamma values */
961 gamma_i=params->gamma_m;
965 ci2=params->c2_mixed;
966 di2=params->d2_mixed;
967 ci2di2=params->c2d2_m;
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);
975 /* store data for second k loop */
976 exchange->dist_ik[kcount]=dist_ik;
977 exchange->d_ik2[kcount]=d_ik2;
979 /* return if not within cutoff */
989 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
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..
999 h_cos=h_i+cos_theta; // + in albe formalism
1000 d2_h_cos2=di2+(h_cos*h_cos);
1002 g=gamma_i*(1.0+ci2di2-frac);
1003 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
1005 /* zeta sum += f_c_ik * g_ijk */
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));
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;
1027 /* increase k counter */
1035 } while(list_next_f(that)!=\
1041 /* j2 func here ... */
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];
1058 r0=params->r0_mixed;
1060 lambda=params->lambda_m;
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));
1076 f_a=-B*exp(-mu*(d_ij-r0));
1080 f_r=A*exp(-lambda*(d_ij-r0));
1081 df_r=lambda*f_r/d_ij;
1089 b=1.0/sqrt(1.0+zeta_ij);
1090 db=-0.5*b/(1.0+zeta_ij);
1093 /* force contribution for atom i */
1095 scale=-0.5*(f_c*(df_r-b*df_a)); // - in albe formalism
1097 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
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");
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");
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");
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),
1133 printf(" %f %f %f\n",zeta_ij,.0,.0);
1138 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
1139 pre_dzeta=0.5*f_a*f_c*db;
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");
1151 if(pthread_mutex_unlock(&(amutex[ai->tag])))
1152 perror("[albe fast] mutex unlock (4)\n");
1154 /* reset k counter for second k loop */
1158 /* second loop over atoms k */
1165 while(neighbour_i[k][q]!=-1) {
1167 ktom=&(itom[neighbour_i[k][q]]);
1175 q=lc->subcell->list[q];
1177 that=&(neighbour_i2[k]);
1180 if(that->start==NULL)
1184 ktom=that->current->data;
1187 if(!(ktom->attr&ATOM_ATTR_3BP))
1193 if(ktom==&(itom[i]))
1197 /* k2 func here ... */
1198 /* albe 3 body potential function (second k loop) */
1200 if(kcount>ALBE_MAXN)
1201 printf("FATAL: neighbours!\n");
1204 d_ik2=exchange->d_ik2[kcount];
1206 if(brand_i==ktom->brand)
1207 Sk2=params->S2[brand_i];
1209 Sk2=params->S2mixed;
1211 /* return if d_ik > S */
1218 dist_ik=exchange->dist_ik[kcount];
1219 d_ik=exchange->d_ik[kcount];
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];
1225 /* g, dg, cos_theta */
1226 g=exchange->g[kcount];
1227 dg=exchange->dg[kcount];
1228 cos_theta=exchange->cos_theta[kcount];
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);
1239 /* f_c_ik * dg, df_c_ik * g */
1243 /* derivative wrt j */
1244 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
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");
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);
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);
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");
1277 /* derivative wrt k */
1279 v3_scale(&tmp,&dcosdrk,fcdg);
1280 v3_scale(&force,&tmp,pre_dzeta);
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);
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");
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);
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);
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");
1319 /* increase k counter */
1329 } while(list_next_f(that)!=\
1340 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
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) {
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);
1369 int albe_potential_force_calc(t_moldyn *moldyn) {
1372 t_pft_data pft_data[MAX_THREADS];
1374 pthread_t pft_thread[MAX_THREADS];
1378 count=moldyn->count;
1384 /* reset global virial */
1385 memset(&(moldyn->gvir),0,sizeof(t_virial));
1387 /* reset force, site energy and virial of every atom */
1388 for(i=0;i<count;i++) {
1391 v3_zero(&(itom[i].f));
1394 virial=&(itom[i].virial);
1402 /* reset site energy */
1408 for(j=0;j<MAX_THREADS;j++) {
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;
1417 pft_data[j].end=pft_data[j].start;
1418 pft_data[j].end+=count/MAX_THREADS;
1421 ret=pthread_create(&(pft_thread[j]),NULL,
1422 potential_force_thread,
1425 perror("[albe fast] pf thread create");
1431 for(j=0;j<MAX_THREADS;j++) {
1433 ret=pthread_join(pft_thread[j],NULL);
1435 perror("[albe fast] pf thread join");
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;
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",