4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
12 #include <sys/types.h>
28 #include "../moldyn.h"
29 #include "../math/math.h"
33 typedef struct s_kdata {
43 #define albe_v_calc(a,f,d) a->virial.xx+=f->x*d->x; \
44 a->virial.yy+=f->y*d->y; \
45 a->virial.zz+=f->z*d->z; \
46 a->virial.xy+=f->x*d->y; \
47 a->virial.xz+=f->x*d->z; \
48 a->virial.yz+=f->y*d->z
52 void *k1_calc(void *ptr) {
54 /* albe 3 body potential function (first k loop) */
56 t_albe_mult_params *params;
57 t_albe_exchange *exchange;
58 unsigned char brand_i;
59 double Rk,Sk,Sk2,gamma_i,c_i,d_i,h_i,ci2,di2,ci2di2;
60 t_atom *ai,*jtom,*ktom;
63 if(kcount>ALBE_MAXN) {
64 printf("FATAL: neighbours = %d\n",kcount);
65 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
69 if(brand_i==ktom->brand) {
70 Rk=params->R[brand_i];
71 Sk=params->S[brand_i];
72 Sk2=params->S2[brand_i];
73 /* albe needs i,k depending c,d,h and gamma values */
74 gamma_i=params->gamma[brand_i];
75 c_i=params->c[brand_i];
76 d_i=params->d[brand_i];
77 h_i=params->h[brand_i];
78 ci2=params->c2[brand_i];
79 di2=params->d2[brand_i];
80 ci2di2=params->c2d2[brand_i];
86 /* albe needs i,k depending c,d,h and gamma values */
87 gamma_i=params->gamma_m;
93 ci2di2=params->c2d2_m;
97 v3_sub(&dist_ik,&(ktom->r),&(ai->r));
98 if(bc_ik) check_per_bound(moldyn,&dist_ik);
99 d_ik2=v3_absolute_square(&dist_ik);
101 /* store data for second k loop */
102 exchange->dist_ik[kcount]=dist_ik;
103 exchange->d_ik2[kcount]=d_ik2;
105 /* return if not within cutoff */
115 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
118 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
119 d2_h_cos2=exchange->di2+(h_cos*h_cos);
120 frac=exchange->ci2/d2_h_cos2;
121 g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
122 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
125 h_cos=h_i+cos_theta; // + in albe formalism
126 d2_h_cos2=di2+(h_cos*h_cos);
128 g=gamma_i*(1.0+ci2di2-frac);
129 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
131 /* zeta sum += f_c_ik * g_ijk */
139 arg=M_PI*(d_ik-Rk)/s_r;
140 f_c_ik=0.5+0.5*cos(arg);
141 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
145 /* store even more data for second k loop */
146 exchange->g[kcount]=g;
147 exchange->dg[kcount]=dg;
148 exchange->d_ik[kcount]=d_ik;
149 exchange->cos_theta[kcount]=cos_theta;
150 exchange->f_c_ik[kcount]=f_c_ik;
151 exchange->df_c_ik[kcount]=df_c_ik;
153 /* increase k counter */
160 int albe_potential_force_calc(t_moldyn *moldyn) {
163 t_atom *itom,*jtom,*ktom;
167 int *neighbour_i[27];
173 t_list neighbour_i[27];
174 t_list neighbour_i2[27];
182 pthread_t kthread[27];
189 t_albe_mult_params *params;
190 t_albe_exchange *exchange;
204 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
205 double f_c_ik,df_c_ik;
208 double f_a,df_a,b,db,f_c,df_c;
217 double dijdik_inv,fcdg,dfcg;
218 t_3dvec dcosdrj,dcosdrk;
235 params=moldyn->pot_params;
236 exchange=&(params->exchange);
242 /* reset global virial */
243 memset(&(moldyn->gvir),0,sizeof(t_virial));
245 /* reset force, site energy and virial of every atom */
247 #pragma omp parallel for private(virial)
249 for(i=0;i<count;i++) {
252 v3_zero(&(itom[i].f));
255 virial=(&(itom[i].virial));
263 /* reset site energy */
268 /* get energy, force and virial of every atom */
270 /* first (and only) loop over atoms i */
271 for(i=0;i<count;i++) {
273 if(!(itom[i].attr&ATOM_ATTR_3BP))
276 link_cell_neighbour_index(moldyn,
277 (itom[i].r.x+moldyn->dim.x/2)/lc->x,
278 (itom[i].r.y+moldyn->dim.y/2)/lc->y,
279 (itom[i].r.z+moldyn->dim.z/2)/lc->z,
284 /* copy the neighbour lists */
288 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
294 /* loop over atoms j */
301 while(neighbour_i[j][p]!=-1) {
303 jtom=&(itom[neighbour_i[j][p]]);
311 p=lc->subcell->list[p];
313 this=&(neighbour_i[j]);
316 if(this->start==NULL)
321 jtom=this->current->data;
327 if(!(jtom->attr&ATOM_ATTR_3BP))
334 /* j1 func here ... */
335 /* albe 3 body potential function (first ij loop) */
341 * set ij depending values
344 if(brand_i==jtom->brand) {
345 S2=params->S2[brand_i];
352 v3_sub(&dist_ij,&(jtom->r),&(ai->r));
353 if(bc_ij) check_per_bound(moldyn,&dist_ij);
354 d_ij2=v3_absolute_square(&dist_ij);
356 /* if d_ij2 > S2 => no force & potential energy contribution */
363 /* reset k counter for first k loop */
366 /* first loop over atoms k */
373 kdata.itom=&(itom[i]);
374 ret=pthread_create(&(kthread[k]),NULL,k1_thread,&(kdata[k]));
376 perror("[albe fast] thread create");
384 while(neighbour_i[k][q]!=-1) {
386 ktom=&(itom[neighbour_i[k][q]]);
394 q=lc->subcell->list[q];
396 that=&(neighbour_i2[k]);
399 if(that->start==NULL)
403 ktom=that->current->data;
406 if(!(ktom->attr&ATOM_ATTR_3BP))
417 ret=pthread_create(&(k_thread[k]),NULL,k1_calc,k_data);
419 perror("[albe fast] create thread\n");
424 /* k1 func here ... */
425 /* albe 3 body potential function (first k loop) */
427 if(kcount>ALBE_MAXN) {
428 printf("FATAL: neighbours = %d\n",kcount);
429 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
433 if(brand_i==ktom->brand) {
434 Rk=params->R[brand_i];
435 Sk=params->S[brand_i];
436 Sk2=params->S2[brand_i];
437 /* albe needs i,k depending c,d,h and gamma values */
438 gamma_i=params->gamma[brand_i];
439 c_i=params->c[brand_i];
440 d_i=params->d[brand_i];
441 h_i=params->h[brand_i];
442 ci2=params->c2[brand_i];
443 di2=params->d2[brand_i];
444 ci2di2=params->c2d2[brand_i];
450 /* albe needs i,k depending c,d,h and gamma values */
451 gamma_i=params->gamma_m;
455 ci2=params->c2_mixed;
456 di2=params->d2_mixed;
457 ci2di2=params->c2d2_m;
461 v3_sub(&dist_ik,&(ktom->r),&(ai->r));
462 if(bc_ik) check_per_bound(moldyn,&dist_ik);
463 d_ik2=v3_absolute_square(&dist_ik);
465 /* store data for second k loop */
466 exchange->dist_ik[kcount]=dist_ik;
467 exchange->d_ik2[kcount]=d_ik2;
469 /* return if not within cutoff */
479 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
482 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
483 d2_h_cos2=exchange->di2+(h_cos*h_cos);
484 frac=exchange->ci2/d2_h_cos2;
485 g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
486 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
489 h_cos=h_i+cos_theta; // + in albe formalism
490 d2_h_cos2=di2+(h_cos*h_cos);
492 g=gamma_i*(1.0+ci2di2-frac);
493 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
495 /* zeta sum += f_c_ik * g_ijk */
503 arg=M_PI*(d_ik-Rk)/s_r;
504 f_c_ik=0.5+0.5*cos(arg);
505 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
509 /* store even more data for second k loop */
510 exchange->g[kcount]=g;
511 exchange->dg[kcount]=dg;
512 exchange->d_ik[kcount]=d_ik;
513 exchange->cos_theta[kcount]=cos_theta;
514 exchange->f_c_ik[kcount]=f_c_ik;
515 exchange->df_c_ik[kcount]=df_c_ik;
517 /* increase k counter */
525 } while(list_next_f(that)!=\
536 ret=pthread_join(kthread[k],NULL);
538 perror("[albe fast] join thread");
545 /* j2 func here ... */
548 if(brand_i==jtom->brand) {
549 S=params->S[brand_i];
550 R=params->R[brand_i];
551 B=params->B[brand_i];
552 A=params->A[brand_i];
553 r0=params->r0[brand_i];
554 mu=params->mu[brand_i];
555 lambda=params->lambda[brand_i];
564 lambda=params->lambda_m;
574 arg=M_PI*(d_ij-R)/s_r;
575 f_c=0.5+0.5*cos(arg);
576 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
580 f_a=-B*exp(-mu*(d_ij-r0));
584 f_r=A*exp(-lambda*(d_ij-r0));
585 df_r=lambda*f_r/d_ij;
593 b=1.0/sqrt(1.0+zeta_ij);
594 db=-0.5*b/(1.0+zeta_ij);
597 /* force contribution for atom i */
598 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
599 v3_scale(&force,&(dist_ij),scale);
600 v3_add(&(ai->f),&(ai->f),&force);
602 /* force contribution for atom j */
603 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
604 v3_add(&(jtom->f),&(jtom->f),&force);
607 virial_calc(ai,&force,&(dist_ij));
610 if(moldyn->time>DSTART&&moldyn->time<DEND) {
611 if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
612 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
613 printf(" adding %f %f %f\n",force.x,force.y,force.z);
614 if(ai==&(moldyn->atom[0]))
615 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
616 if(jtom==&(moldyn->atom[0]))
617 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
618 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
620 printf(" %f %f %f\n",zeta_ij,.0,.0);
625 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
626 pre_dzeta=0.5*f_a*f_c*db;
628 /* energy contribution */
629 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
630 moldyn->energy+=energy;
633 /* reset k counter for second k loop */
637 /* second loop over atoms k */
644 while(neighbour_i[k][q]!=-1) {
646 ktom=&(itom[neighbour_i[k][q]]);
654 q=lc->subcell->list[q];
656 that=&(neighbour_i2[k]);
659 if(that->start==NULL)
663 ktom=that->current->data;
666 if(!(ktom->attr&ATOM_ATTR_3BP))
676 /* k2 func here ... */
677 /* albe 3 body potential function (second k loop) */
680 printf("FATAL: neighbours!\n");
683 d_ik2=exchange->d_ik2[kcount];
685 if(brand_i==ktom->brand)
686 Sk2=params->S2[brand_i];
690 /* return if d_ik > S */
697 dist_ik=exchange->dist_ik[kcount];
698 d_ik=exchange->d_ik[kcount];
700 /* f_c_ik, df_c_ik */
701 f_c_ik=exchange->f_c_ik[kcount];
702 df_c_ik=exchange->df_c_ik[kcount];
704 /* g, dg, cos_theta */
705 g=exchange->g[kcount];
706 dg=exchange->dg[kcount];
707 cos_theta=exchange->cos_theta[kcount];
709 /* cos_theta derivatives wrt j,k */
710 dijdik_inv=1.0/(d_ij*d_ik);
711 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
712 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
713 v3_add(&dcosdrj,&dcosdrj,&tmp);
714 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
715 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
716 v3_add(&dcosdrk,&dcosdrk,&tmp);
718 /* f_c_ik * dg, df_c_ik * g */
722 /* derivative wrt j */
723 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
725 /* force contribution */
726 v3_add(&(jtom->f),&(jtom->f),&force);
729 if(moldyn->time>DSTART&&moldyn->time<DEND) {
730 if(jtom==&(moldyn->atom[DATOM])) {
731 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
732 printf(" adding %f %f %f\n",force.x,force.y,force.z);
733 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
734 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
735 printf(" d ij ik = %f %f\n",d_ij,d_ik);
741 virial_calc(ai,&force,&dist_ij);
743 /* force contribution to atom i */
744 v3_scale(&force,&force,-1.0);
745 v3_add(&(ai->f),&(ai->f),&force);
747 /* derivative wrt k */
748 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
749 v3_scale(&tmp,&dcosdrk,fcdg);
750 v3_add(&force,&force,&tmp);
751 v3_scale(&force,&force,pre_dzeta);
753 /* force contribution */
754 v3_add(&(ktom->f),&(ktom->f),&force);
757 if(moldyn->time>DSTART&&moldyn->time<DEND) {
758 if(ktom==&(moldyn->atom[DATOM])) {
759 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
760 printf(" adding %f %f %f\n",force.x,force.y,force.z);
761 printf(" total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
762 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
763 printf(" d ij ik = %f %f\n",d_ij,d_ik);
769 virial_calc(ai,&force,&dist_ik);
771 /* force contribution to atom i */
772 v3_scale(&force,&force,-1.0);
773 v3_add(&(ai->f),&(ai->f),&force);
775 /* increase k counter */
785 } while(list_next_f(that)!=\
796 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
811 //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
812 if(moldyn->time>DSTART&&moldyn->time<DEND) {
814 printf(" x: %0.40f\n",moldyn->atom[DATOM].f.x);
815 printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y);
816 printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z);
820 /* some postprocessing */
822 #pragma omp parallel for
824 for(i=0;i<count;i++) {
825 /* calculate global virial */
826 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
827 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
828 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
829 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
830 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
831 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
833 /* check forces regarding the given timestep */
834 if(v3_norm(&(itom[i].f))>\
835 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
836 printf("[moldyn] WARNING: pfc (high force: atom %d)\n",