From: hackbard Date: Mon, 1 Sep 2008 12:28:15 +0000 (+0200) Subject: supply c2, d2, c2/d2 X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=a198ee2eb93c0f14bdaea14722f4d98e5890935f;p=physik%2Fposic.git supply c2, d2, c2/d2 --- diff --git a/potentials/albe.c b/potentials/albe.c index 02c254f..08ce4fa 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -96,6 +96,15 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) { p->S2[0]=p->S[0]*p->S[0]; p->S2[1]=p->S[1]*p->S[1]; p->S2mixed=p->Smixed*p->Smixed; + p->c2[0]=p->c[0]*p->c[0]; + p->c2[1]=p->c[1]*p->c[1]; + p->c2_mixed=p->c_mixed*p->c_mixed; + p->d2[0]=p->d[0]*p->d[0]; + p->d2[1]=p->d[1]*p->d[1]; + p->d2_mixed=p->d_mixed*p->d_mixed; + p->c2d2[0]=p->c2[0]/p->d2[0]; + p->c2d2[1]=p->c2[1]/p->d2[1]; + p->c2d2_m=p->c2_mixed/p->d2_mixed; printf("[albe] mult parameter info:\n"); printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed); diff --git a/potentials/albe.h b/potentials/albe.h index 15386e7..87f45a0 100644 --- a/potentials/albe.h +++ b/potentials/albe.h @@ -66,10 +66,16 @@ typedef struct s_albe_mult_params { double gamma_m; double c[2]; double c_mixed; + double c2[2]; + double c2_mixed; double d[2]; double d_mixed; + double d2[2]; + double d2_mixed; double h[2]; double h_mixed; + double c2d2[2]; + double c2d2_m; t_albe_exchange exchange; /* exchange between 2bp and 3bp calc */ } t_albe_mult_params; diff --git a/potentials/albe_fast.c b/potentials/albe_fast.c index 3b5f681..6331064 100644 --- a/potentials/albe_fast.c +++ b/potentials/albe_fast.c @@ -90,6 +90,14 @@ int albe_potential_force_calc(t_moldyn *moldyn) { t_3dvec dcosdrj,dcosdrk; t_3dvec tmp; + // even more ... + double gamma_i; + double c_i; + double d_i; + double h_i; + double ci2; + double di2; + double ci2di2; count=moldyn->count; itom=moldyn->atom; @@ -267,24 +275,27 @@ int albe_potential_force_calc(t_moldyn *moldyn) { Sk=params->S[brand_i]; Sk2=params->S2[brand_i]; /* albe needs i,k depending c,d,h and gamma values */ - exchange->gamma_i=&(params->gamma[brand_i]); - exchange->c_i=&(params->c[brand_i]); - exchange->d_i=&(params->d[brand_i]); - exchange->h_i=&(params->h[brand_i]); + gamma_i=params->gamma[brand_i]; + c_i=params->c[brand_i]; + d_i=params->d[brand_i]; + h_i=params->h[brand_i]; + ci2=params->c2[brand_i]; + di2=params->d2[brand_i]; + ci2di2=params->c2d2[brand_i]; } else { Rk=params->Rmixed; Sk=params->Smixed; Sk2=params->S2mixed; /* albe needs i,k depending c,d,h and gamma values */ - exchange->gamma_i=&(params->gamma_m); - exchange->c_i=&(params->c_mixed); - exchange->d_i=&(params->d_mixed); - exchange->h_i=&(params->h_mixed); + gamma_i=params->gamma_m; + c_i=params->c_mixed; + d_i=params->d_mixed; + h_i=params->h_mixed; + ci2=params->c2_mixed; + di2=params->d2_mixed; + ci2di2=params->c2d2_m; } - exchange->ci2=*(exchange->c_i)**(exchange->c_i); - exchange->di2=*(exchange->d_i)**(exchange->d_i); - exchange->ci2di2=exchange->ci2/exchange->di2; /* dist_ik, d_ik2 */ v3_sub(&dist_ik,&(ktom->r),&(ai->r)); @@ -307,12 +318,19 @@ int albe_potential_force_calc(t_moldyn *moldyn) { /* cos theta */ cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); - /* g_ijk */ + /* g_ijk h_cos=*(exchange->h_i)+cos_theta; // + in albe formalism d2_h_cos2=exchange->di2+(h_cos*h_cos); frac=exchange->ci2/d2_h_cos2; g=*(exchange->gamma_i)*(1.0+exchange->ci2di2-frac); dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f.. + */ + + h_cos=h_i+cos_theta; // + in albe formalism + d2_h_cos2=di2+(h_cos*h_cos); + frac=ci2/d2_h_cos2; + g=gamma_i*(1.0+ci2di2-frac); + dg=2.0*frac*gamma_i*h_cos/d2_h_cos2; // + in albe f.. /* zeta sum += f_c_ik * g_ijk */ if(d_ik<=Rk) {