From e763f6e3072fbbbc5ecf7345b989d7d29929a38f Mon Sep 17 00:00:00 2001 From: hackbard Date: Thu, 17 Aug 2006 15:41:43 +0000 Subject: [PATCH] more .. --- moldyn.c | 50 +++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 45 insertions(+), 5 deletions(-) diff --git a/moldyn.c b/moldyn.c index 1b52a56..ad0f9a3 100644 --- a/moldyn.c +++ b/moldyn.c @@ -1005,6 +1005,9 @@ int tersoff(t_moldyn *moldyn) { lambda=; B=; mu=; + chi=; + beta=; + betaN=; if(d_ij<=R) { f_c=1.0; @@ -1032,6 +1035,10 @@ int tersoff(t_moldyn *moldyn) { ck=link_cell_neighbour_index(moldyn,ki,kj,kk, neighbourk); + /* go for zeta - 3 body stuff! */ + zeta=0.0; + d_ij2=d_ij*d_ij; + /* cell of btom */ thisk=&(neighbourk[0]); list_reset(thisk); @@ -1044,12 +1051,45 @@ int tersoff(t_moldyn *moldyn) { /* 3 body stuff (1) */ - theta_ijk=; - sin_theta=; - cos_theta=; - hi_cos=; - hi_cos_square=; + v3_sub(&dist_ik,ktom,&(atom[i])); + d_ik=v3_norm(&dist_ik); + if(d_ik<=Sik) { + + Rik=; + Sik=; + Aik=; + lambda_ik=; + Bik=; + mu_ik=; + omega_ik=; + c_i=; + d_i=; + h_i=; + + + if(d_ik<=Rik) { + f_cik=1.0; + df_cik=0.0; + } + else { + sik_rik=Sik-Rik; + arg1ik=PI*(d_ik-Rik)/sik_rik; + f_cik=0.5+0.5*cos(arg1ik); + df_cik=-0.5*sin(arg1ik)* \ + (PI/(sik_rik*d_ik)); + f_rik=Aik*exp(-lambda_ik*d_ik); + f_aik=-Bik*exp(-mu_ik*d_ik); + } + + v3_sub(&distance_jk,ktom,btom); + cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/\ + (2*d_ij*d_ik); + theta=arccos(cos_theta); + + } + else + continue; /* end 3 body stuff (1) */ -- 2.20.1