From 2b6a06cdf69b0d8f1dafc0e21b89c6bb4bfc192c Mon Sep 17 00:00:00 2001 From: hackbard Date: Thu, 24 Jul 2008 13:55:41 +0200 Subject: [PATCH] new albe potential + new force calc routine (old potentials need to be adapted!) --- moldyn.c | 170 +++++++++++++++++++++------------------------- moldyn.h | 24 +++---- potentials/albe.c | 151 ++++++++++++++++++++++++++++++---------- potentials/albe.h | 14 ++-- 4 files changed, 216 insertions(+), 143 deletions(-) diff --git a/moldyn.c b/moldyn.c index bf7b9f0..1aff5a8 100644 --- a/moldyn.c +++ b/moldyn.c @@ -215,11 +215,11 @@ int set_potential(t_moldyn *moldyn,u8 type) { switch(type) { case MOLDYN_POTENTIAL_TM: - moldyn->func1b=tersoff_mult_1bp; - moldyn->func3b_j1=tersoff_mult_3bp_j1; - moldyn->func3b_k1=tersoff_mult_3bp_k1; - moldyn->func3b_j2=tersoff_mult_3bp_j2; - moldyn->func3b_k2=tersoff_mult_3bp_k2; + moldyn->func_i0=tersoff_mult_1bp; + moldyn->func_j1=tersoff_mult_3bp_j1; + moldyn->func_j1_k0=tersoff_mult_3bp_k1; + moldyn->func_j1c=tersoff_mult_3bp_j2; + moldyn->func_j1_k1=tersoff_mult_3bp_k2; moldyn->check_2b_bond=tersoff_mult_check_2b_bond; break; /* @@ -232,21 +232,21 @@ int set_potential(t_moldyn *moldyn,u8 type) { break; */ case MOLDYN_POTENTIAL_AM: - moldyn->func1b=albe_mult_i0; - moldyn->func2b=albe_mult_i0_j0; - moldyn->func3b_0=albe_mult_i0_j0_k0; - moldyn->func3b_1=albe_mult_i0_j1; - moldyn->func3b_j1=albe_mult_i0_j2; - moldyn->func3b_k2=albe_mult_i0_j2_k0; - moldyn->func3b_j2=albe_mult_i0_j3; + moldyn->func_i0=albe_mult_i0; + moldyn->func_j0=albe_mult_i0_j0; + moldyn->func_j0_k0=albe_mult_i0_j0_k0; + moldyn->func_j0e=albe_mult_i0_j1; + moldyn->func_j1=albe_mult_i0_j2; + moldyn->func_j1_k0=albe_mult_i0_j2_k0; + moldyn->func_j1c=albe_mult_i0_j3; moldyn->check_2b_bond=albe_mult_check_2b_bond; break; case MOLDYN_POTENTIAL_HO: - moldyn->func2b=harmonic_oscillator; + moldyn->func_j0=harmonic_oscillator; moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond; break; case MOLDYN_POTENTIAL_LJ: - moldyn->func2b=lennard_jones; + moldyn->func_j0=lennard_jones; moldyn->check_2b_bond=lennard_jones_check_2b_bond; break; default: @@ -1860,8 +1860,8 @@ int potential_force_calc(t_moldyn *moldyn) { /* single particle potential/force */ if(itom[i].attr&ATOM_ATTR_1BP) - if(moldyn->func1b) - moldyn->func1b(moldyn,&(itom[i])); + if(moldyn->func_i0) + moldyn->func_i0(moldyn,&(itom[i])); if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP))) continue; @@ -1876,8 +1876,14 @@ int potential_force_calc(t_moldyn *moldyn) { dnlc=lc->dnlc; +#ifndef STATIC_LISTS + /* check for later 3 body interaction */ + if(itom[i].attr&ATOM_ATTR_3BP) + memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); +#endif + /* first loop over atoms j */ - if(moldyn->func2b) { + if(moldyn->func_j0) { for(j=0;j<27;j++) { bc_ij=(jattr&ATOM_ATTR_2BP)& - (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); - } - } #else this=&(neighbour_i[j]); list_reset_f(this); @@ -1909,6 +1903,7 @@ int potential_force_calc(t_moldyn *moldyn) { do { jtom=this->current->data; +#endif if(jtom==&(itom[i])) continue; @@ -1918,13 +1913,12 @@ int potential_force_calc(t_moldyn *moldyn) { if((jtom->attr&ATOM_ATTR_2BP)& (itom[i].attr&ATOM_ATTR_2BP)) { - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); + moldyn->func_j0(moldyn, + &(itom[i]), + jtom, + bc_ij); } -// REWRITE ONCE WORKING!!! /* 3 body potential/force */ /* in j loop, 3bp run can be skipped */ @@ -1937,9 +1931,10 @@ int potential_force_calc(t_moldyn *moldyn) { if(!(jtom->attr&ATOM_ATTR_3BP)) continue; - if(moldyn->func3b_0==NULL) + if(moldyn->func_j0_k0==NULL) continue; + /* first loop over atoms k in first j loop */ for(k=0;k<27;k++) { bc_ik=(kattr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; - moldyn->func3b_0(moldyn, - &(itom[i]), - jtom, - ktom, - bc_ik|bc_ij); + moldyn->func_j0_k0(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); #ifdef STATIC_LISTS } #else @@ -1984,35 +1979,27 @@ int potential_force_calc(t_moldyn *moldyn) { } - - if(moldyn->func3b_1) - moldyn->func3b_1(moldyn, + /* finish of first j loop after first k loop */ + if(moldyn->func_j0e) + moldyn->func_j0e(moldyn, &(itom[i]), jtom, bc_ij); -// UNTIL HERE + change function names! - - - +#ifdef STATIC_LISTS + } +#else } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); #endif } } - /* 3 body potential/force */ + /* continued 3 body potential/force */ if(!(itom[i].attr&ATOM_ATTR_3BP)) continue; - /* copy the neighbour lists */ -#ifdef STATIC_LISTS - /* no copy needed for static lists */ -#else - memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list)); -#endif - /* second loop over atoms j */ for(j=0;j<27;j++) { @@ -2045,18 +2032,18 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset 3bp run */ moldyn->run3bp=1; - if(moldyn->func3b_j1) - moldyn->func3b_j1(moldyn, - &(itom[i]), - jtom, - bc_ij); + if(moldyn->func_j1) + moldyn->func_j1(moldyn, + &(itom[i]), + jtom, + bc_ij); - /* in first j loop, 3bp run can be skipped */ + /* in j loop, 3bp run can be skipped */ if(!(moldyn->run3bp)) continue; - /* first loop over atoms k */ - if(moldyn->func3b_k1) { + /* first loop over atoms k in second j loop */ + if(moldyn->func_j1_k0) { for(k=0;k<27;k++) { @@ -2082,17 +2069,17 @@ int potential_force_calc(t_moldyn *moldyn) { if(!(ktom->attr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; - moldyn->func3b_k1(moldyn, - &(itom[i]), - jtom, - ktom, - bc_ik|bc_ij); + moldyn->func_j1_k0(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); #ifdef STATIC_LISTS } #else @@ -2104,14 +2091,15 @@ int potential_force_calc(t_moldyn *moldyn) { } - if(moldyn->func3b_j2) - moldyn->func3b_j2(moldyn, - &(itom[i]), - jtom, - bc_ij); + /* continued j loop after first k loop */ + if(moldyn->func_j1c) + moldyn->func_j1c(moldyn, + &(itom[i]), + jtom, + bc_ij); /* second loop over atoms k */ - if(moldyn->func3b_k2) { + if(moldyn->func_j1_k1) { for(k=0;k<27;k++) { @@ -2137,17 +2125,17 @@ int potential_force_calc(t_moldyn *moldyn) { if(!(ktom->attr&ATOM_ATTR_3BP)) continue; - if(ktom==jtom) - continue; + //if(ktom==jtom) + // continue; if(ktom==&(itom[i])) continue; - moldyn->func3b_k2(moldyn, - &(itom[i]), - jtom, - ktom, - bc_ik|bc_ij); + moldyn->func_j1_k1(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); #ifdef STATIC_LISTS } @@ -2160,11 +2148,11 @@ int potential_force_calc(t_moldyn *moldyn) { } - /* 2bp post function */ - if(moldyn->func3b_j3) { - moldyn->func3b_j3(moldyn, - &(itom[i]), - jtom,bc_ij); + /* finish of j loop after second k loop */ + if(moldyn->func_j1e) { + moldyn->func_j1e(moldyn, + &(itom[i]), + jtom,bc_ij); } #ifdef STATIC_LISTS } @@ -2173,7 +2161,7 @@ int potential_force_calc(t_moldyn *moldyn) { #endif } - + #ifdef DEBUG //printf("\n\n"); #endif diff --git a/moldyn.h b/moldyn.h index 6345f3f..066ad36 100644 --- a/moldyn.h +++ b/moldyn.h @@ -104,18 +104,18 @@ typedef struct s_moldyn { double volume; /* volume of sim cell (dim.x*dim.y*dim.z) */ /* potential force function and parameter pointers */ - int (*func1b)(struct s_moldyn *moldyn,t_atom *ai); - int (*func2b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_0)(struct s_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); - int (*func3b_1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j2)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j3)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_k1)(struct s_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); - int (*func3b_k2)(struct s_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_i0)(struct s_moldyn *moldyn,t_atom *ai); + int (*func_j0)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j0_k0)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j0e)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1_k0)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j1c)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1_k1)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j1e)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); void *pot_params; unsigned char run3bp; diff --git a/potentials/albe.c b/potentials/albe.c index 9dfba13..e8aee74 100644 --- a/potentials/albe.c +++ b/potentials/albe.c @@ -53,7 +53,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) { p->mu[0]=ALBE_MU_SI; p->gamma[0]=ALBE_GAMMA_SI; p->c[0]=ALBE_C_SI; + p->c2[0]=p->c[0]*p->c[0]; p->d[0]=ALBE_D_SI; + p->d2[0]=p->d[0]*p->d[0]; + p->c2d2[0]=p->c2[0]/p->d2[0]; p->h[0]=ALBE_H_SI; switch(element2) { case C: @@ -67,7 +70,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) { p->mu[1]=ALBE_MU_C; p->gamma[1]=ALBE_GAMMA_C; p->c[1]=ALBE_C_C; + p->c2[1]=p->c[1]*p->c[1]; p->d[1]=ALBE_D_C; + p->d2[1]=p->d[1]*p->d[1]; + p->c2d2[1]=p->c2[1]/p->d2[1]; p->h[1]=ALBE_H_C; /* mixed type: silicon carbide */ p->Smixed=ALBE_S_SIC; @@ -79,7 +85,10 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) { p->mu_m=ALBE_MU_SIC; p->gamma_m=ALBE_GAMMA_SIC; p->c_mixed=ALBE_C_SIC; + p->c2_mixed=p->c_mixed*p->c_mixed; p->d_mixed=ALBE_D_SIC; + p->d2_mixed=p->d_mixed*p->d_mixed; + p->c2d2_m=p->c2_mixed/p->d2_mixed; p->h_mixed=ALBE_H_SIC; break; default: @@ -105,10 +114,13 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) { printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1], p->lambda_m); printf(" mu | %f | %f | %f\n",p->mu[0],p->mu[1],p->mu_m); - printf(" gamma | %f | %f\n",p->gamma[0],p->gamma[1]); - printf(" c | %f | %f\n",p->c[0],p->c[1]); - printf(" d | %f | %f\n",p->d[0],p->d[1]); - printf(" h | %f | %f\n",p->h[0],p->h[1]); + printf(" gamma | %f | %f | %f\n",p->gamma[0],p->gamma[1],p->gamma_m); + printf(" c | %f | %f | %f\n",p->c[0],p->c[1],p->c_mixed); + printf(" d | %f | %f | %f\n",p->d[0],p->d[1],p->d_mixed); + printf(" c2 | %f | %f | %f\n",p->c2[0],p->c2[1],p->c2_mixed); + printf(" d2 | %f | %f | %f\n",p->d2[0],p->d2[1],p->d2_mixed); + printf(" c2d2 | %f | %f | %f\n",p->c2d2[0],p->c2d2[1],p->c2d2_m); + printf(" h | %f | %f | %f\n",p->h[0],p->h[1],p->h_mixed); return 0; } @@ -116,21 +128,18 @@ int albe_mult_set_params(t_moldyn *moldyn,int element1,int element2) { /* first i loop */ int albe_mult_i0(t_moldyn *moldyn,t_atom *ai) { - int i; t_albe_mult_params *params; t_albe_exchange *exchange; + int i; + params=moldyn->pot_params; exchange=&(params->exchange); /* 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)); + memset(exchange->dzeta[i],0,ALBE_MAXN*sizeof(t_3dvec)); exchange->jcnt=0; exchange->j2cnt=0; @@ -142,6 +151,7 @@ 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; @@ -164,8 +174,8 @@ int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* dist_ij, d_ij2 */ v3_sub(&dist,&(aj->r),&(ai->r)); - exchange->dist[j]=dist; if(bc) check_per_bound(moldyn,&dist); + exchange->dist[j]=dist; d2=v3_absolute_square(&dist); exchange->d2[j]=d2; @@ -176,16 +186,20 @@ int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { exchange->jcnt+=1; return 0; } + exchange->skip[j]=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 */ + /* 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]); + exchange->c2_[j]=&(params->c2[brand]); + exchange->d2_[j]=&(params->d2[brand]); + exchange->c2d2_[j]=&(params->c2d2[brand]); } else { R=params->Rmixed; @@ -195,13 +209,13 @@ int albe_mult_i0_j0(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { exchange->c_[j]=&(params->c_mixed); exchange->d_[j]=&(params->d_mixed); exchange->h_[j]=&(params->h_mixed); + exchange->c2_[j]=&(params->c2_mixed); + exchange->d2_[j]=&(params->d2_mixed); + exchange->c2d2_[j]=&(params->c2d2_m); } - 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=sqrt(exchange->d2[j]); + d=sqrt(d2); exchange->d[j]=d; /* f_c, df_c */ @@ -228,16 +242,23 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn, t_albe_mult_params *params; t_albe_exchange *exchange; + 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; + t_3dvec dcosdrj,dcosdrk,tmp; + t_3dvec *dzjj,*dzkk,*dzjk,*dzkj; params=moldyn->pot_params; exchange=&(params->exchange); + if(aj==ak) { + exchange->kcnt+=1; + return 0; + } + /* kjcnt; k=exchange->kcnt; @@ -249,7 +270,7 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn, exchange->kcnt+=1; return 0; } - + /* distances */ distj=exchange->dist[j]; distk=exchange->dist[k]; @@ -260,11 +281,11 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn, /* cos theta */ cos_theta=v3_scalar_product(&distj,&distk)*djdk_inv; - /* g(cos(theta)) - in albe: ik-depending values! */ + /* g(cos(theta)) ij and ik 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; + 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) { @@ -273,14 +294,19 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn, } else { 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; + 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; } - /* zeta */ + /* zeta - for albe: ik depending g function */ +//if(ai->tag==0) { +// printf("------> %.15f %.15f\n",dj,dk); +// printf("------> %.15f %.15f\n",dj,dk); +//} + exchange->zeta[j]+=(exchange->f_c[k]*gk); exchange->zeta[k]+=(exchange->f_c[j]*gj); @@ -293,19 +319,22 @@ int albe_mult_i0_j0_k0(t_moldyn *moldyn, v3_add(&dcosdrk,&dcosdrk,&tmp); /* zeta derivatives */ - dzeta=exchange->dzeta; + dzjj=&(exchange->dzeta[j][j]); + dzkk=&(exchange->dzeta[k][k]); + dzjk=&(exchange->dzeta[j][k]); + dzkj=&(exchange->dzeta[k][j]); v3_scale(&tmp,&dcosdrj,exchange->f_c[k]*dgk); - v3_add(&dzeta[j][j],&dzeta[j][j],&tmp); // j j + v3_add(dzjj,dzjj,&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_add(dzkk,dzkk,&tmp); // k k + v3_scale(&tmp,&distk,-exchange->df_c[k]*gk); // dk rik = - di rik + v3_add(dzjk,dzjk,&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_add(dzjk,dzjk,&tmp); // j k + v3_scale(&tmp,&distj,-exchange->df_c[j]*gj); // dj rij = - di rij + v3_add(dzkj,dzkj,&tmp); v3_scale(&tmp,&dcosdrj,exchange->f_c[j]*dgj); - v3_add(&dzeta[k][j],&dzeta[k][j],&tmp); // k j + v3_add(dzkj,dzkj,&tmp); // k j /* increase k counter */ exchange->kcnt+=1; @@ -400,22 +429,51 @@ int albe_mult_i0_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { v3_scale(&force,dist,scale); v3_add(&(ai->f),&(ai->f),&force); +#ifdef NDEBUG +if(ai->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d (ij) %.15f\n",force.x,force.y,force.z,ai->tag,aj->tag,exchange->zeta[j]); +printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z); +} +#endif + /* 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); +#ifdef NDEBUG +if(aj->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d (ji) %.15f\n",force.x,force.y,force.z,aj->tag,ai->tag,exchange->zeta[j]); +printf(" t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z); +} +#endif + /* virial */ 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*exchange->f_c[j]*db; + exchange->pre_dzeta=0.5*f_a*f_c*db; /* force contribution (drj derivative) */ v3_scale(&force,&(exchange->dzeta[j][j]),exchange->pre_dzeta); v3_add(&(aj->f),&(aj->f),&force); + +#ifdef NDEBUG +if(aj->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d (j der)\n",force.x,force.y,force.z,aj->tag,ai->tag); +printf(" t: %.15f %.15f %.15f\n",aj->f.x,aj->f.y,aj->f.z); +} +#endif + v3_scale(&force,&force,-1.0); v3_add(&(ai->f),&(ai->f),&force); +#ifdef NDEBUG +if(ai->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d (i contr j der)\n",force.x,force.y,force.z,ai->tag,aj->tag); +printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z); +} +#endif + /* virial */ virial_calc(ai,&force,dist); @@ -438,6 +496,11 @@ int albe_mult_i0_j2_k0(t_moldyn *moldyn, params=moldyn->pot_params; exchange=&(params->exchange); + if(aj==ak) { + exchange->kcnt+=1; + return 0; + } + /* k!=j & check whether to run k */ k=exchange->kcnt; j=exchange->j2cnt; @@ -449,9 +512,25 @@ int albe_mult_i0_j2_k0(t_moldyn *moldyn, /* force contribution (drk derivative) */ v3_scale(&force,&(exchange->dzeta[j][k]),exchange->pre_dzeta); v3_add(&(ak->f),&(ak->f),&force); + +#ifdef NDEBUG +if(ak->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d %d (k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag); +printf(" t: %.15f %.15f %.15f\n",ak->f.x,ak->f.y,ak->f.z); +} +#endif + v3_scale(&force,&force,-1.0); v3_add(&(ai->f),&(ai->f),&force); +#ifdef NDEBUG +if(ai->tag==0) { +printf("force: %.15f %.15f %.15f | %d %d %d -- %d(i contr k der)\n",force.x,force.y,force.z,ai->tag,aj->tag,ak->tag,k); +printf(" t: %.15f %.15f %.15f\n",ai->f.x,ai->f.y,ai->f.z); +printf(" ## %f\n",exchange->d[k]); +} +#endif + /* virial */ virial_calc(ai,&force,&(exchange->dist[k])); diff --git a/potentials/albe.h b/potentials/albe.h index 0e4d4ee..39c4df8 100644 --- a/potentials/albe.h +++ b/potentials/albe.h @@ -8,7 +8,7 @@ #ifndef ALBE_H #define ALBE_H -#define ALBE_MAXN 16*27 +#define ALBE_MAXN (4*27) /* albe exchange type */ typedef struct s_albe_exchange { @@ -28,9 +28,9 @@ typedef struct s_albe_exchange { 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 *c2_[ALBE_MAXN]; + double *d2_[ALBE_MAXN]; + double *c2d2_[ALBE_MAXN]; double *h_[ALBE_MAXN]; double pre_dzeta; @@ -62,9 +62,15 @@ typedef struct s_albe_mult_params { double gamma[2]; double gamma_m; double c[2]; + double c2[2]; double c_mixed; + double c2_mixed; double d[2]; + double d2[2]; double d_mixed; + double d2_mixed; + double c2d2[2]; + double c2d2_m; double h[2]; double h_mixed; -- 2.20.1