X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=potentials%2Falbe.c;h=fea7ba818fdeb9562e895166ef181c078a740929;hb=b5b47daaa3718c4dec2056fe5147668023575b8e;hp=c2acf2f3cfe9902a8f781fe7920dad060fa637cc;hpb=a70de3dccbf0a01c39c2643818ec86c0b465caba;p=physik%2Fposic.git diff --git a/potentials/albe.c b/potentials/albe.c index c2acf2f..fea7ba8 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -43,34 +43,6 @@ int albe_mult_complete_params(t_albe_mult_params *p) { return 0; } -/* albe 1 body part */ -int albe_mult_1bp(t_moldyn *moldyn,t_atom *ai) { - - int brand; - t_albe_mult_params *params; - t_albe_exchange *exchange; - - brand=ai->brand; - params=moldyn->pot_params; - exchange=&(params->exchange); - - /* - * simple: point constant parameters only depending on atom i to - * their right values - */ - - exchange->gamma_i=&(params->gamma[brand]); - exchange->c_i=&(params->c[brand]); - exchange->d_i=&(params->d[brand]); - exchange->h_i=&(params->h[brand]); - - exchange->ci2=params->c[brand]*params->c[brand]; - exchange->di2=params->d[brand]*params->d[brand]; - exchange->ci2di2=exchange->ci2/exchange->di2; - - return 0; -} - /* albe 3 body potential function (first ij loop) */ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { @@ -92,11 +64,12 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { */ brand=ai->brand; - - if(brand==aj->brand) + if(brand==aj->brand) { S2=params->S2[brand]; - else + } + else { S2=params->S2mixed; + } /* dist_ij, d_ij2 */ v3_sub(&dist_ij,&(aj->r),&(ai->r)); @@ -108,6 +81,8 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { moldyn->run3bp=0; return 0; } +//if(ai==&(moldyn->atom[5832])) +// printf(" %d: %d (%f)\n",ai->tag,aj->tag,moldyn->time); /* d_ij */ d_ij=sqrt(d_ij2); @@ -125,7 +100,7 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* albe 3 body potential function (first k loop) */ int albe_mult_3bp_k1(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; @@ -152,12 +127,25 @@ int albe_mult_3bp_k1(t_moldyn *moldyn, R=params->R[brand]; S=params->S[brand]; S2=params->S2[brand]; + /* albe needs i,k depending c,d,h and gamma values */ + exchange->gamma_i=&(params->gamma[brand]); + exchange->c_i=&(params->c[brand]); + exchange->d_i=&(params->d[brand]); + exchange->h_i=&(params->h[brand]); } else { R=params->Rmixed; S=params->Smixed; S2=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); } + 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,&(ak->r),&(ai->r)); @@ -185,11 +173,11 @@ int albe_mult_3bp_k1(t_moldyn *moldyn, cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); /* g_ijk */ - h_cos=*(exchange->h_i)-cos_theta; + 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; + dg=2.0*frac**(exchange->gamma_i)*h_cos/d2_h_cos2; // + in albe f.. /* zeta sum += f_c_ik * g_ijk */ if(d_ik<=R) { @@ -232,6 +220,7 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { double d_ij,r0; unsigned char brand; double S,R,s_r,arg; + double energy; params=moldyn->pot_params; exchange=&(params->exchange); @@ -289,13 +278,14 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { } /* force contribution */ - scale=-0.5*(f_c*(df_r+b*df_a)+df_c*(f_r+b*df_a)); + scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); v3_scale(&force,&(exchange->dist_ij),scale); v3_add(&(ai->f),&(ai->f),&force); v3_sub(&(aj->f),&(aj->f),&force); // dri rij = - drj rij #ifdef DEBUG - if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) { +if(moldyn->time>DSTART&&moldyn->timeatom[5832]))|(aj==&(moldyn->atom[5832]))) { printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag); printf("adding %f %f %f\n",force.x,force.y,force.z); if(ai==&(moldyn->atom[0])) @@ -304,18 +294,21 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r), f_c,b,f_a,f_r); + printf(" %f %f %f\n",exchange->zeta_ij,.0,.0); } +} #endif /* virial */ - if(ajdist_ij)); + virial_calc(ai,&force,&(exchange->dist_ij)); - /* dzeta prefactor = - 0.5 f_c f_a db */ - exchange->pre_dzeta=-0.5*f_a*f_c*db; + /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */ + exchange->pre_dzeta=0.5*f_a*f_c*db; /* energy contribution */ - moldyn->energy+=0.5*f_c*(f_r+b*f_a); + energy=0.5*f_c*(f_r-b*f_a); + moldyn->energy+=energy; + ai->e+=energy; /* reset k counter for second k loop */ exchange->kcount=0; @@ -325,7 +318,7 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* albe 3 body potential function (second k loop) */ int albe_mult_3bp_k2(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; @@ -386,13 +379,13 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, /* cos_theta derivatives wrt i,j,k */ dijdik_inv=1.0/(d_ij*d_ik); - v3_scale(&dcosdrj,&dist_ik,dijdik_inv); + v3_scale(&dcosdrj,&dist_ik,dijdik_inv); // j v3_scale(&tmp,&dist_ij,-cos_theta/d_ij2); v3_add(&dcosdrj,&dcosdrj,&tmp); - v3_scale(&dcosdrk,&dist_ij,dijdik_inv); + v3_scale(&dcosdrk,&dist_ij,dijdik_inv); // k v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2); v3_add(&dcosdrk,&dcosdrk,&tmp); - v3_add(&dcosdri,&dcosdrj,&dcosdrk); + v3_add(&dcosdri,&dcosdrj,&dcosdrk); // i v3_scale(&dcosdri,&dcosdri,-1.0); /* f_c_ik * dg, df_c_ik * g */ @@ -409,11 +402,13 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, v3_add(&(ai->f),&(ai->f),&force); #ifdef DEBUG - if(ai==&(moldyn->atom[0])) { +if(moldyn->time>DSTART&&moldyn->timeatom[5832])) { printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag); printf("adding %f %f %f\n",force.x,force.y,force.z); printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); } +} #endif /* virial */ @@ -426,17 +421,18 @@ int albe_mult_3bp_k2(t_moldyn *moldyn, v3_add(&(aj->f),&(aj->f),&force); #ifdef DEBUG - if(aj==&(moldyn->atom[0])) { +if(moldyn->time>DSTART&&moldyn->timeatom[5832])) { printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag); printf("adding %f %f %f\n",force.x,force.y,force.z); printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z); } +} #endif /* virial */ - //v3_scale(&force,&force,-1.0); - if(ajf),&(ak->f),&force); #ifdef DEBUG - if(ak==&(moldyn->atom[0])) { +if(moldyn->time>DSTART&&moldyn->timeatom[5832])) { printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag); printf("adding %f %f %f\n",force.x,force.y,force.z); printf("total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z); } +} #endif /* virial */ - //v3_scale(&force,&force,-1.0); - if(ajkcount++;