tersoff MAYBE fixed now!

diff --git a/moldyn.c b/moldyn.c
index dae975c..ca6f84a 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -162,28 +162,28 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
        return 0;
 }
 
-int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params) {
+int set_potential1b(t_moldyn *moldyn,pf_func1b func) {
 
        moldyn->func1b=func;
 
        return 0;
 }
 
-int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) {
+int set_potential2b(t_moldyn *moldyn,pf_func2b func) {
 
        moldyn->func2b=func;
 
        return 0;
 }
 
-int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) {
+int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func) {
 
        moldyn->func2b_post=func;
 
        return 0;
 }
 
-int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) {
+int set_potential3b(t_moldyn *moldyn,pf_func3b func) {
 
        moldyn->func3b=func;
 

diff --git a/moldyn.h b/moldyn.h
index 1300998..b140372 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -230,9 +230,9 @@ typedef struct s_moldyn {
 #define LC_SI                  (0.543105e-9*METER)             /* A */
 #define M_SI                   28.08553                        /* amu */
 
-//#define LJ_SIGMA_SI          ((0.25*sqrt(3.0)*LC_SI)/1.122462)       /* A */
+#define LJ_SIGMA_SI            ((0.25*sqrt(3.0)*LC_SI)/1.122462)       /* A */
 //#define LJ_SIGMA_SI          (LC_SI/1.122462)                        /* A */
-#define LJ_SIGMA_SI            (0.5*sqrt(2.0)*LC_SI/1.122462)                  /* A */
+//#define LJ_SIGMA_SI          (0.5*sqrt(2.0)*LC_SI/1.122462)                  /* A */
 #define LJ_EPSILON_SI          (2.1678*EV)                             /* NA */
 
 #define TM_R_SI                        (2.7e-10*METER)                 /* A */
@@ -292,10 +292,11 @@ int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc);
 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize);
 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 set_potential1b(t_moldyn *moldyn,pf_func1b func);
+int set_potential2b(t_moldyn *moldyn,pf_func2b func);
+int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func);
+int set_potential3b(t_moldyn *moldyn,pf_func3b func);
+int set_potential_params(t_moldyn *moldyn,void *params);
 
 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir);
 int moldyn_set_report(t_moldyn *moldyn,char *author,char *title);

diff --git a/potentials/lennard_jones.c b/potentials/lennard_jones.c
index 1738ea5..1277c46 100644 (file)
--- a/potentials/lennard_jones.c
+++ b/potentials/lennard_jones.c
@@ -35,10 +35,11 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        v3_sub(&distance,&(aj->r),&(ai->r));
        if(bc) check_per_bound(moldyn,&distance);
-       d=v3_absolute_square(&distance);        /* 1/r^2 */
+       d=v3_absolute_square(&distance);        /* r^2 */
        if(d<=moldyn->cutoff_square) {
                d=1.0/d;                        /* 1/r^2 */
                h2=d*d;                         /* 1/r^4 */
+               h2*=d;                          /* 1/r^6 */
                h1=h2*h2;                       /* 1/r^12 */
                moldyn->energy+=(eps*(sig12*h1-sig6*h2)-params->uc);
                h2*=d;                          /* 1/r^8 */

diff --git a/potentials/tersoff.c b/potentials/tersoff.c
index 972075d..01364ea 100644 (file)
--- a/potentials/tersoff.c
+++ b/potentials/tersoff.c
@@ -95,6 +95,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        double s_r;
        double arg;
 
+       /* use newtons third law */
+       //if(ai<aj) return 0;
+
        params=moldyn->pot_params;
        brand=aj->brand;
        exchange=&(params->exchange);
@@ -119,27 +122,11 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
         *
         */
 
-       /* constants */
-       if(brand==ai->brand) {
-               S=params->S[brand];
+       /* determine cutoff square */
+       if(brand==ai->brand)
                S2=params->S2[brand];
-               R=params->R[brand];
-               A=params->A[brand];
-               B=params->B[brand];
-               lambda=params->lambda[brand];
-               mu=params->mu[brand];
-               exchange->chi=1.0;
-       }
-       else {
-               S=params->Smixed;
+       else
                S2=params->S2mixed;
-               R=params->Rmixed;
-               A=params->Amixed;
-               B=params->Bmixed;
-               lambda=params->lambda_m;
-               mu=params->mu_m;
-               params->exchange.chi=params->chi;
-       }
 
        /* dist_ij, d_ij */
        v3_sub(&dist_ij,&(aj->r),&(ai->r));
@@ -166,12 +153,26 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        exchange->d_j=&(params->d[brand]);
        exchange->h_j=&(params->h[brand]);
        if(brand==ai->brand) {
+               S=params->S[brand];
+               R=params->R[brand];
+               A=params->A[brand];
+               B=params->B[brand];
+               lambda=params->lambda[brand];
+               mu=params->mu[brand];
+               exchange->chi=1.0;
                exchange->betajnj=exchange->betaini;
                exchange->cj2=exchange->ci2;
                exchange->dj2=exchange->di2;
                exchange->cj2dj2=exchange->ci2di2;
        }
        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;
                exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
                exchange->cj2=params->c[brand]*params->c[brand];
                exchange->dj2=params->d[brand]*params->d[brand];

diff --git a/sic.c b/sic.c
index e148709..167735f 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -51,6 +51,8 @@ int main(int argc,char **argv) {
 
        /* testing location & velocity vector */
        t_3dvec r,v;
+       memset(&r,0,sizeof(t_3dvec));
+       memset(&v,0,sizeof(t_3dvec));
 
        /* initialize moldyn */
        moldyn_init(&md,argc,argv);
@@ -59,17 +61,20 @@ int main(int argc,char **argv) {
        set_int_alg(&md,MOLDYN_INTEGRATE_VERLET);
 
        /* choose potential */
-       //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);
-       //set_potential2b(&md,harmonic_oscillator,&ho);
+       set_potential1b(&md,tersoff_mult_1bp);
+       set_potential2b(&md,tersoff_mult_2bp);
+       set_potential2b_post(&md,tersoff_mult_post_2bp);
+       set_potential3b(&md,tersoff_mult_3bp);
+       //set_potential2b(&md,lennard_jones);
+       //set_potential2b(&md,harmonic_oscillator);
+       set_potential_params(&md,&tp);
+       //set_potential_params(&md,&lj);
+       //set_potential_params(&md,&ho);
 
        /* cutoff radius */
-       //set_cutoff(&md,TM_S_SI);
-       //set_cutoff(&md,2*LC_SI*0.5*sqrt(1.5));
-       set_cutoff(&md,2.0*LC_SI);
+       set_cutoff(&md,TM_S_SI);
+       //set_cutoff(&md,LC_SI*sqrt(3.0));
+       //set_cutoff(&md,2.0*LC_SI);
 
        /*
         * potential parameters
@@ -83,8 +88,8 @@ int main(int argc,char **argv) {
        lj.uc=lj.epsilon4*(lj.sigma12/pow(md.cutoff,12.0)-lj.sigma6/pow(md.cutoff,6));
 
        /* harmonic oscillator */
-       //ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
-       ho.equilibrium_distance=LC_SI;
+       ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
+       //ho.equilibrium_distance=LC_SI;
        ho.spring_constant=LJ_EPSILON_SI;
 
        /*
@@ -126,16 +131,16 @@ int main(int argc,char **argv) {
 
        /* create the lattice / place atoms */
        //create_lattice(&md,CUBIC,LC_SI,SI,M_SI,
-       create_lattice(&md,FCC,LC_SI,SI,M_SI,
-       //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,
+       //create_lattice(&md,FCC,LC_SI,SI,M_SI,
+       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,6,6,6);
        moldyn_bc_check(&md);
 
        /* testing configuration */
-       //r.x=0.28*sqrt(3)*LC_SI/2;     v.x=0;
-       //r.x=1.75*LC_SI;       v.x=-0.01;
+       //r.x=0.27*sqrt(3.0)*LC_SI/2.0; v.x=0;
+       //r.x=(TM_S_SI+TM_R_SI)/4.0;    v.x=0;
        //r.y=0;                v.y=0;
        //r.z=0;                v.z=0;
        //add_atom(&md,SI,M_SI,0,
@@ -149,6 +154,13 @@ int main(int argc,char **argv) {
        //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
        //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
        //           &r,&v);
+       //r.x=0;                v.x=0;
+       //r.y=0;        v.y=0;
+       //r.z=sin(M_PI*60.0/180.0)*(TM_S_SI+TM_R_SI)/4.0;       v.z=0;
+       //add_atom(&md,SI,M_SI,0,
+       //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+       //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
+       //           &r,&v);
 
        /* set temperature & pressure */
        set_temperature(&md,atof(argv[2])+273.0);
@@ -164,7 +176,7 @@ int main(int argc,char **argv) {
        thermal_init(&md,TRUE);
 
        /* create the simulation schedule */
-       moldyn_add_schedule(&md,101,1.0);
+       moldyn_add_schedule(&md,10001,1.0);
        //moldyn_add_schedule(&md,501,1.0);
        //moldyn_add_schedule(&md,501,1.0);
 
@@ -175,7 +187,7 @@ int main(int argc,char **argv) {
        moldyn_set_log_dir(&md,argv[1]);
        moldyn_set_report(&md,"Frank Zirkelbach","Test 1");
        moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
-       moldyn_set_log(&md,VISUAL_STEP,1);
+       moldyn_set_log(&md,VISUAL_STEP,100);
        moldyn_set_log(&md,CREATE_REPORT,0);
 
        /*