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 link_cell_neighbour_index(moldyn,
808 (itom[i].r.x+moldyn->dim.x/2)/lc->x,
809 (itom[i].r.y+moldyn->dim.y/2)/lc->y,
810 (itom[i].r.z+moldyn->dim.z/2)/lc->z,
815 /* copy the neighbour lists */
819 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
825 /* loop over atoms j */
832 while(neighbour_i[j][p]!=-1) {
834 jtom=&(itom[neighbour_i[j][p]]);
842 p=lc->subcell->list[p];
844 this=&(neighbour_i[j]);
847 if(this->start==NULL)
852 jtom=this->current->data;
858 if(!(jtom->attr&ATOM_ATTR_3BP))
865 /* j1 func here ... */
866 /* albe 3 body potential function (first ij loop) */
872 * set ij depending values
875 if(brand_i==jtom->brand) {
876 S2=params->S2[brand_i];
883 v3_sub(&dist_ij,&(jtom->r),&(ai->r));
884 if(bc_ij) check_per_bound(moldyn,&dist_ij);
885 d_ij2=v3_absolute_square(&dist_ij);
887 /* if d_ij2 > S2 => no force & potential energy contribution */
894 /* reset k counter for first k loop */
897 /* first loop over atoms k */
904 while(neighbour_i[k][q]!=-1) {
906 ktom=&(itom[neighbour_i[k][q]]);
914 q=lc->subcell->list[q];
916 that=&(neighbour_i2[k]);
919 if(that->start==NULL)
923 ktom=that->current->data;
926 if(!(ktom->attr&ATOM_ATTR_3BP))
935 /* k1 func here ... */
936 /* albe 3 body potential function (first k loop) */
938 if(kcount>ALBE_MAXN) {
939 printf("FATAL: neighbours = %d\n",kcount);
940 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
944 if(brand_i==ktom->brand) {
945 Rk=params->R[brand_i];
946 Sk=params->S[brand_i];
947 Sk2=params->S2[brand_i];
948 /* albe needs i,k depending c,d,h and gamma values */
949 gamma_i=params->gamma[brand_i];
950 c_i=params->c[brand_i];
951 d_i=params->d[brand_i];
952 h_i=params->h[brand_i];
953 ci2=params->c2[brand_i];
954 di2=params->d2[brand_i];
955 ci2di2=params->c2d2[brand_i];
961 /* albe needs i,k depending c,d,h and gamma values */
962 gamma_i=params->gamma_m;
966 ci2=params->c2_mixed;
967 di2=params->d2_mixed;
968 ci2di2=params->c2d2_m;
972 v3_sub(&dist_ik,&(ktom->r),&(ai->r));
973 if(bc_ik) check_per_bound(moldyn,&dist_ik);
974 d_ik2=v3_absolute_square(&dist_ik);
976 /* store data for second k loop */
977 exchange->dist_ik[kcount]=dist_ik;
978 exchange->d_ik2[kcount]=d_ik2;
980 /* return if not within cutoff */
990 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
993 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
994 d2_h_cos2=exchange->di2+(h_cos*h_cos);
995 frac=exchange->ci2/d2_h_cos2;
996 g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
997 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
1000 h_cos=h_i+cos_theta; // + in albe formalism
1001 d2_h_cos2=di2+(h_cos*h_cos);
1003 g=gamma_i*(1.0+ci2di2-frac);
1004 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
1006 /* zeta sum += f_c_ik * g_ijk */
1014 arg=M_PI*(d_ik-Rk)/s_r;
1015 f_c_ik=0.5+0.5*cos(arg);
1016 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
1020 /* store even more data for second k loop */
1021 exchange->g[kcount]=g;
1022 exchange->dg[kcount]=dg;
1023 exchange->d_ik[kcount]=d_ik;
1024 exchange->cos_theta[kcount]=cos_theta;
1025 exchange->f_c_ik[kcount]=f_c_ik;
1026 exchange->df_c_ik[kcount]=df_c_ik;
1028 /* increase k counter */
1036 } while(list_next_f(that)!=\
1042 /* j2 func here ... */
1045 if(brand_i==jtom->brand) {
1046 S=params->S[brand_i];
1047 R=params->R[brand_i];
1048 B=params->B[brand_i];
1049 A=params->A[brand_i];
1050 r0=params->r0[brand_i];
1051 mu=params->mu[brand_i];
1052 lambda=params->lambda[brand_i];
1059 r0=params->r0_mixed;
1061 lambda=params->lambda_m;
1071 arg=M_PI*(d_ij-R)/s_r;
1072 f_c=0.5+0.5*cos(arg);
1073 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
1077 f_a=-B*exp(-mu*(d_ij-r0));
1081 f_r=A*exp(-lambda*(d_ij-r0));
1082 df_r=lambda*f_r/d_ij;
1090 b=1.0/sqrt(1.0+zeta_ij);
1091 db=-0.5*b/(1.0+zeta_ij);
1094 /* force contribution for atom i */
1096 scale=-0.5*(f_c*(df_r-b*df_a)); // - in albe formalism
1098 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
1100 v3_scale(&force,&(dist_ij),scale);
1101 pthread_mutex_lock(&(amutex[ai->tag]));
1102 v3_add(&(ai->f),&(ai->f),&force);
1103 pthread_mutex_unlock(&(amutex[ai->tag]));
1105 /* force contribution for atom j */
1106 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
1107 pthread_mutex_lock(&(amutex[jtom->tag]));
1108 v3_add(&(jtom->f),&(jtom->f),&force);
1109 pthread_mutex_unlock(&(amutex[jtom->tag]));
1112 pthread_mutex_lock(&(amutex[ai->tag]));
1113 albe_v_calc(ai,&force,&(dist_ij));
1114 //virial_calc(ai,&force,&(dist_ij));
1115 pthread_mutex_unlock(&(amutex[ai->tag]));
1118 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1119 if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
1120 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
1121 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1122 if(ai==&(moldyn->atom[0]))
1123 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
1124 if(jtom==&(moldyn->atom[0]))
1125 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1126 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
1128 printf(" %f %f %f\n",zeta_ij,.0,.0);
1133 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
1134 pre_dzeta=0.5*f_a*f_c*db;
1136 /* energy contribution */
1137 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
1138 pthread_mutex_lock(&emutex);
1139 moldyn->energy+=energy;
1140 pthread_mutex_unlock(&emutex);
1141 pthread_mutex_lock(&(amutex[ai->tag]));
1143 pthread_mutex_unlock(&(amutex[ai->tag]));
1145 /* reset k counter for second k loop */
1149 /* second loop over atoms k */
1156 while(neighbour_i[k][q]!=-1) {
1158 ktom=&(itom[neighbour_i[k][q]]);
1166 q=lc->subcell->list[q];
1168 that=&(neighbour_i2[k]);
1171 if(that->start==NULL)
1175 ktom=that->current->data;
1178 if(!(ktom->attr&ATOM_ATTR_3BP))
1184 if(ktom==&(itom[i]))
1188 /* k2 func here ... */
1189 /* albe 3 body potential function (second k loop) */
1191 if(kcount>ALBE_MAXN)
1192 printf("FATAL: neighbours!\n");
1195 d_ik2=exchange->d_ik2[kcount];
1197 if(brand_i==ktom->brand)
1198 Sk2=params->S2[brand_i];
1200 Sk2=params->S2mixed;
1202 /* return if d_ik > S */
1209 dist_ik=exchange->dist_ik[kcount];
1210 d_ik=exchange->d_ik[kcount];
1212 /* f_c_ik, df_c_ik */
1213 f_c_ik=exchange->f_c_ik[kcount];
1214 df_c_ik=exchange->df_c_ik[kcount];
1216 /* g, dg, cos_theta */
1217 g=exchange->g[kcount];
1218 dg=exchange->dg[kcount];
1219 cos_theta=exchange->cos_theta[kcount];
1221 /* cos_theta derivatives wrt j,k */
1222 dijdik_inv=1.0/(d_ij*d_ik);
1223 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
1224 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
1225 v3_add(&dcosdrj,&dcosdrj,&tmp);
1226 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
1227 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
1228 v3_add(&dcosdrk,&dcosdrk,&tmp);
1230 /* f_c_ik * dg, df_c_ik * g */
1234 /* derivative wrt j */
1235 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
1237 /* force contribution */
1238 pthread_mutex_lock(&(amutex[jtom->tag]));
1239 v3_add(&(jtom->f),&(jtom->f),&force);
1240 pthread_mutex_unlock(&(amutex[jtom->tag]));
1243 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1244 if(jtom==&(moldyn->atom[DATOM])) {
1245 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1246 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1247 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
1248 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1249 printf(" d ij ik = %f %f\n",d_ij,d_ik);
1255 pthread_mutex_lock(&(amutex[ai->tag]));
1256 albe_v_calc(ai,&force,&dist_ij);
1257 //virial_calc(ai,&force,&dist_ij);
1259 /* force contribution to atom i */
1260 v3_scale(&force,&force,-1.0);
1261 v3_add(&(ai->f),&(ai->f),&force);
1262 pthread_mutex_unlock(&(amutex[ai->tag]));
1264 /* derivative wrt k */
1266 v3_scale(&tmp,&dcosdrk,fcdg);
1267 v3_scale(&force,&tmp,pre_dzeta);
1269 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
1270 v3_scale(&tmp,&dcosdrk,fcdg);
1271 v3_add(&force,&force,&tmp);
1272 v3_scale(&force,&force,pre_dzeta);
1275 /* force contribution */
1276 pthread_mutex_lock(&(amutex[ktom->tag]));
1277 v3_add(&(ktom->f),&(ktom->f),&force);
1278 pthread_mutex_unlock(&(amutex[ktom->tag]));
1281 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1282 if(ktom==&(moldyn->atom[DATOM])) {
1283 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
1284 printf(" adding %f %f %f\n",force.x,force.y,force.z);
1285 printf(" total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
1286 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
1287 printf(" d ij ik = %f %f\n",d_ij,d_ik);
1293 pthread_mutex_lock(&(amutex[ai->tag]));
1294 albe_v_calc(ai,&force,&dist_ik);
1295 //virial_calc(ai,&force,&dist_ik);
1297 /* force contribution to atom i */
1298 v3_scale(&force,&force,-1.0);
1299 v3_add(&(ai->f),&(ai->f),&force);
1300 pthread_mutex_unlock(&(amutex[ai->tag]));
1302 /* increase k counter */
1312 } while(list_next_f(that)!=\
1323 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
1338 //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
1339 if(moldyn->time>DSTART&&moldyn->time<DEND) {
1341 printf(" x: %0.40f\n",moldyn->atom[DATOM].f.x);
1342 printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y);
1343 printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z);
1352 int albe_potential_force_calc(t_moldyn *moldyn) {
1355 t_pft_data pft_data[MAX_THREADS];
1357 pthread_t pft_thread[MAX_THREADS];
1361 count=moldyn->count;
1367 /* reset global virial */
1368 memset(&(moldyn->gvir),0,sizeof(t_virial));
1370 /* reset force, site energy and virial of every atom */
1371 for(i=0;i<count;i++) {
1374 v3_zero(&(itom[i].f));
1377 virial=&(itom[i].virial);
1385 /* reset site energy */
1391 for(j=0;j<MAX_THREADS;j++) {
1393 /* prepare thread data */
1394 pft_data[j].moldyn=moldyn;
1395 pft_data[j].start=j*(count/MAX_THREADS);
1396 if(j==MAX_THREADS-1) {
1397 pft_data[j].end=count;
1400 pft_data[j].end=pft_data[j].start;
1401 pft_data[j].end+=count/MAX_THREADS;
1404 ret=pthread_create(&(pft_thread[j]),NULL,
1405 potential_force_thread,
1408 perror("[albe fast] pf thread create");
1414 for(j=0;j<MAX_THREADS;j++) {
1416 ret=pthread_join(pft_thread[j],NULL);
1418 perror("[albe fast] pf thread join");
1423 /* some postprocessing */
1424 for(i=0;i<count;i++) {
1425 /* calculate global virial */
1426 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
1427 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
1428 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
1429 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
1430 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
1431 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
1433 /* check forces regarding the given timestep */
1434 if(v3_norm(&(itom[i].f))>\
1435 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
1436 printf("[moldyn] WARNING: pfc (high force: atom %d)\n",