- /*
- * indirect neighbour cells of atom i
- */
- for(j=c;j<27;j++) {
- this=&(neighbour[j]);
- list_reset(this);
- if(this->start!=NULL) {
-
- do {
- btom=this->current->data;
-
- /* 2 body stuff */
-
-
- /* determine cell neighbours of btom */
- ki=(btom->r.x+(moldyn->dim.x/2))/lc->x;
- kj=(btom->r.y+(moldyn->dim.y/2))/lc->y;
- kk=(btom->r.z+(moldyn->dim.z/2))/lc->z;
- ck=link_cell_neighbour_index(moldyn,ki,kj,kk,
- neighbourk);
-
- /* cell of btom */
- thisk=&(neighbourk[0]);
- list_reset(thisk);
- do {
- ktom=thisk->current->data;
- if(ktom==btom)
- continue;
- if(ktom==&(atom[i]))
- continue;
-
- /* 3 body stuff (1) */
+ t_tersoff_mult_params *params;
+ t_tersoff_exchange *exchange;
+
+ t_3dvec force,temp;
+ t_3dvec *dist_ij;
+ double b,db,tmp;
+ double f_c,df_c,f_a,df_a;
+ double chi,ni,betaini,nj,betajnj;
+ double zeta;
+
+ params=moldyn->pot2b_params;
+ exchange=&(params->exchange);
+
+ /* we do not run if f_c_ij was detected to be 0! */
+ if(!(exchange->run2bp_post))
+ return 0;
+
+ f_c=exchange->f_c;
+ df_c=exchange->df_c;
+ f_a=exchange->f_a;
+ df_a=exchange->df_a;
+ betaini=exchange->betaini;
+ betajnj=exchange->betajnj;
+ ni=*(exchange->n_i);
+ nj=*(exchange->n_j);
+ chi=exchange->chi;
+ dist_ij=&(exchange->dist_ij);
+
+ /* Vij and dVij */
+ zeta=exchange->zeta_ij;
+ tmp=betaini*pow(zeta,ni-1.0); /* beta^n * zeta^n-1 */
+ b=(1+zeta*tmp); /* 1 + beta^n * zeta^n */
+ db=chi*pow(b,-1.0/(2*ni)-1); /* chi * (...)^(-1/2n - 1) */
+ b=db*b; /* b_ij */
+ db*=-0.5*tmp; /* db_ij */
+ v3_scale(&force,&(exchange->dzeta_ij),f_a*db);
+ v3_scale(&temp,dist_ij,df_a*b);
+ v3_add(&force,&force,&temp);
+ v3_scale(&force,&force,f_c);
+ v3_scale(&temp,dist_ij,df_c*b*f_a);
+ v3_add(&force,&force,&temp);
+
+ /* add energy of 3bp sum */
+ moldyn->energy+=(0.5*f_c*b*f_a);
+
+ /* add force (sub, as F = - dVij) */
+ v3_sub(&(ai->f),&(ai->f),&force);
+
+ /* dVji */
+ zeta=exchange->zeta_ji;
+ tmp=betajnj*pow(zeta,nj-1.0); /* beta^n * zeta^n-1 */
+ b=(1+zeta*tmp); /* 1 + beta^n * zeta^n */
+ db=chi*pow(b,-1.0/(2*nj)-1); /* chi * (...)^(-1/2n - 1) */
+ b=db*b; /* b_ij */
+ db*=-0.5*tmp; /* db_ij */
+ v3_scale(&force,&(exchange->dzeta_ji),f_a*db);
+ v3_scale(&temp,dist_ij,df_a*b);
+ v3_add(&force,&force,&temp);
+ v3_scale(&force,&force,f_c);
+ v3_scale(&temp,dist_ij,df_c*b*f_a);
+ v3_add(&force,&force,&temp);
+
+ /* add force (sub, as F = - dVji) */
+ v3_sub(&(ai->f),&(ai->f),&force);