4 * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
12 #include <sys/types.h>
28 #include "../moldyn.h"
29 #include "../math/math.h"
36 #define albe_v_calc(a,f,d) a->virial.xx+=f->x*d->x; \
37 a->virial.yy+=f->y*d->y; \
38 a->virial.zz+=f->z*d->z; \
39 a->virial.xy+=f->x*d->y; \
40 a->virial.xz+=f->x*d->z; \
41 a->virial.yz+=f->y*d->z
45 void *k1_calc(void *ptr) {
47 /* albe 3 body potential function (first k loop) */
49 t_albe_mult_params *params;
50 t_albe_exchange *exchange;
51 unsigned char brand_i;
52 double Rk,Sk,Sk2,gamma_i,c_i,d_i,h_i,ci2,di2,ci2di2;
53 t_atom *ai,*jtom,*ktom;
56 if(kcount>ALBE_MAXN) {
57 printf("FATAL: neighbours = %d\n",kcount);
58 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
62 if(brand_i==ktom->brand) {
63 Rk=params->R[brand_i];
64 Sk=params->S[brand_i];
65 Sk2=params->S2[brand_i];
66 /* albe needs i,k depending c,d,h and gamma values */
67 gamma_i=params->gamma[brand_i];
68 c_i=params->c[brand_i];
69 d_i=params->d[brand_i];
70 h_i=params->h[brand_i];
71 ci2=params->c2[brand_i];
72 di2=params->d2[brand_i];
73 ci2di2=params->c2d2[brand_i];
79 /* albe needs i,k depending c,d,h and gamma values */
80 gamma_i=params->gamma_m;
86 ci2di2=params->c2d2_m;
90 v3_sub(&dist_ik,&(ktom->r),&(ai->r));
91 if(bc_ik) check_per_bound(moldyn,&dist_ik);
92 d_ik2=v3_absolute_square(&dist_ik);
94 /* store data for second k loop */
95 exchange->dist_ik[kcount]=dist_ik;
96 exchange->d_ik2[kcount]=d_ik2;
98 /* return if not within cutoff */
108 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
111 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
112 d2_h_cos2=exchange->di2+(h_cos*h_cos);
113 frac=exchange->ci2/d2_h_cos2;
114 g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
115 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
118 h_cos=h_i+cos_theta; // + in albe formalism
119 d2_h_cos2=di2+(h_cos*h_cos);
121 g=gamma_i*(1.0+ci2di2-frac);
122 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
124 /* zeta sum += f_c_ik * g_ijk */
132 arg=M_PI*(d_ik-Rk)/s_r;
133 f_c_ik=0.5+0.5*cos(arg);
134 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
138 /* store even more data for second k loop */
139 exchange->g[kcount]=g;
140 exchange->dg[kcount]=dg;
141 exchange->d_ik[kcount]=d_ik;
142 exchange->cos_theta[kcount]=cos_theta;
143 exchange->f_c_ik[kcount]=f_c_ik;
144 exchange->df_c_ik[kcount]=df_c_ik;
146 /* increase k counter */
153 int albe_potential_force_calc(t_moldyn *moldyn) {
156 t_atom *itom,*jtom,*ktom;
160 int *neighbour_i[27];
166 t_list neighbour_i[27];
167 t_list neighbour_i2[27];
177 t_albe_mult_params *params;
178 t_albe_exchange *exchange;
192 double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg;
193 double f_c_ik,df_c_ik;
196 double f_a,df_a,b,db,f_c,df_c;
205 double dijdik_inv,fcdg,dfcg;
206 t_3dvec dcosdrj,dcosdrk;
223 params=moldyn->pot_params;
224 exchange=&(params->exchange);
230 /* reset global virial */
231 memset(&(moldyn->gvir),0,sizeof(t_virial));
233 /* reset force, site energy and virial of every atom */
235 #pragma omp parallel for private(virial)
237 for(i=0;i<count;i++) {
240 v3_zero(&(itom[i].f));
243 virial=(&(itom[i].virial));
251 /* reset site energy */
256 /* get energy, force and virial of every atom */
258 /* first (and only) loop over atoms i */
259 for(i=0;i<count;i++) {
261 if(!(itom[i].attr&ATOM_ATTR_3BP))
264 link_cell_neighbour_index(moldyn,
265 (itom[i].r.x+moldyn->dim.x/2)/lc->x,
266 (itom[i].r.y+moldyn->dim.y/2)/lc->y,
267 (itom[i].r.z+moldyn->dim.z/2)/lc->z,
272 /* copy the neighbour lists */
276 memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
282 /* loop over atoms j */
289 while(neighbour_i[j][p]!=-1) {
291 jtom=&(itom[neighbour_i[j][p]]);
299 p=lc->subcell->list[p];
301 this=&(neighbour_i[j]);
304 if(this->start==NULL)
309 jtom=this->current->data;
315 if(!(jtom->attr&ATOM_ATTR_3BP))
322 /* j1 func here ... */
323 /* albe 3 body potential function (first ij loop) */
329 * set ij depending values
332 if(brand_i==jtom->brand) {
333 S2=params->S2[brand_i];
340 v3_sub(&dist_ij,&(jtom->r),&(ai->r));
341 if(bc_ij) check_per_bound(moldyn,&dist_ij);
342 d_ij2=v3_absolute_square(&dist_ij);
344 /* if d_ij2 > S2 => no force & potential energy contribution */
351 /* reset k counter for first k loop */
354 /* first loop over atoms k */
361 while(neighbour_i[k][q]!=-1) {
363 ktom=&(itom[neighbour_i[k][q]]);
371 q=lc->subcell->list[q];
373 that=&(neighbour_i2[k]);
376 if(that->start==NULL)
380 ktom=that->current->data;
383 if(!(ktom->attr&ATOM_ATTR_3BP))
394 ret=pthread_create(&(k_thread[k]),NULL,k1_calc,k_data);
396 perror("[albe fast] create thread\n");
401 /* k1 func here ... */
402 /* albe 3 body potential function (first k loop) */
404 if(kcount>ALBE_MAXN) {
405 printf("FATAL: neighbours = %d\n",kcount);
406 printf(" -> %d %d %d\n",ai->tag,jtom->tag,ktom->tag);
410 if(brand_i==ktom->brand) {
411 Rk=params->R[brand_i];
412 Sk=params->S[brand_i];
413 Sk2=params->S2[brand_i];
414 /* albe needs i,k depending c,d,h and gamma values */
415 gamma_i=params->gamma[brand_i];
416 c_i=params->c[brand_i];
417 d_i=params->d[brand_i];
418 h_i=params->h[brand_i];
419 ci2=params->c2[brand_i];
420 di2=params->d2[brand_i];
421 ci2di2=params->c2d2[brand_i];
427 /* albe needs i,k depending c,d,h and gamma values */
428 gamma_i=params->gamma_m;
432 ci2=params->c2_mixed;
433 di2=params->d2_mixed;
434 ci2di2=params->c2d2_m;
438 v3_sub(&dist_ik,&(ktom->r),&(ai->r));
439 if(bc_ik) check_per_bound(moldyn,&dist_ik);
440 d_ik2=v3_absolute_square(&dist_ik);
442 /* store data for second k loop */
443 exchange->dist_ik[kcount]=dist_ik;
444 exchange->d_ik2[kcount]=d_ik2;
446 /* return if not within cutoff */
456 cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
459 h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism
460 d2_h_cos2=exchange->di2+(h_cos*h_cos);
461 frac=exchange->ci2/d2_h_cos2;
462 g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac);
463 dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f..
466 h_cos=h_i+cos_theta; // + in albe formalism
467 d2_h_cos2=di2+(h_cos*h_cos);
469 g=gamma_i*(1.0+ci2di2-frac);
470 dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f..
472 /* zeta sum += f_c_ik * g_ijk */
480 arg=M_PI*(d_ik-Rk)/s_r;
481 f_c_ik=0.5+0.5*cos(arg);
482 df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik));
486 /* store even more data for second k loop */
487 exchange->g[kcount]=g;
488 exchange->dg[kcount]=dg;
489 exchange->d_ik[kcount]=d_ik;
490 exchange->cos_theta[kcount]=cos_theta;
491 exchange->f_c_ik[kcount]=f_c_ik;
492 exchange->df_c_ik[kcount]=df_c_ik;
494 /* increase k counter */
504 } while(list_next_f(that)!=\
511 /* j2 func here ... */
514 if(brand_i==jtom->brand) {
515 S=params->S[brand_i];
516 R=params->R[brand_i];
517 B=params->B[brand_i];
518 A=params->A[brand_i];
519 r0=params->r0[brand_i];
520 mu=params->mu[brand_i];
521 lambda=params->lambda[brand_i];
530 lambda=params->lambda_m;
540 arg=M_PI*(d_ij-R)/s_r;
541 f_c=0.5+0.5*cos(arg);
542 df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij));
546 f_a=-B*exp(-mu*(d_ij-r0));
550 f_r=A*exp(-lambda*(d_ij-r0));
551 df_r=lambda*f_r/d_ij;
559 b=1.0/sqrt(1.0+zeta_ij);
560 db=-0.5*b/(1.0+zeta_ij);
563 /* force contribution for atom i */
564 scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism
565 v3_scale(&force,&(dist_ij),scale);
566 v3_add(&(ai->f),&(ai->f),&force);
568 /* force contribution for atom j */
569 v3_scale(&force,&force,-1.0); // dri rij = - drj rij
570 v3_add(&(jtom->f),&(jtom->f),&force);
573 virial_calc(ai,&force,&(dist_ij));
576 if(moldyn->time>DSTART&&moldyn->time<DEND) {
577 if((ai==&(moldyn->atom[DATOM]))|(jtom==&(moldyn->atom[DATOM]))) {
578 printf("force 3bp (j2): [%d %d sum]\n",ai->tag,jtom->tag);
579 printf(" adding %f %f %f\n",force.x,force.y,force.z);
580 if(ai==&(moldyn->atom[0]))
581 printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
582 if(jtom==&(moldyn->atom[0]))
583 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
584 printf(" energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
586 printf(" %f %f %f\n",zeta_ij,.0,.0);
591 /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
592 pre_dzeta=0.5*f_a*f_c*db;
594 /* energy contribution */
595 energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism
596 moldyn->energy+=energy;
599 /* reset k counter for second k loop */
603 /* second loop over atoms k */
610 while(neighbour_i[k][q]!=-1) {
612 ktom=&(itom[neighbour_i[k][q]]);
620 q=lc->subcell->list[q];
622 that=&(neighbour_i2[k]);
625 if(that->start==NULL)
629 ktom=that->current->data;
632 if(!(ktom->attr&ATOM_ATTR_3BP))
642 /* k2 func here ... */
643 /* albe 3 body potential function (second k loop) */
646 printf("FATAL: neighbours!\n");
649 d_ik2=exchange->d_ik2[kcount];
651 if(brand_i==ktom->brand)
652 Sk2=params->S2[brand_i];
656 /* return if d_ik > S */
663 dist_ik=exchange->dist_ik[kcount];
664 d_ik=exchange->d_ik[kcount];
666 /* f_c_ik, df_c_ik */
667 f_c_ik=exchange->f_c_ik[kcount];
668 df_c_ik=exchange->df_c_ik[kcount];
670 /* g, dg, cos_theta */
671 g=exchange->g[kcount];
672 dg=exchange->dg[kcount];
673 cos_theta=exchange->cos_theta[kcount];
675 /* cos_theta derivatives wrt j,k */
676 dijdik_inv=1.0/(d_ij*d_ik);
677 v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j
678 v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2);
679 v3_add(&dcosdrj,&dcosdrj,&tmp);
680 v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k
681 v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
682 v3_add(&dcosdrk,&dcosdrk,&tmp);
684 /* f_c_ik * dg, df_c_ik * g */
688 /* derivative wrt j */
689 v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
691 /* force contribution */
692 v3_add(&(jtom->f),&(jtom->f),&force);
695 if(moldyn->time>DSTART&&moldyn->time<DEND) {
696 if(jtom==&(moldyn->atom[DATOM])) {
697 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
698 printf(" adding %f %f %f\n",force.x,force.y,force.z);
699 printf(" total j: %f %f %f\n",jtom->f.x,jtom->f.y,jtom->f.z);
700 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
701 printf(" d ij ik = %f %f\n",d_ij,d_ik);
707 virial_calc(ai,&force,&dist_ij);
709 /* force contribution to atom i */
710 v3_scale(&force,&force,-1.0);
711 v3_add(&(ai->f),&(ai->f),&force);
713 /* derivative wrt k */
714 v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
715 v3_scale(&tmp,&dcosdrk,fcdg);
716 v3_add(&force,&force,&tmp);
717 v3_scale(&force,&force,pre_dzeta);
719 /* force contribution */
720 v3_add(&(ktom->f),&(ktom->f),&force);
723 if(moldyn->time>DSTART&&moldyn->time<DEND) {
724 if(ktom==&(moldyn->atom[DATOM])) {
725 printf("force 3bp (k2): [%d %d %d]\n",ai->tag,jtom->tag,ktom->tag);
726 printf(" adding %f %f %f\n",force.x,force.y,force.z);
727 printf(" total k: %f %f %f\n",ktom->f.x,ktom->f.y,ktom->f.z);
728 printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
729 printf(" d ij ik = %f %f\n",d_ij,d_ik);
735 virial_calc(ai,&force,&dist_ik);
737 /* force contribution to atom i */
738 v3_scale(&force,&force,-1.0);
739 v3_add(&(ai->f),&(ai->f),&force);
741 /* increase k counter */
751 } while(list_next_f(that)!=\
762 } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
777 //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
778 if(moldyn->time>DSTART&&moldyn->time<DEND) {
780 printf(" x: %0.40f\n",moldyn->atom[DATOM].f.x);
781 printf(" y: %0.40f\n",moldyn->atom[DATOM].f.y);
782 printf(" z: %0.40f\n",moldyn->atom[DATOM].f.z);
786 /* some postprocessing */
788 #pragma omp parallel for
790 for(i=0;i<count;i++) {
791 /* calculate global virial */
792 moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
793 moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
794 moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
795 moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
796 moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
797 moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
799 /* check forces regarding the given timestep */
800 if(v3_norm(&(itom[i].f))>\
801 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
802 printf("[moldyn] WARNING: pfc (high force: atom %d)\n",