added 2bp post function .. this is getting sick!

diff --git a/moldyn.c b/moldyn.c
index f504e68..52a716e 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -950,6 +950,13 @@ printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[
                                        }
                                        
                                } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+               
+                               /* 2bp post function */
+                               if(moldyn->func2b_post)
+                                       mlodyn->func2b_post(moldyn,
+                                                           &(itom[i]),
+                                                           jtom,bc_ij);
+
                        }
                }
        }
@@ -1072,6 +1079,21 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
        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;
 }
 
@@ -1200,6 +1222,57 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        /* enable the run of 3bp function */
        exchange->run3bp=1;
 
+       /* reset 3bp sums */
+       exchange->3bp_sum1=0.0;
+       exchange->3bp_sum2=0.0;
+
+       return 0;
+}
+
+/* tersoff 2 body post part */
+
+int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,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;
+       double db_ij_scale1,db_ij_scale2;
+       double b_ij;
+       double f_c,df_c,f_a,df_a;
+
+       params=moldyn->pot2b_params;
+       exchange=&(moldyn->exchange);
+
+       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;
+
+       db_ij_scale1=(1+betan*3bp_sum1);
+       db_ij_scale2=(n*betan*3bp_sum2);
+       help=pow(db_ij_scale1,-1.0/(2*n)-1);
+       b_ij=chi*db_ij_scale1*help;
+       db_ij_scale1=-chi/(2*n)*help;
+
+       v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2));
+       v3_scale(db_ij,db_ij,f_a);
+
+       v3_scale(&temp,dist_ij,b_ij*df_a);
+
+       v3_add(&force,&temp,db_ij);
+       v3_scale(&force,&force,f_c);
+
+       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 */      
+       v3_add(&(ai->f),&temp,&force);
+
        return 0;
 }
 
@@ -1242,9 +1315,13 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
        d_ij=exchange->d_ij;
        d_ij2=exchange->d_ij2;
+       dist_ij=exchange->dist_ij;
 
        f_a=params->exchange.f_a;
        df_a=params->exchange.df_a;
+
+       f_c=exchange->f_c;
+       df_c=exchange->df_c;
        
        /* d_ij is <= S, as we didn't return so far! */
 
@@ -1305,8 +1382,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
        numer=d_ij2+d_ik*d_ik-d_jk*d_jk;
        denom=2*d_ij*d_ik;
-       //cos_theta=numer/denom;
-       cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
+       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);
@@ -1317,10 +1394,10 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
 
        h_cos=(h-cos_theta);
        h_cos2=h_cos*h_cos;
-       d2_h_cos2=d2-h_cos2;
+       d2_h_cos2=d2+h_cos2;
 
        /* some usefull expressions */
-       frac1=c2/(d2-h_cos2);
+       frac1=c2/(d2_h_cos2);
        bracket1=1+c2d2-frac1;
        if(f_c_ik==0.0) {
                bracket2=0.0;
@@ -1337,8 +1414,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
                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;
+printf("%.15f %.15f %.15f\n",bracket2_n_1,bracket2_n);
 
        /* now go on with calc of b_ij and derivation of b_ij */
        b_ij=chi*bracket3_pow;

diff --git a/moldyn.h b/moldyn.h
index 813af68..41b711a 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -36,9 +36,9 @@ typedef struct s_atom {
 #define ATOM_ATTR_FP   0x01    /* fixed position (bulk material) */
 #define ATOM_ATTR_HB   0x02    /* coupled to heat bath (velocity scaling) */
 
-#define ATOM_ATTR_1BP  0x10    /* single paricle potential */
-#define ATOM_ATTR_2BP  0x20    /* pair potential */
-#define ATOM_ATTR_3BP  0x40    /* 3 body potential */ 
+#define ATOM_ATTR_1BP          0x10    /* single paricle potential */
+#define ATOM_ATTR_2BP          0x20    /* pair potential */
+#define ATOM_ATTR_3BP          0x40    /* 3 body potential */ 
 
 /* cell lists */
 typedef struct s_linkcell {
@@ -73,6 +73,8 @@ typedef struct s_moldyn {
        int (*func1b)(struct s_moldyn *moldyn,t_atom *ai);
        void *pot1b_params;
        int (*func2b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+       void (*func2b_post)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,
+                           u8 bc);
        void *pot2b_params;
        int (*func3b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,
                      u8 bck);
@@ -194,6 +196,10 @@ typedef struct s_tersoff_exchange {
        double betan;
 
        u8 run3bp;
+
+       t_3dvec db_ij;
+       double 3bp_sum1;
+       double 3bp_sum2;
 } t_tersoff_exchange;
 
 /* tersoff multi (2!) potential parameters */
@@ -285,7 +291,7 @@ typedef struct s_tersoff_mult_params {
 #define TM_BETA_SI             1.1000e-6
 #define TM_N_SI                        0.78734
 #define TM_C_SI                        1.0039e5
-#define TM_D_SI                        1.62170
+#define TM_D_SI                        16.217
 #define TM_H_SI                        (-0.59825)
 
 #define TM_R_C                 1.8e-10                         /* m */
@@ -311,6 +317,7 @@ typedef struct s_tersoff_mult_params {
 
 typedef int (*pf_func1b)(t_moldyn *,t_atom *ai);
 typedef int (*pf_func2b)(t_moldyn *,t_atom *,t_atom *,u8 bc);
+typedef int (*pf_func2b_post)(t_moldyn *,t_atom *,t_atom *,u8 bc);
 typedef int (*pf_func3b)(t_moldyn *,t_atom *,t_atom *,t_atom *,u8 bc);
 
 int moldyn_init(t_moldyn *moldyn,int argc,char **argv);
@@ -325,6 +332,7 @@ int set_nn_dist(t_moldyn *moldyn,double dist);
 int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z);
 int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params);
 int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params);
+int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params);
 int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params);
 
 int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer);

diff --git a/sic.c b/sic.c
index b10c681..3408895 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -83,6 +83,7 @@ int main(int argc,char **argv) {
        tp.n[1]=TM_N_C;
        tp.c[1]=TM_C_C;
        tp.d[1]=TM_D_C;
+       tp.h[1]=TM_H_C;
 
        tp.chi=TM_CHI_SIC;
 
@@ -94,7 +95,7 @@ int main(int argc,char **argv) {
 
        /* set (initial) dimensions of simulation volume */
        printf("[sic] setting dimensions\n");
-       set_dim(&md,4*LC_SI,4*LC_SI,4*LC_SI,TRUE);
+       set_dim(&md,3*LC_SI,3*LC_SI,3*LC_SI,TRUE);
 
        /* set periodic boundary conditions in all directions */
        printf("[sic] setting periodic boundary conditions\n");
@@ -105,7 +106,7 @@ int main(int argc,char **argv) {
        create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
                       //ATOM_ATTR_2BP|ATOM_ATTR_HB,
-                      0,4,4,4);
+                      0,3,3,3);
 
        /* setting a nearest neighbour distance for the moldyn checks */
        set_nn_dist(&md,sqrt(3.0)*LC_SI/4.0); /* diamond ! */
@@ -124,7 +125,7 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        printf("[sic] adding schedule\n");
-       moldyn_add_schedule(&md,100,1.0e-15);
+       moldyn_add_schedule(&md,10,1.0e-15);
 
        /* activate logging */
        printf("[sic] activate logging\n");