no segfaulst anymore BUT take care, definetly calculating bullshit

diff --git a/moldyn.c b/moldyn.c
index 52a716e..716a087 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -139,6 +139,14 @@ int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) {
        return 0;
 }
 
+int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) {
+
+       moldyn->func2b=func;
+       moldyn->pot2b_params=params;
+
+       return 0;
+}
+
 int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) {
 
        moldyn->func3b=func;
@@ -940,9 +948,9 @@ int potential_force_calc(t_moldyn *moldyn) {
                        if(ktom==&(itom[i]))
                                continue;
 
-printf("Debug: atom %d before 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
+//printf("Debug: atom %d before 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
                        moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
-printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
+//printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[i],jtom,ktom,itom[i].r.x,itom[i].f.x,itom[i].v.x);
 
                                                } while(list_next(that)!=\
                                                        L_NO_NEXT_ELEMENT);
@@ -953,7 +961,7 @@ printf("Debug: atom %d after 3bp: %08x %08x %08x | %.15f %.15f %.15f\n",i,&itom[
                
                                /* 2bp post function */
                                if(moldyn->func2b_post)
-                                       mlodyn->func2b_post(moldyn,
+                                       moldyn->func2b_post(moldyn,
                                                            &(itom[i]),
                                                            jtom,bc_ij);
 
@@ -1120,6 +1128,7 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
        exchange->h=&(params->h[num]);
 
        exchange->betan=pow(*(exchange->beta),*(exchange->n));
+       exchange->n_betan=*(exchange->n)*exchange->betan;
        exchange->c2=params->c[num]*params->c[num];
        exchange->d2=params->d[num]*params->d[num];
        exchange->c2d2=exchange->c2/exchange->d2;
@@ -1223,37 +1232,44 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        exchange->run3bp=1;
 
        /* reset 3bp sums */
-       exchange->3bp_sum1=0.0;
-       exchange->3bp_sum2=0.0;
+       exchange->sum1_3bp=0.0;
+       exchange->sum2_3bp=0.0;
+       v3_zero(&(exchange->db_ij));
 
        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) {
+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;
+       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,betan;
+       double help;
+       double n;
 
        params=moldyn->pot2b_params;
-       exchange=&(moldyn->exchange);
+       exchange=&(params->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;
+       betan=exchange->betan;
+       n=*(exchange->n);
+       dist_ij=&(exchange->dist_ij);
 
-       db_ij_scale1=(1+betan*3bp_sum1);
-       db_ij_scale2=(n*betan*3bp_sum2);
+       db_ij_scale1=(1+betan*exchange->sum1_3bp);
+       db_ij_scale2=(exchange->n_betan*exchange->sum2_3bp);
        help=pow(db_ij_scale1,-1.0/(2*n)-1);
        b_ij=chi*db_ij_scale1*help;
        db_ij_scale1=-chi/(2*n)*help;
@@ -1266,7 +1282,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        v3_add(&force,&temp,db_ij);
        v3_scale(&force,&force,f_c);
 
-       v3_scale(&temp,&dist_ij,f_a*b_ij*df_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);
@@ -1286,18 +1302,16 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        t_3dvec temp,force;
        double R,S,s_r;
        double d_ij,d_ij2,d_ik,d_jk;
-       double f_c,df_c,b_ij,f_a,df_a;
+       double f_c,df_c,f_a,df_a;
        double f_c_ik,df_c_ik,arg;
-       double scale;
-       double chi;
-       double n,c,d,h,beta,betan;
+       double n,c,d,h;
        double c2,d2,c2d2;
        double numer,denom;
        double theta,cos_theta,sin_theta;
        double d_theta,d_theta1,d_theta2;
-       double h_cos,h_cos2,d2_h_cos2;
-       double frac1,bracket1,bracket2,bracket2_n_1,bracket2_n;
-       double bracket3,bracket3_pow_1,bracket3_pow;
+       double h_cos,d2_h_cos2;
+       double frac,bracket,bracket_n_1,bracket_n;
+       double g;
        int num;
 
        params=moldyn->pot3b_params;
@@ -1369,13 +1383,11 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        if(bc) check_per_bound(moldyn,&dist_jk);
        d_jk=v3_norm(&dist_jk);
 
-       beta=*(exchange->beta);
-       betan=exchange->betan;
+       /* get exchange data */
        n=*(exchange->n);
        c=*(exchange->c);
        d=*(exchange->d);
        h=*(exchange->h);
-       chi=exchange->chi;
        c2=exchange->c2;
        d2=exchange->d2;
        c2d2=exchange->c2d2;
@@ -1383,6 +1395,7 @@ 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;
+       /* prefere law of cosines, dot product -> nan (often) */
        //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);
@@ -1393,61 +1406,43 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        d_theta2*=d_theta;
 
        h_cos=(h-cos_theta);
-       h_cos2=h_cos*h_cos;
-       d2_h_cos2=d2+h_cos2;
+       d2_h_cos2=d2+(h_cos*h_cos);
+
+       frac=c2/(d2_h_cos2);
+       g=1.0+c2d2-frac;
 
-       /* some usefull expressions */
-       frac1=c2/(d2_h_cos2);
-       bracket1=1+c2d2-frac1;
        if(f_c_ik==0.0) {
-               bracket2=0.0;
-               bracket2_n_1=0.0;
-               bracket2_n=0.0;
-               bracket3=1.0;
+               bracket=0.0;
+               bracket_n_1=0.0;
+               bracket_n=0.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;
+               bracket=f_c_ik*g;
+               bracket_n_1=pow(bracket,n-1.0);
+               bracket_n=bracket_n_1*bracket;
                printf("Foo -> 1: ");
        }
-       bracket3_pow_1=pow(bracket3,(-1.0/(2.0*n))-1.0);
-       bracket3_pow=bracket3_pow_1*bracket3;
-printf("%.15f %.15f %.15f\n",bracket2_n_1,bracket2_n);
+//printf("%.15f %.15f %.15f\n",bracket_n_1,bracket_n,bracket);
 
-       /* now go on with calc of b_ij and derivation of b_ij */
-       b_ij=chi*bracket3_pow;
+       /* calc of db_ij and the 2 sums */
+       exchange->sum1_3bp+=bracket_n;
+       exchange->sum2_3bp+=bracket_n_1;
 
        /* derivation of theta */
        v3_scale(&force,&dist_ij,d_theta1);
        v3_scale(&temp,&dist_ik,d_theta2);
        v3_add(&force,&force,&temp);
 
-       /* part 1 of derivation of b_ij */
-       v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac1);
-
-       /* part 2 of derivation of b_ij */
-       v3_scale(&temp,&dist_ik,df_c_ik*bracket1);
+       /* part 1 of db_ij */
+       v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac/d2_h_cos2);
 
-       /* sum up and scale ... */
-       v3_add(&temp,&temp,&force);
-       scale=bracket2_n_1*n*betan*(1+betan*bracket3_pow_1)*chi*(1.0/(2.0*n));
-       v3_scale(&temp,&temp,scale);
+       /* part 2 of db_ij */
+       v3_scale(&temp,&dist_ik,df_c_ik*g);
 
-       /* now construct an energy and a force out of that */
-       v3_scale(&temp,&temp,f_a);
-       v3_scale(&force,&dist_ij,df_a*b_ij);
+       /* sum up and add to db_ij */
        v3_add(&temp,&temp,&force);
-       v3_scale(&temp,&temp,f_c);
-       v3_scale(&force,&dist_ij,df_c*b_ij*f_a);
-       v3_add(&force,&force,&temp);
-
-       /* add forces */
-       v3_add(&(ai->f),&(ai->f),&force);
-       /* energy is 0.5 f_r f_c */
-       moldyn->energy+=(0.5*f_a*b_ij*f_c);
+       v3_add(&(exchange->db_ij),&(exchange->db_ij),&temp);
                                
        return 0;
 }

diff --git a/moldyn.h b/moldyn.h
index 41b711a..09471a1 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -194,12 +194,13 @@ typedef struct s_tersoff_exchange {
        double d2;
        double c2d2;
        double betan;
+       double n_betan;
 
        u8 run3bp;
 
        t_3dvec db_ij;
-       double 3bp_sum1;
-       double 3bp_sum2;
+       double sum1_3bp;
+       double sum2_3bp;
 } t_tersoff_exchange;
 
 /* tersoff multi (2!) potential parameters */
@@ -372,6 +373,7 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
 int tersoff_mult_complete_params(t_tersoff_mult_params *p);
 int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai);
 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
 int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
 
 #endif

diff --git a/sic.c b/sic.c
index 3408895..5f50bd8 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -41,6 +41,7 @@ int main(int argc,char **argv) {
        printf("[sic] selecting potential\n");
        set_potential1b(&md,tersoff_mult_1bp,&tp);
        set_potential2b(&md,tersoff_mult_2bp,&tp);
+       set_potential2b_post(&md,tersoff_mult_post_2bp,&tp);
        set_potential3b(&md,tersoff_mult_3bp,&tp);
        //set_potential2b(&md,lennard_jones,&lj);
 
@@ -125,7 +126,7 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        printf("[sic] adding schedule\n");
-       moldyn_add_schedule(&md,10,1.0e-15);
+       moldyn_add_schedule(&md,100,1.0e-15);
 
        /* activate logging */
        printf("[sic] activate logging\n");