From 959801961cfdd51e080e6500f51647b3397d336c Mon Sep 17 00:00:00 2001 From: hackbard Date: Thu, 10 Jul 2008 22:13:33 +0200 Subject: [PATCH] first init of new albe implementation (still warnings!) --- potentials/albe.c | 532 +++++++++++++++++++++------------------------- potentials/albe.h | 54 ++--- 2 files changed, 270 insertions(+), 316 deletions(-) diff --git a/potentials/albe.c b/potentials/albe.c index c01ab15..9dfba13 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -113,26 +113,47 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) { 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) { +/* first i loop */ +int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) { + int i; t_albe_mult_params *params; t_albe_exchange *exchange; - unsigned char brand; - double S2; - t_3dvec dist_ij; - double d_ij2,d_ij; - + params=moldyn->pot_params; exchange=&(params->exchange); - /* reset zeta sum */ - exchange->zeta_ij=0.0; + /* zero exchange values */ + memset(exchange->dist,0,ALBE_MAXN*sizeof(t_3dvec)); + memset(exchange->d2,0,ALBE_MAXN*sizeof(double)); + memset(exchange->d,0,ALBE_MAXN*sizeof(double)); + memset(exchange->zeta,0,ALBE_MAXN*sizeof(double)); + for(i=0;idzeta[i],0,ALBE_MAXN*sizeof(double)); + memset(exchange->skip,0,ALBE_MAXN*sizeof(u8)); + exchange->jcnt=0; + exchange->j2cnt=0; + + return 0; +} + +/* first j loop within first i loop */ +int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_albe_mult_params *params; + t_albe_exchange *exchange; + double S2,S,R,d2,d,s_r,arg; + t_3dvec dist; + int j; + u8 brand; + + params=moldyn->pot_params; + exchange=&(params->exchange); - /* - * set ij depending values - */ + /* get j counter */ + j=exchange->jcnt; + /* set ij depending values */ brand=ai->brand; if(brand==aj->brand) { S2=params->S2[brand]; @@ -142,161 +163,207 @@ int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { } /* dist_ij, d_ij2 */ - v3_sub(&dist_ij,&(aj->r),&(ai->r)); - if(bc) check_per_bound(moldyn,&dist_ij); - d_ij2=v3_absolute_square(&dist_ij); + v3_sub(&dist,&(aj->r),&(ai->r)); + exchange->dist[j]=dist; + if(bc) check_per_bound(moldyn,&dist); + d2=v3_absolute_square(&dist); + exchange->d2[j]=d2; /* if d_ij2 > S2 => no force & potential energy contribution */ - if(d_ij2>S2) { + if(d2>S2) { moldyn->run3bp=0; + exchange->skip[j]=1; + exchange->jcnt+=1; return 0; } + /* more ij depending values */ + if(brand==aj->brand) { + R=params->R[brand]; + S=params->S[brand]; + /* set albe needs i,(j/k) depending c,d,h and gamma values */ + exchange->gamma_[j]=&(params->gamma[brand]); + exchange->c_[j]=&(params->c[brand]); + exchange->d_[j]=&(params->d[brand]); + exchange->h_[j]=&(params->h[brand]); + } + else { + R=params->Rmixed; + S=params->Smixed; + /* albe needs i,(j/k) depending c,d,h and gamma values */ + exchange->gamma_[j]=&(params->gamma_m); + exchange->c_[j]=&(params->c_mixed); + exchange->d_[j]=&(params->d_mixed); + exchange->h_[j]=&(params->h_mixed); + } + exchange->c2_[j]=*exchange->c_[j]**exchange->c_[j]; + exchange->d2_[j]=*exchange->d_[j]**exchange->d_[j]; + exchange->c2d2_[j]=exchange->c2_[j]/exchange->d2_[j]; + /* d_ij */ - d_ij=sqrt(d_ij2); + d=sqrt(exchange->d2[j]); + exchange->d[j]=d; + + /* f_c, df_c */ + if(df_c[j]=1.0; + exchange->df_c[j]=0.0; + } + else { + s_r=S-R; + arg=M_PI*(d-R)/s_r; + exchange->f_c[j]=0.5+0.5*cos(arg); + exchange->df_c[j]=0.5*sin(arg)*(M_PI/(s_r*d)); + } - /* store values */ - exchange->dist_ij=dist_ij; - exchange->d_ij2=d_ij2; - exchange->d_ij=d_ij; + /* reset k counter */ + exchange->kcnt=0; - /* reset k counter for first k loop */ - exchange->kcount=0; - return 0; } -/* 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) { +/* first k loop within first j loop within first i loop */ +int albe_mult_i0_j0_k0(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; - unsigned char brand; - double R,S,S2; - t_3dvec dist_ij,dist_ik; - double d_ik2,d_ik,d_ij; - double cos_theta,h_cos,d2_h_cos2,frac,g,dg,s_r,arg; - double f_c_ik,df_c_ik; - int kcount; + int j,k; + t_3dvec distj,distk; + double dj,dk,djdk_inv,cos_theta; + double gj,dgj,h_cos_j,d2_h_cos2_j,frac_j; + double gk,dgk,h_cos_k,d2_h_cos2_k,frac_k; + t_3dvec dcosdrj,dcosdrk,tmp,**dzeta; params=moldyn->pot_params; exchange=&(params->exchange); - kcount=exchange->kcount; - if(kcount>ALBE_MAXN) { - printf("FATAL: neighbours = %d\n",kcount); - printf(" -> %d %d %d\n",ai->tag,aj->tag,ak->tag); - } - - /* ik constants */ - brand=ai->brand; - if(brand==ak->brand) { - 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]); + /* kjcnt; + k=exchange->kcnt; + if(k>=ALBE_MAXN) { + printf("FATAL: too many neighbours! (%d)\n",k); + printf(" atom i:%d | j:%d | k:%d\n",ai->tag,aj->tag,ak->tag); } - 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)); - if(bc) check_per_bound(moldyn,&dist_ik); - d_ik2=v3_absolute_square(&dist_ik); - - /* store data for second k loop */ - exchange->dist_ik[kcount]=dist_ik; - exchange->d_ik2[kcount]=d_ik2; - - /* return if not within cutoff */ - if(d_ik2>S2) { - exchange->kcount++; + if((k>=j)|(exchange->skip[k])) { + exchange->kcnt+=1; return 0; } - - /* d_ik */ - d_ik=sqrt(d_ik2); - - /* dist_ij, d_ij */ - dist_ij=exchange->dist_ij; - d_ij=exchange->d_ij; + + /* distances */ + distj=exchange->dist[j]; + distk=exchange->dist[k]; + dj=exchange->d[j]; + dk=exchange->d[k]; + djdk_inv=1.0/(dj*dk); /* cos theta */ - cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); - - /* 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.. - - /* zeta sum += f_c_ik * g_ijk */ - if(d_ik<=R) { - exchange->zeta_ij+=g; - f_c_ik=1.0; - df_c_ik=0.0; + cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv; + + /* g(cos(theta)) - in albe: ik-depending values! */ + h_cos_j=*(exchange->h_[j])+cos_theta; // + in albe formalism + d2_h_cos2_j=exchange->d2_[j]+(h_cos_j*h_cos_j); + frac_j=exchange->c2_[j]/d2_h_cos2_j; + gj=1.0+exchange->c2d2_[j]-frac_j; + gj*=*(exchange->gamma_[j]); + dgj=*(exchange->gamma_[j])*2.0*frac_j*h_cos_j/d2_h_cos2_j; // + in albe + if(ak->brand==aj->brand) { + gk=gj; + dgk=dgj; } else { - s_r=S-R; - arg=M_PI*(d_ik-R)/s_r; - f_c_ik=0.5+0.5*cos(arg); - df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik)); - exchange->zeta_ij+=f_c_ik*g; + h_cos_k=*(exchange->h_[k])+cos_theta; + d2_h_cos2_k=exchange->d2_[k]+(h_cos_k*h_cos_k); + frac_k=exchange->c2_[k]/d2_h_cos2_k; + gk=1.0+exchange->c2d2_[k]-frac_k; + gk*=*(exchange->gamma_[k]); + dgk=*(exchange->gamma_[k])*2.0*frac_k*h_cos_k/d2_h_cos2_k; } - /* store even more data for second k loop */ - exchange->g[kcount]=g; - exchange->dg[kcount]=dg; - exchange->d_ik[kcount]=d_ik; - exchange->cos_theta[kcount]=cos_theta; - exchange->f_c_ik[kcount]=f_c_ik; - exchange->df_c_ik[kcount]=df_c_ik; + /* zeta */ + exchange->zeta[j]+=(exchange->f_c[k]*gk); + exchange->zeta[k]+=(exchange->f_c[j]*gj); + + /* cos theta derivatives */ + v3_scale(&dcosdrj,&distk,djdk_inv); // j + v3_scale(&tmp,&distj,-cos_theta/exchange->d2[j]); + v3_add(&dcosdrj,&dcosdrj,&tmp); + v3_scale(&dcosdrk,&distj,djdk_inv); // k + v3_scale(&tmp,&distk,-cos_theta/exchange->d2[k]); + v3_add(&dcosdrk,&dcosdrk,&tmp); + + /* zeta derivatives */ + dzeta=exchange->dzeta; + v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk); + v3_add(&dzeta[j][j],&dzeta[j][j],&tmp); // j j + v3_scale(&tmp,&dcosdrk,exchange->f_c[j]*dgj); + v3_add(&dzeta[k][k],&dzeta[k][k],&tmp); // k k + v3_scale(&tmp,&distk,exchange->df_c[k]*gk/dk); + v3_add(&dzeta[j][k],&dzeta[j][k],&tmp); + v3_scale(&tmp,&dcosdrk,exchange->f_c[k]*dgk); + v3_add(&dzeta[j][k],&dzeta[j][k],&tmp); // j k + v3_scale(&tmp,&distj,exchange->df_c[j]*gj/dj); + v3_add(&dzeta[k][j],&dzeta[k][j],&tmp); + v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj); + v3_add(&dzeta[k][j],&dzeta[k][j],&tmp); // k j /* increase k counter */ - exchange->kcount++; + exchange->kcnt+=1; + + return 0; +} + +/* first j loop within first i loop */ +int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { + + t_albe_mult_params *params; + t_albe_exchange *exchange; + + params=moldyn->pot_params; + exchange=&(params->exchange); + + /* increase j counter */ + exchange->jcnt+=1; return 0; } -int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { +/* second j loop within first i loop */ +int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; - t_3dvec force; - double f_a,df_a,b,db,f_c,df_c; - double f_r,df_r; - double scale; - double mu,B; - double lambda,A; - double d_ij,r0; - unsigned char brand; - double S,R,s_r,arg; + + int j; + double d,f_a,df_a,f_r,df_r,f_c,df_c,b,db; + double A,B,mu,lambda,r0; double energy; + t_3dvec *dist,force; + double scale; + u8 brand; params=moldyn->pot_params; exchange=&(params->exchange); + /* get j counter */ + j=exchange->j2cnt; + + /* skip if j not within cutoff */ + if(exchange->skip[j]) { + moldyn->run3bp=0; + exchange->j2cnt+=1; + return 0; + } + + /* distance */ + d=exchange->d[j]; + dist=&(exchange->dist[j]); + f_c=exchange->f_c[j]; + df_c=exchange->df_c[j]; + + /* determine parameters to calculate fa, dfa, fr, dfr */ brand=aj->brand; if(brand==ai->brand) { - S=params->S[brand]; - R=params->R[brand]; B=params->B[brand]; A=params->A[brand]; r0=params->r0[brand]; @@ -304,8 +371,6 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { lambda=params->lambda[brand]; } else { - S=params->Smixed; - R=params->Rmixed; B=params->Bmixed; A=params->Amixed; r0=params->r0_mixed; @@ -313,212 +378,101 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { lambda=params->lambda_m; } - d_ij=exchange->d_ij; - - /* f_c, df_c */ - if(d_ijzeta_ij==0.0) { - b=1.0; - db=0.0; - } - else { - b=1.0/sqrt(1.0+exchange->zeta_ij); - db=-0.5*b/(1.0+exchange->zeta_ij); - } + b=1.0/sqrt(1.0+exchange->zeta[j]); + db=-0.5*b/(1.0+exchange->zeta[j]); + + /* energy contribution */ + energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism + moldyn->energy+=energy; + ai->e+=energy; - /* force contribution for atom i */ + /* force contribution for atom i due to ij bond */ scale=-0.5*(f_c*(df_r-b*df_a)+df_c*(f_r-b*f_a)); // - in albe formalism - v3_scale(&force,&(exchange->dist_ij),scale); + v3_scale(&force,dist,scale); v3_add(&(ai->f),&(ai->f),&force); - /* force contribution for atom j */ + /* force contribution for atom j due to ij bond */ v3_scale(&force,&force,-1.0); // dri rij = - drj rij v3_add(&(aj->f),&(aj->f),&force); /* virial */ - virial_calc(aj,&force,&(exchange->dist_ij)); - -#ifdef DEBUG -if(moldyn->time>DSTART&&moldyn->timeatom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) { - 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])) - printf(" total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z); - if(aj==&(moldyn->atom[0])) - 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_calc(aj,&force,dist); /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */ - exchange->pre_dzeta=0.5*f_a*f_c*db; + exchange->pre_dzeta=0.5*f_a*exchange->f_c[j]*db; - /* energy contribution */ - energy=0.5*f_c*(f_r-b*f_a); // - in albe formalism - moldyn->energy+=energy; - ai->e+=energy; + /* force contribution (drj derivative) */ + v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta); + v3_add(&(aj->f),&(aj->f),&force); + v3_scale(&force,&force,-1.0); + v3_add(&(ai->f),&(ai->f),&force); + + /* virial */ + virial_calc(ai,&force,dist); /* reset k counter for second k loop */ - exchange->kcount=0; + exchange->kcnt=0; return 0; } -/* 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) { +/* second k loop within second j loop within first i loop */ +int albe_mult_i0_j2_k0(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { t_albe_mult_params *params; t_albe_exchange *exchange; - int kcount; - t_3dvec dist_ik,dist_ij; - double d_ik2,d_ik,d_ij2,d_ij; - unsigned char brand; - double S2; - double g,dg,cos_theta; - double pre_dzeta; - double f_c_ik,df_c_ik; - double dijdik_inv,fcdg,dfcg; - t_3dvec dcosdrj,dcosdrk; - t_3dvec force,tmp; + + int j,k; + t_3dvec force; params=moldyn->pot_params; exchange=&(params->exchange); - kcount=exchange->kcount; - - if(kcount>ALBE_MAXN) - printf("FATAL: neighbours!\n"); - - /* d_ik2 */ - d_ik2=exchange->d_ik2[kcount]; - - brand=ak->brand; - if(brand==ai->brand) - S2=params->S2[brand]; - else - S2=params->S2mixed; - /* return if d_ik > S */ - if(d_ik2>S2) { - exchange->kcount++; + /* k!=j & check whether to run k */ + k=exchange->kcnt; + j=exchange->j2cnt; + if((k==j)|(exchange->skip[k])) { + exchange->kcnt+=1; return 0; } - - /* prefactor dzeta */ - pre_dzeta=exchange->pre_dzeta; - - /* dist_ik, d_ik */ - dist_ik=exchange->dist_ik[kcount]; - d_ik=exchange->d_ik[kcount]; - - /* f_c_ik, df_c_ik */ - f_c_ik=exchange->f_c_ik[kcount]; - df_c_ik=exchange->df_c_ik[kcount]; - - /* dist_ij, d_ij2, d_ij */ - dist_ij=exchange->dist_ij; - d_ij2=exchange->d_ij2; - d_ij=exchange->d_ij; - - /* g, dg, cos_theta */ - g=exchange->g[kcount]; - dg=exchange->dg[kcount]; - cos_theta=exchange->cos_theta[kcount]; - - /* cos_theta derivatives wrt j,k */ - dijdik_inv=1.0/(d_ij*d_ik); - 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); // k - v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2); - v3_add(&dcosdrk,&dcosdrk,&tmp); - - /* f_c_ik * dg, df_c_ik * g */ - fcdg=f_c_ik*dg; - dfcg=df_c_ik*g; - - /* derivative wrt j */ - v3_scale(&force,&dcosdrj,fcdg*pre_dzeta); - - /* force contribution */ - v3_add(&(aj->f),&(aj->f),&force); - -#ifdef DEBUG -if(moldyn->time>DSTART&&moldyn->timeatom[DATOM])) { - 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); - printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI)); - printf(" d ij ik = %f %f\n",d_ij,d_ik); - } -} -#endif - - /* force contribution to atom i */ + + /* force contribution (drk derivative) */ + v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta); + v3_add(&(ak->f),&(ak->f),&force); v3_scale(&force,&force,-1.0); v3_add(&(ai->f),&(ai->f),&force); /* virial */ - virial_calc(ai,&force,&dist_ij); - - /* derivative wrt k */ - v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik - v3_scale(&tmp,&dcosdrk,fcdg); - v3_add(&force,&force,&tmp); - v3_scale(&force,&force,pre_dzeta); + virial_calc(ai,&force,&(exchange->dist[k])); - /* force contribution */ - v3_add(&(ak->f),&(ak->f),&force); + /* increase k counter */ + exchange->kcnt+=1; -#ifdef DEBUG -if(moldyn->time>DSTART&&moldyn->timeatom[DATOM])) { - 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); - printf(" angle: %f\n",acos(cos_theta)*360.0/(2*M_PI)); - printf(" d ij ik = %f %f\n",d_ij,d_ik); - } + return 0; } -#endif - /* force contribution to atom i */ - v3_scale(&force,&force,-1.0); - v3_add(&(ai->f),&(ai->f),&force); +int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { - /* virial */ - virial_calc(ai,&force,&dist_ik); - - /* increase k counter */ - exchange->kcount++; + t_albe_mult_params *params; + t_albe_exchange *exchange; - return 0; + params=moldyn->pot_params; + exchange=&(params->exchange); + /* increase j counter */ + exchange->j2cnt+=1; + + return 0; } int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc) { diff --git a/potentials/albe.h b/potentials/albe.h index 15386e7..0e4d4ee 100644 --- a/potentials/albe.h +++ b/potentials/albe.h @@ -13,34 +13,31 @@ /* albe exchange type */ typedef struct s_albe_exchange { - t_3dvec dist_ij; - double d_ij2; - double d_ij; + t_3dvec dist[ALBE_MAXN]; + double d2[ALBE_MAXN]; + double d[ALBE_MAXN]; - t_3dvec dist_ik[ALBE_MAXN]; - double d_ik2[ALBE_MAXN]; - double d_ik[ALBE_MAXN]; + double f_c[ALBE_MAXN]; + double df_c[ALBE_MAXN]; - double f_c_ik[ALBE_MAXN]; - double df_c_ik[ALBE_MAXN]; + double zeta[ALBE_MAXN]; + t_3dvec dzeta[ALBE_MAXN][ALBE_MAXN]; - double g[ALBE_MAXN]; - double dg[ALBE_MAXN]; - double cos_theta[ALBE_MAXN]; + u8 skip[ALBE_MAXN]; - double *gamma_i; - double *c_i; - double *d_i; - double *h_i; + double *gamma_[ALBE_MAXN]; + double *c_[ALBE_MAXN]; + double *d_[ALBE_MAXN]; + double c2_[ALBE_MAXN]; + double d2_[ALBE_MAXN]; + double c2d2_[ALBE_MAXN]; + double *h_[ALBE_MAXN]; - double ci2; - double di2; - double ci2di2; - - double zeta_ij; double pre_dzeta; - int kcount; + int jcnt; + int j2cnt; + int kcnt; } t_albe_exchange; /* albe mult (2!) potential parameters */ @@ -76,12 +73,15 @@ typedef struct s_albe_mult_params { /* function prototypes */ int albe_mult_set_params(t_moldyn *moldyn,int element1,int elemnt2); -int albe_mult_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); -int albe_mult_3bp_k1(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); -int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); -int albe_mult_3bp_k2(t_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); +int albe_mult_i0(t_moldyn *moldyn,t_atom *ai); +int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); +int albe_mult_i0_j0_k0(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); +int albe_mult_i0_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); +int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); +int albe_mult_i0_j2_k0(t_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bc); +int albe_mult_i0_j3(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); int albe_mult_check_2b_bond(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,u8 bc); /* albe potential parameter defines */ -- 2.20.1