- /* Vij and dVij */
- zeta=exchange->zeta_ij;
- if(zeta==0.0) {
- moldyn->debug++; /* just for debugging ... */
- db=0.0;
- b=chi;
- v3_scale(&force,dist_ij,df_a*b*f_c);
- }
- else {
- 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); /* x(...)^(-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);
- v3_scale(&force,&force,-0.5);
-
- /* add force */
- v3_add(&(ai->f),&(ai->f),&force);
-
- /* add energy of 3bp sum */
- moldyn->energy+=(0.5*f_c*b*f_a);
-
- /* dVji */
- zeta=exchange->zeta_ji;
- if(zeta==0.0) {
- moldyn->debug++;
- b=chi;
- v3_scale(&force,dist_ij,df_a*b*f_c);
- }
- else {
- 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); /* x(...)^(-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);
- v3_scale(&force,&force,-0.5);
-
- /* add force */
- v3_add(&(ai->f),&(ai->f),&force);
-
- return 0;
-}
-
-/* tersoff 3 body part */
-
-int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
-
- t_tersoff_mult_params *params;
- t_tersoff_exchange *exchange;
- t_3dvec dist_ij,dist_ik,dist_jk;
- t_3dvec temp1,temp2;
- t_3dvec *dzeta;
- double R,S,s_r;
- double B,mu;
- double d_ij,d_ik,d_jk;
- double rr,dd;
- double f_c,df_c;
- double f_c_ik,df_c_ik,arg;
- double f_c_jk;
- double n,c,d,h;
- double c2,d2,c2d2;
- double cos_theta,d_costheta1,d_costheta2;
- double h_cos,d2_h_cos2;
- double frac,g,zeta,chi;
- double tmp;
- int num;
-
- params=moldyn->pot3b_params;
- exchange=&(params->exchange);
-
- if(!(exchange->run3bp))
- return 0;
-
- /*
- * calc of 3bp contribution of V_ij and dV_ij/ji/jk &
- * 2bp contribution of dV_jk
- *
- * for Vij and dV_ij we still need:
- * - b_ij, db_ij (zeta_ij)
- * - f_c_ik, df_c_ik, constants_i, cos_theta_ijk, d_costheta_ijk
- *
- * for dV_ji we still need:
- * - b_ji, db_ji (zeta_ji)
- * - f_c_jk, d_c_jk, constants_j, cos_theta_jik, d_costheta_jik
- *
- * for dV_jk we need:
- * - f_c_jk
- * - f_a_jk
- * - db_jk (zeta_jk)
- * - f_c_ji, df_c_ji, constants_j, cos_theta_jki, d_costheta_jki
- *
- */
-
- /*
- * get exchange data
- */
-
- /* dist_ij, d_ij - this is < S_ij ! */
- dist_ij=exchange->dist_ij;
- d_ij=exchange->d_ij;
-
- /* f_c_ij, df_c_ij (same for ji) */
- f_c=exchange->f_c;
- df_c=exchange->df_c;
-
- /*
- * calculate unknown values now ...
- */
-
- /* V_ij and dV_ij stuff (in b_ij there is f_c_ik) */
-
- /* dist_ik, d_ik */
- v3_sub(&dist_ik,&(ak->r),&(ai->r));
- if(bc) check_per_bound(moldyn,&dist_ik);
- d_ik=v3_norm(&dist_ik);