From 000e702c884e337f27ae66b6506ff9f9b613f8a8 Mon Sep 17 00:00:00 2001 From: hackbard Date: Sun, 3 Dec 2006 01:54:29 +0000 Subject: [PATCH] tests ... --- math/math.c | 5 +++++ math/math.h | 1 + moldyn.c | 49 +++++++++++++++++++++++++++++-------------------- moldyn.h | 3 --- 4 files changed, 35 insertions(+), 23 deletions(-) diff --git a/math/math.c b/math/math.c index 6f4345e..54d79fd 100644 --- a/math/math.c +++ b/math/math.c @@ -66,6 +66,11 @@ int v3_cmp(t_3dvec *a,t_3dvec *b) { return(memcmp(a,b,sizeof(t_3dvec))); } +double v3_scalar_product(t_3dvec *a,t_3dvec *b) { + + return(a->x*b->x+a->y*b->y+a->z*b->z); +} + double v3_absolute_square(t_3dvec *a) { return(a->x*a->x+a->y*a->y+a->z*a->z); diff --git a/math/math.h b/math/math.h index 77100f3..8bbcc6c 100644 --- a/math/math.h +++ b/math/math.h @@ -25,6 +25,7 @@ int v3_zero(t_3dvec *vec); int v3_set(t_3dvec *vec,double *ptr); int v3_copy(t_3dvec *trg,t_3dvec *src); int v3_cmp(t_3dvec *a,t_3dvec *b); +double v3_scalar_product(t_3dvec *a,t_3dvec *b); double v3_absolute_square(t_3dvec *a); double v3_norm(t_3dvec *a); int v3_per_bound(t_3dvec *a,t_3dvec *dim); diff --git a/moldyn.c b/moldyn.c index af0132d..f504e68 100644 --- a/moldyn.c +++ b/moldyn.c @@ -1139,18 +1139,17 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { d_ij=v3_norm(&dist_ij); /* save for use in 3bp */ - exchange->dist_ij=dist_ij; /* <- needed ? */ exchange->d_ij=d_ij; + exchange->dist_ij=dist_ij; + exchange->d_ij2=d_ij*d_ij; /* constants */ if(num==aj->bnum) { S=params->S[num]; R=params->R[num]; A=params->A[num]; + B=params->B[num]; lambda=params->lambda[num]; - /* more constants depending of atoms i and j, needed in 3bp */ - params->exchange.B=&(params->B[num]); - params->exchange.mu=&(params->mu[num]); mu=params->mu[num]; params->exchange.chi=1.0; } @@ -1158,10 +1157,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { S=params->Smixed; R=params->Rmixed; A=params->Amixed; + B=params->Bmixed; lambda=params->lambda_m; - /* more constants depending of atoms i and j, needed in 3bp */ - params->exchange.B=&(params->Bmixed); - params->exchange.mu=&(params->mu_m); mu=params->mu_m; params->exchange.chi=params->chi; } @@ -1261,7 +1258,6 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { * */ - v3_sub(&dist_ik,&(ai->r),&(ak->r)); if(bc) check_per_bound(moldyn,&dist_ik); d_ik=v3_norm(&dist_ik); @@ -1277,11 +1273,11 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { } /* calc of f_c_ik */ - if(d_ik>S) - return 0; - - if(d_ikS) { + f_c_ik=0.0; + df_c_ik=0.0; + } + else if(d_ikc); d=*(exchange->d); h=*(exchange->h); + chi=exchange->chi; c2=exchange->c2; d2=exchange->d2; c2d2=exchange->c2d2; numer=d_ij2+d_ik*d_ik-d_jk*d_jk; denom=2*d_ij*d_ik; - cos_theta=numer/denom; + //cos_theta=numer/denom; + cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik); sin_theta=sqrt(1.0-(cos_theta*cos_theta)); theta=acos(cos_theta); d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom); @@ -1324,11 +1322,22 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* some usefull expressions */ frac1=c2/(d2-h_cos2); bracket1=1+c2d2-frac1; - bracket2=f_c_ik*bracket1; - bracket2_n_1=pow(bracket2,n-1.0); - bracket2_n=bracket2_n_1*bracket2; - bracket3=1+betan*bracket2_n; + if(f_c_ik==0.0) { + bracket2=0.0; + bracket2_n_1=0.0; + bracket2_n=0.0; + bracket3=1.0; + printf("Foo -> 0: "); + } + else { + bracket2=f_c_ik*bracket1; + bracket2_n_1=pow(bracket2,n-1.0); + bracket2_n=bracket2_n_1*bracket2; + bracket3=1.0+betan*bracket2_n; + printf("Foo -> 1: "); + } bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0); +printf("THETA: %.15f %.15f\n",cos_theta,theta*180/(2*M_PI)); bracket3_pow=bracket3_pow_1*bracket3; /* now go on with calc of b_ij and derivation of b_ij */ @@ -1360,8 +1369,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* add forces */ v3_add(&(ai->f),&(ai->f),&force); - /* energy is 0.5 f_r f_c, but we will sum it up twice ... */ - moldyn->energy+=(0.25*f_a*b_ij*f_c); + /* energy is 0.5 f_r f_c */ + moldyn->energy+=(0.5*f_a*b_ij*f_c); return 0; } diff --git a/moldyn.h b/moldyn.h index 8b40d4c..813af68 100644 --- a/moldyn.h +++ b/moldyn.h @@ -182,9 +182,6 @@ typedef struct s_tersoff_exchange { double chi; - double *B; - double *mu; - double *beta; double *n; double *c; -- 2.39.2