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 albe_v_calc(ai,&force,&(dist_ij));
465 //virial_calc(ai,&force,&(dist_ij));
468 if(moldyn->time>DSTART&&moldyn->time<DEND) {
469 if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
470 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
471 printf(" adding %f %f %f\n",force.x,force.y,force.z);
472 if(ai==&(moldyn->atom[0]))
473 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
474 if(jtom==&(moldyn->atom[0]))
475 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
476 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
478 printf(" %f %f %f\n",zeta_ij,.0,.0);
483 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
484 pre_dzeta=0.5*f_a*f_c*db;
486 /* energy contribution */
487 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
488 moldyn->energy+=energy;
491 /* reset k counter for second k loop */
495 /* second loop over atoms k */
502 while(neighbour_i[k][q]!=-1) {
504 ktom=&(itom[neighbour_i[k][q]]);
512 q=lc->subcell->list[q];
514 that=&(neighbour_i2[k]);
517 if(that->start==NULL)
521 ktom=that->current->data;
524 if(!(ktom->attr&ATOM_ATTR_3BP))
534 /* k2 func here ... */
535 /* albe 3 body potential function (second k loop) */
538 printf("FATAL: neighbours!\n");
541 d_ik2=exchange->d_ik2[kcount];
543 if(brand_i==ktom->brand)
544 Sk2=params->S2[brand_i];
548 /* return if d_ik > S */
555 dist_ik=exchange->dist_ik[kcount];
556 d_ik=exchange->d_ik[kcount];
558 /* f_c_ik, df_c_ik */
559 f_c_ik=exchange->f_c_ik[kcount];
560 df_c_ik=exchange->df_c_ik[kcount];
562 /* g, dg, cos_theta */
563 g=exchange->g[kcount];
564 dg=exchange->dg[kcount];
565 cos_theta=exchange->cos_theta[kcount];
567 /* cos_theta derivatives wrt j,k */
568 dijdik_inv=1.0/(d_ij*d_ik);
569 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
570 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
571 v3_add(&dcosdrj,&dcosdrj,&tmp);
572 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
573 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
574 v3_add(&dcosdrk,&dcosdrk,&tmp);
576 /* f_c_ik * dg, df_c_ik * g */
580 /* derivative wrt j */
581 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
583 /* force contribution */
584 v3_add(&(jtom->f),&(jtom->f),&force);
587 if(moldyn->time>DSTART&&moldyn->time<DEND) {
588 if(jtom==&(moldyn->atom[DATOM])) {
589 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
590 printf(" adding %f %f %f\n",force.x,force.y,force.z);
591 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
592 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
593 printf(" d ij ik = %f %f\n",d_ij,d_ik);
599 albe_v_calc(ai,&force,&dist_ij);
600 //virial_calc(ai,&force,&dist_ij);
602 /* force contribution to atom i */
603 v3_scale(&force,&force,-1.0);
604 v3_add(&(ai->f),&(ai->f),&force);
606 /* derivative wrt k */
607 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
608 v3_scale(&tmp,&dcosdrk,fcdg);
609 v3_add(&force,&force,&tmp);
610 v3_scale(&force,&force,pre_dzeta);
612 /* force contribution */
613 v3_add(&(ktom->f),&(ktom->f),&force);
616 if(moldyn->time>DSTART&&moldyn->time<DEND) {
617 if(ktom==&(moldyn->atom[DATOM])) {
618 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
619 printf(" adding %f %f %f\n",force.x,force.y,force.z);
620 printf(" total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
621 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
622 printf(" d ij ik = %f %f\n",d_ij,d_ik);
628 albe_v_calc(ai,&force,&dist_ik);
629 //virial_calc(ai,&force,&dist_ik);
631 /* force contribution to atom i */
632 v3_scale(&force,&force,-1.0);
633 v3_add(&(ai->f),&(ai->f),&force);
635 /* increase k counter */
645 } while(list_next_f(that)!=\
656 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
671 //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
672 if(moldyn->time>DSTART&&moldyn->time<DEND) {
674 printf(" x: %0.40f\n",moldyn->atom[DATOM].f.x);
675 printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y);
676 printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z);
680 /* some postprocessing */
682 #pragma omp parallel for
684 for(i=0;i<count;i++) {
685 /* calculate global virial */
686 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
687 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
688 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
689 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
690 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
691 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
693 /* check forces regarding the given timestep */
694 if(v3_norm(&(itom[i].f))>\
695 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
696 printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
707 typedef struct s_pft_data {
712 void *potential_force_thread(void *ptr) {
714 t_pft_data *pft_data;
719 t_atom *itom,*jtom,*ktom;
722 int *neighbour_i[27];
728 t_list neighbour_i[27];
729 t_list neighbour_i2[27];
739 t_albe_mult_params *params;
740 t_albe_exchange *exchange;
754 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
755 double f_c_ik,df_c_ik;
758 double f_a,df_a,b,db,f_c,df_c;
767 double dijdik_inv,fcdg,dfcg;
768 t_3dvec dcosdrj,dcosdrk;
781 moldyn=pft_data->moldyn;
789 params=moldyn->pot_params;
791 /* get energy, force and virial for atoms */
793 for(i=pft_data->start;i<pft_data->end;i++) {
795 if(!(itom[i].attr&ATOM_ATTR_3BP))
798 link_cell_neighbour_index(moldyn,
799 (itom[i].r.x+moldyn->dim.x/2)/lc->x,
800 (itom[i].r.y+moldyn->dim.y/2)/lc->y,
801 (itom[i].r.z+moldyn->dim.z/2)/lc->z,
806 /* copy the neighbour lists */
810 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
816 /* loop over atoms j */
823 while(neighbour_i[j][p]!=-1) {
825 jtom=&(itom[neighbour_i[j][p]]);
833 p=lc->subcell->list[p];
835 this=&(neighbour_i[j]);
838 if(this->start==NULL)
843 jtom=this->current->data;
849 if(!(jtom->attr&ATOM_ATTR_3BP))
856 /* j1 func here ... */
857 /* albe 3 body potential function (first ij loop) */
863 * set ij depending values
866 if(brand_i==jtom->brand) {
867 S2=params->S2[brand_i];
874 v3_sub(&dist_ij,&(jtom->r),&(ai->r));
875 if(bc_ij) check_per_bound(moldyn,&dist_ij);
876 d_ij2=v3_absolute_square(&dist_ij);
878 /* if d_ij2 > S2 => no force & potential energy contribution */
885 /* reset k counter for first k loop */
888 /* first loop over atoms k */
895 while(neighbour_i[k][q]!=-1) {
897 ktom=&(itom[neighbour_i[k][q]]);
905 q=lc->subcell->list[q];
907 that=&(neighbour_i2[k]);
910 if(that->start==NULL)
914 ktom=that->current->data;
917 if(!(ktom->attr&ATOM_ATTR_3BP))
926 /* k1 func here ... */
927 /* albe 3 body potential function (first k loop) */
929 if(kcount>ALBE_MAXN) {
930 printf("FATAL: neighbours = %d\n",kcount);
931 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
935 if(brand_i==ktom->brand) {
936 Rk=params->R[brand_i];
937 Sk=params->S[brand_i];
938 Sk2=params->S2[brand_i];
939 /* albe needs i,k depending c,d,h and gamma values */
940 gamma_i=params->gamma[brand_i];
941 c_i=params->c[brand_i];
942 d_i=params->d[brand_i];
943 h_i=params->h[brand_i];
944 ci2=params->c2[brand_i];
945 di2=params->d2[brand_i];
946 ci2di2=params->c2d2[brand_i];
952 /* albe needs i,k depending c,d,h and gamma values */
953 gamma_i=params->gamma_m;
957 ci2=params->c2_mixed;
958 di2=params->d2_mixed;
959 ci2di2=params->c2d2_m;
963 v3_sub(&dist_ik,&(ktom->r),&(ai->r));
964 if(bc_ik) check_per_bound(moldyn,&dist_ik);
965 d_ik2=v3_absolute_square(&dist_ik);
967 /* store data for second k loop */
968 exchange->dist_ik[kcount]=dist_ik;
969 exchange->d_ik2[kcount]=d_ik2;
971 /* return if not within cutoff */
981 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
984 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
985 d2_h_cos2=exchange->di2+(h_cos*h_cos);
986 frac=exchange->ci2/d2_h_cos2;
987 g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
988 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
991 h_cos=h_i+cos_theta; // + in albe formalism
992 d2_h_cos2=di2+(h_cos*h_cos);
994 g=gamma_i*(1.0+ci2di2-frac);
995 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
997 /* zeta sum += f_c_ik * g_ijk */
1005 arg=M_PI*(d_ik-Rk)/s_r;
1006 f_c_ik=0.5+0.5*cos(arg);
1007 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
1011 /* store even more data for second k loop */
1012 exchange->g[kcount]=g;
1013 exchange->dg[kcount]=dg;
1014 exchange->d_ik[kcount]=d_ik;
1015 exchange->cos_theta[kcount]=cos_theta;
1016 exchange->f_c_ik[kcount]=f_c_ik;
1017 exchange->df_c_ik[kcount]=df_c_ik;
1019 /* increase k counter */
1027 } while(list_next_f(that)!=\
1033 /* j2 func here ... */
1036 if(brand_i==jtom->brand) {
1037 S=params->S[brand_i];
1038 R=params->R[brand_i];
1039 B=params->B[brand_i];
1040 A=params->A[brand_i];
1041 r0=params->r0[brand_i];
1042 mu=params->mu[brand_i];
1043 lambda=params->lambda[brand_i];
1050 r0=params->r0_mixed;
1052 lambda=params->lambda_m;
1062 arg=M_PI*(d_ij-R)/s_r;
1063 f_c=0.5+0.5*cos(arg);
1064 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
1068 f_a=-B*exp(-mu*(d_ij-r0));
1072 f_r=A*exp(-lambda*(d_ij-r0));
1073 df_r=lambda*f_r/d_ij;
1081 b=1.0/sqrt(1.0+zeta_ij);
1082 db=-0.5*b/(1.0+zeta_ij);
1085 /* force contribution for atom i */
1086 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
1087 v3_scale(&force,&(dist_ij),scale);
1088 pthread_mutex_lock(&(amutex[ai->tag]));
1089 v3_add(&(ai->f),&(ai->f),&force);
1090 pthread_mutex_unlock(&(amutex[ai->tag]));
1092 /* force contribution for atom j */
1093 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
1094 pthread_mutex_lock(&(amutex[jtom->tag]));
1095 v3_add(&(jtom->f),&(jtom->f),&force);
1096 pthread_mutex_unlock(&(amutex[jtom->tag]));
1099 pthread_mutex_lock(&(amutex[ai->tag]));
1100 albe_v_calc(ai,&force,&(dist_ij));
1101 //virial_calc(ai,&force,&(dist_ij));
1102 pthread_mutex_unlock(&(amutex[ai->tag]));
1105 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1106 if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
1107 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
1108 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1109 if(ai==&(moldyn->atom[0]))
1110 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
1111 if(jtom==&(moldyn->atom[0]))
1112 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1113 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
1115 printf(" %f %f %f\n",zeta_ij,.0,.0);
1120 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
1121 pre_dzeta=0.5*f_a*f_c*db;
1123 /* energy contribution */
1124 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
1125 pthread_mutex_lock(&emutex);
1126 moldyn->energy+=energy;
1127 pthread_mutex_unlock(&emutex);
1128 pthread_mutex_lock(&(amutex[ai->tag]));
1130 pthread_mutex_unlock(&(amutex[ai->tag]));
1132 /* reset k counter for second k loop */
1136 /* second loop over atoms k */
1143 while(neighbour_i[k][q]!=-1) {
1145 ktom=&(itom[neighbour_i[k][q]]);
1153 q=lc->subcell->list[q];
1155 that=&(neighbour_i2[k]);
1158 if(that->start==NULL)
1162 ktom=that->current->data;
1165 if(!(ktom->attr&ATOM_ATTR_3BP))
1171 if(ktom==&(itom[i]))
1175 /* k2 func here ... */
1176 /* albe 3 body potential function (second k loop) */
1178 if(kcount>ALBE_MAXN)
1179 printf("FATAL: neighbours!\n");
1182 d_ik2=exchange->d_ik2[kcount];
1184 if(brand_i==ktom->brand)
1185 Sk2=params->S2[brand_i];
1187 Sk2=params->S2mixed;
1189 /* return if d_ik > S */
1196 dist_ik=exchange->dist_ik[kcount];
1197 d_ik=exchange->d_ik[kcount];
1199 /* f_c_ik, df_c_ik */
1200 f_c_ik=exchange->f_c_ik[kcount];
1201 df_c_ik=exchange->df_c_ik[kcount];
1203 /* g, dg, cos_theta */
1204 g=exchange->g[kcount];
1205 dg=exchange->dg[kcount];
1206 cos_theta=exchange->cos_theta[kcount];
1208 /* cos_theta derivatives wrt j,k */
1209 dijdik_inv=1.0/(d_ij*d_ik);
1210 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
1211 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
1212 v3_add(&dcosdrj,&dcosdrj,&tmp);
1213 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
1214 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
1215 v3_add(&dcosdrk,&dcosdrk,&tmp);
1217 /* f_c_ik * dg, df_c_ik * g */
1221 /* derivative wrt j */
1222 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
1224 /* force contribution */
1225 pthread_mutex_lock(&(amutex[jtom->tag]));
1226 v3_add(&(jtom->f),&(jtom->f),&force);
1227 pthread_mutex_unlock(&(amutex[jtom->tag]));
1230 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1231 if(jtom==&(moldyn->atom[DATOM])) {
1232 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1233 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1234 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1235 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1236 printf(" d ij ik = %f %f\n",d_ij,d_ik);
1242 pthread_mutex_lock(&(amutex[ai->tag]));
1243 albe_v_calc(ai,&force,&dist_ij);
1244 //virial_calc(ai,&force,&dist_ij);
1246 /* force contribution to atom i */
1247 v3_scale(&force,&force,-1.0);
1248 v3_add(&(ai->f),&(ai->f),&force);
1249 pthread_mutex_unlock(&(amutex[ai->tag]));
1251 /* derivative wrt k */
1252 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
1253 v3_scale(&tmp,&dcosdrk,fcdg);
1254 v3_add(&force,&force,&tmp);
1255 v3_scale(&force,&force,pre_dzeta);
1257 /* force contribution */
1258 pthread_mutex_lock(&(amutex[ktom->tag]));
1259 v3_add(&(ktom->f),&(ktom->f),&force);
1260 pthread_mutex_unlock(&(amutex[ktom->tag]));
1263 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1264 if(ktom==&(moldyn->atom[DATOM])) {
1265 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1266 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1267 printf(" total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
1268 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1269 printf(" d ij ik = %f %f\n",d_ij,d_ik);
1275 pthread_mutex_lock(&(amutex[ai->tag]));
1276 albe_v_calc(ai,&force,&dist_ik);
1277 //virial_calc(ai,&force,&dist_ik);
1279 /* force contribution to atom i */
1280 v3_scale(&force,&force,-1.0);
1281 v3_add(&(ai->f),&(ai->f),&force);
1282 pthread_mutex_unlock(&(amutex[ai->tag]));
1284 /* increase k counter */
1294 } while(list_next_f(that)!=\
1305 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1320 //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
1321 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1323 printf(" x: %0.40f\n",moldyn->atom[DATOM].f.x);
1324 printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y);
1325 printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z);
1334 int albe_potential_force_calc(t_moldyn *moldyn) {
1337 t_pft_data pft_data[MAX_THREADS];
1339 pthread_t pft_thread[MAX_THREADS];
1343 count=moldyn->count;
1349 /* reset global virial */
1350 memset(&(moldyn->gvir),0,sizeof(t_virial));
1352 /* reset force, site energy and virial of every atom */
1353 for(i=0;i<count;i++) {
1356 v3_zero(&(itom[i].f));
1359 virial=&(itom[i].virial);
1367 /* reset site energy */
1373 for(j=0;j<MAX_THREADS;j++) {
1375 /* prepare thread data */
1376 pft_data[j].moldyn=moldyn;
1377 pft_data[j].start=j*(count/MAX_THREADS);
1378 if(j==MAX_THREADS-1) {
1379 pft_data[j].end=count;
1382 pft_data[j].end=pft_data[j].start;
1383 pft_data[j].end+=count/MAX_THREADS;
1386 ret=pthread_create(&(pft_thread[j]),NULL,
1387 potential_force_thread,
1390 perror("[albe fast] pf thread create");
1396 for(j=0;j<MAX_THREADS;j++) {
1398 ret=pthread_join(pft_thread[j],NULL);
1400 perror("[albe fast] pf thread join");
1405 /* some postprocessing */
1406 for(i=0;i<count;i++) {
1407 /* calculate global virial */
1408 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
1409 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
1410 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
1411 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
1412 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
1413 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
1415 /* check forces regarding the given timestep */
1416 if(v3_norm(&(itom[i].f))>\
1417 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
1418 printf("[moldyn] WARNING: pfc (high force: atom %d)\n",