-/* create mixed terms from parameters and set them */
-int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
-
- printf("[moldyn] tersoff parameter completion\n");
- p->Smixed=sqrt(p->S[0]*p->S[1]);
- p->Rmixed=sqrt(p->R[0]*p->R[1]);
- p->Amixed=sqrt(p->A[0]*p->A[1]);
- p->Bmixed=sqrt(p->B[0]*p->B[1]);
- p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
- p->mu_m=0.5*(p->mu[0]+p->mu[1]);
-
- printf("[moldyn] tersoff mult parameter info:\n");
- printf(" S (m) | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed);
- printf(" R (m) | %.12f | %.12f | %.12f\n",p->R[0],p->R[1],p->Rmixed);
- printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
- printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
- 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(" beta | %.10f | %.10f\n",p->beta[0],p->beta[1]);
- printf(" n | %f | %f\n",p->n[0],p->n[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(" chi | %f \n",p->chi);
-
- return 0;
-}
-
-/* tersoff 1 body part */
-int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
-
- int num;
- t_tersoff_mult_params *params;
- t_tersoff_exchange *exchange;
-
- num=ai->bnum;
- params=moldyn->pot1b_params;
- exchange=&(params->exchange);
-
- /*
- * simple: point constant parameters only depending on atom i to
- * their right values
- */
-
- exchange->beta_i=&(params->beta[num]);
- exchange->n_i=&(params->n[num]);
- exchange->c_i=&(params->c[num]);
- exchange->d_i=&(params->d[num]);
- exchange->h_i=&(params->h[num]);
-
- exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
- exchange->ci2=params->c[num]*params->c[num];
- exchange->di2=params->d[num]*params->d[num];
- exchange->ci2di2=exchange->ci2/exchange->di2;
-
- return 0;
-}
-
-/* tersoff 2 body part */
-int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
-
- t_tersoff_mult_params *params;
- t_tersoff_exchange *exchange;
- t_3dvec dist_ij,force;
- double d_ij;
- double A,B,R,S,lambda,mu;
- double f_r,df_r;
- double f_c,df_c;
- int num;
- double s_r;
- double arg;
- double scale;
-
- params=moldyn->pot2b_params;
- num=aj->bnum;
- exchange=&(params->exchange);
-
- /* clear 3bp and 2bp post run */
- exchange->run3bp_ij=0;
- exchange->run3bp_ji=0;
- exchange->run3bp_jk=0;
- exchange->run2bp_post_ij=0;
- exchange->run2bp_post_ji=0;
- exchange->run2bp_post_jk=0;
-
- /*
- * calc of 2bp contribution of V_ij and dV_ij/ji
- *
- * for Vij and dV_ij we need:
- * - f_c_ij, df_c_ij
- * - f_r_ij, df_r_ij
- *
- * for dV_ji we need:
- * - f_c_ji = f_c_ij, df_c_ji = df_c_ij
- * - f_r_ji = f_r_ij; df_r_ji = df_r_ij
- *
- */
-
- /* dist_ij, d_ij */
- v3_sub(&dist_ij,&(aj->r),&(ai->r));
- if(bc) check_per_bound(moldyn,&dist_ij);
- d_ij=v3_norm(&dist_ij);
-
- /* save for use in 3bp */
- exchange->d_ij=d_ij;
- exchange->dist_ij=dist_ij;
-
- /* constants */
- if(num==ai->bnum) {
- S=params->S[num];
- R=params->R[num];
- A=params->A[num];
- B=params->B[num];
- lambda=params->lambda[num];
- mu=params->mu[num];
- exchange->chi=1.0;
- }
- else {
- S=params->Smixed;
- R=params->Rmixed;
- A=params->Amixed;
- B=params->Bmixed;
- lambda=params->lambda_m;
- mu=params->mu_m;
- params->exchange.chi=params->chi;
- }
-
- /* if d_ij > S => no force & potential energy contribution */
- if(d_ij>S)
- return 0;
-
- /* more constants */
- exchange->beta_j=&(params->beta[num]);
- exchange->n_j=&(params->n[num]);
- exchange->c_j=&(params->c[num]);
- exchange->d_j=&(params->d[num]);
- exchange->h_j=&(params->h[num]);
- if(num==ai->bnum) {
- exchange->betajnj=exchange->betaini;
- exchange->cj2=exchange->ci2;
- exchange->dj2=exchange->di2;
- exchange->cj2dj2=exchange->ci2di2;
- }
- else {
- exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
- exchange->cj2=params->c[num]*params->c[num];
- exchange->dj2=params->d[num]*params->d[num];
- exchange->cj2dj2=exchange->cj2/exchange->dj2;
- }
-
- /* f_r_ij = f_r_ji, df_r_ij = df_r_ji */
- f_r=A*exp(-lambda*d_ij);
- df_r=-lambda*f_r/d_ij;
-
- /* f_a, df_a calc (again, same for ij and ji) | save for later use! */
- exchange->f_a=-B*exp(-mu*d_ij);
- exchange->df_a=-mu*exchange->f_a/d_ij;
-
- /* f_c, df_c calc (again, same for ij and ji) */
- if(d_ij<R) {
- /* f_c = 1, df_c = 0 */
- f_c=1.0;
- df_c=0.0;
- /* two body contribution (ij, ji) */
- v3_scale(&force,&dist_ij,-df_r);
- }
- else {
- s_r=S-R;
- arg=M_PI*(d_ij-R)/s_r;
- f_c=0.5+0.5*cos(arg);
- df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij));
- /* two body contribution (ij, ji) */
- v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c);
- }
-
- /* add forces of 2bp (ij, ji) contribution
- * dVij = dVji and we sum up both: no 1/2) */
- v3_add(&(ai->f),&(ai->f),&force);
-
- /* energy 2bp contribution (ij, ji) is 0.5 f_r f_c ... */
- moldyn->energy+=(0.5*f_r*f_c);
-
- /* save for use in 3bp */
- exchange->f_c=f_c;
- exchange->df_c=df_c;
-
- /* enable the run of 3bp function and 2bp post processing */
- exchange->run3bp=1;
- exchange->run2bp_post=1;
-
- /* reset 3bp sums */
- exchange->zeta_ij=0.0;
- exchange->zeta_ji=0.0;
- exchange->zeta_kl=0.0;
- v3_zero(&(exchange->db_ij));
- v3_zero(&(exchange->db_ji));
- v3_zero(&(exchange->db_jk));
-
- return 0;
-}
-
-/* tersoff 2 body post part */
-
-int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
-
- /* here we have to allow for the 3bp sums */
-
- t_tersoff_mult_params *params;
- t_tersoff_exchange *exchange;
-
- t_3dvec force,temp,*db_ij,*dist_ij;
- double db_ij_scale1,db_ij_scale2;
- double b_ij;
- double f_c,df_c,f_a,df_a;
- double chi,n,n_betan;
- 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;
-
- db_ij=&(exchange->db_ij);
- f_c=exchange->f_c;
- df_c=exchange->df_c;
- f_a=exchange->f_a;
- df_a=exchange->df_a;
- n_betan=exchange->n_betan;
- n=*(exchange->n);
- chi=exchange->chi;
- dist_ij=&(exchange->dist_ij);
- zeta=exchange->zeta;
-
- db_ij_scale2=pow(zeta,n-1.0);
-printf("DEBUG: %.15f %.15f\n",zeta,db_ij_scale2);
- db_ij_scale1=db_ij_scale2*zeta;
- db_ij_scale2*=n_betan;
- db_ij_scale1=pow((1.0+n_betan*db_ij_scale1),-1.0/(2*n)-1);
- b_ij=chi*db_ij_scale1*(1.0+n_betan*db_ij_scale1);
- db_ij_scale1*=(-1.0*chi/(2*n));
-
- /* db_ij part */
- v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2));
- v3_scale(db_ij,db_ij,f_a);
-
- /* df_a part */
- v3_scale(&temp,dist_ij,b_ij*df_a);
-
- /* db_ij + df_a part */
- v3_add(&force,&temp,db_ij);
- v3_scale(&force,&force,f_c);
-
- /* df_c part */
- v3_scale(&temp,dist_ij,f_a*b_ij*df_c);
-
- /* add energy of 3bp sum */
- moldyn->energy+=(0.5*f_c*b_ij*f_a);
-
- /* add force of 3bp calculation (all three parts) */
- v3_add(&(ai->f),&temp,&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;
- double R,S,s_r;
- double d_ij,d_ik,d_jk;
- double rxxryy,dxxdyy;
- double f_c,df_c,f_a,df_a;
- double f_c_ik,df_c_ik,arg;
- double n,c,d,h;
- double c2,d2,c2d2;
- double cos_theta,d_costheta1,d_costheta2;
- double h_cos,d2_h_cos2;
- double frac;
- double g;
- 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) */