removed tersoff 1bp function, added PDEBUG define, check per bound with boundary...

diff --git a/moldyn.c b/moldyn.c
index 7698c24..12cb49b 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -258,7 +258,7 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 
        switch(type) {
                case MOLDYN_POTENTIAL_TM:
-                       moldyn->func1b=tersoff_mult_1bp;
+                       //moldyn->func1b=tersoff_mult_1bp;
                        moldyn->func3b_j1=tersoff_mult_3bp_j1;
                        moldyn->func3b_k1=tersoff_mult_3bp_k1;
                        moldyn->func3b_j2=tersoff_mult_3bp_j2;
@@ -802,6 +802,8 @@ int del_atom(t_moldyn *moldyn,int tag) {
                        case DEFECT_STYPE_DB_Z:\
                                d_o.z=d_params->od;\
                                d_d.z=d_params->dd;\
+d_d.x=0.9;\
+d_d.y=0.9;\
                                break;\
                        case DEFECT_STYPE_DB_R:\
                                break;\
@@ -2071,11 +2073,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
                temperature_calc(moldyn);
                virial_sum(moldyn);
                pressure_calc(moldyn);
-               /*
+#ifdef PDEBUG  
                thermodynamic_pressure_calc(moldyn);
                printf("\n\nDEBUG: numeric pressure calc: %f\n\n",
                       moldyn->tp/BAR);
-               */
+#endif
 
                /* calculate fluctuations + averages */
                average_and_fluctuation_calc(moldyn);
@@ -2182,15 +2184,20 @@ int moldyn_integrate(t_moldyn *moldyn) {
                }
 
                /* display progress */
+#ifndef PDEBUG
                if(!(i%10)) {
+#endif
                        /* get current time */
                        gettimeofday(&t2,NULL);
 
 printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
        sched->count,i,moldyn->total_steps,
        moldyn->t,moldyn->t_avg,
+#ifndef PDEBUG
        moldyn->p/BAR,moldyn->p_avg/BAR,
-       //moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#else
+       moldyn->p/BAR,(moldyn->p-2.0*moldyn->ekin/(3.0*moldyn->volume))/BAR,
+#endif
        moldyn->volume,
        (int)(t2.tv_sec-t1.tv_sec));
 
@@ -2198,7 +2205,9 @@ printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
 
                        /* copy over time */
                        t1=t2;
+#ifndef PDEBUG
                }
+#endif
 
                /* increase absolute time */
                moldyn->time+=moldyn->tau;
@@ -2246,6 +2255,7 @@ int velocity_verlet(t_moldyn *moldyn) {
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
                v3_scale(&delta,&(atom[i].f),h*tau_square);
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
+               //check_per_bound_and_care_for_pbc(moldyn,&(atom[i]));
                check_per_bound(moldyn,&(atom[i].r));
 
                /* velocities [actually v(t+tau/2)] */
@@ -2261,7 +2271,11 @@ int velocity_verlet(t_moldyn *moldyn) {
 
        /* forces depending on chosen potential */
 #ifndef ALBE_FAST
-       potential_force_calc(moldyn);
+       // if albe, use albe force calc routine
+       //if(moldyn->func3b_j1==albe_mult_3bp_j1)
+       //      albe_potential_force_calc(moldyn);
+       //else
+               potential_force_calc(moldyn);
 #else
        albe_potential_force_calc(moldyn);
 #endif
@@ -2713,6 +2727,51 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
        return 0;
 }
         
+int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a) {
+       
+       double x,y,z;
+       t_3dvec *dim;
+
+       dim=&(moldyn->dim);
+
+       x=dim->x/2;
+       y=dim->y/2;
+       z=dim->z/2;
+
+       if(moldyn->status&MOLDYN_STAT_PBX) {
+               if(a->r.x>=x) {
+                       a->pbc[0]+=1;
+                       a->r.x-=dim->x;
+               }
+               else if(-a->r.x>x) {
+                       a->pbc[0]-=1;
+                       a->r.x+=dim->x;
+               }
+       }
+       if(moldyn->status&MOLDYN_STAT_PBY) {
+               if(a->r.y>=y) {
+                       a->pbc[1]+=1;
+                       a->r.y-=dim->y;
+               }
+               else if(-a->r.y>y) {
+                       a->pbc[1]-=1;
+                       a->r.y+=dim->y;
+               }
+       }
+       if(moldyn->status&MOLDYN_STAT_PBZ) {
+               if(a->r.z>=z) {
+                       a->pbc[2]+=1;
+                       a->r.z-=dim->z;
+               }
+               else if(-a->r.z>z) {
+                       a->pbc[2]-=1;
+                       a->r.z+=dim->z;
+               }
+       }
+
+       return 0;
+}
+        
 /*
  * debugging / critical check functions
  */
@@ -3048,6 +3107,7 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
        int i;
        t_atom *atom;
        t_3dvec dist;
+       t_3dvec final_r;
        double d2;
        int a_cnt;
        int b_cnt;
@@ -3061,8 +3121,12 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
 
        for(i=0;i<moldyn->count;i++) {
 
-               v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
-               check_per_bound(moldyn,&dist);
+               /* care for pb crossing */
+               final_r.x=atom[i].r.x+atom[i].pbc[0]*moldyn->dim.x;
+               final_r.y=atom[i].r.y+atom[i].pbc[1]*moldyn->dim.y;
+               final_r.z=atom[i].r.z+atom[i].pbc[2]*moldyn->dim.z;
+               /* calculate distance */
+               v3_sub(&dist,&final_r,&(atom[i].r_0));
                d2=v3_absolute_square(&dist);
 
                if(atom[i].brand) {

diff --git a/moldyn.h b/moldyn.h
index 538e0eb..6310416 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -45,6 +45,7 @@ typedef struct s_atom {
        u8 brand;               /* brand id */
        int tag;                /* atom unique id (number of atom) */
        u8 attr;                /* attributes */
+       int pbc[3];             /* pb crossing in x, y and z direction */
 } t_atom;
 
 #define ATOM_ATTR_FP   0x01    /* fixed position (bulk material) */
@@ -290,6 +291,11 @@ typedef struct s_defect_params {
 #define DEFECT_STYPE_DB_Z      3
 #define DEFECT_STYPE_DB_R      4
 
+typedef struct s_offset_params {
+       t_3dvec o;
+       u8 use;
+} t_offset_params;
+
 /*
  *
  *  defines
@@ -496,6 +502,7 @@ int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d);
 //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d)
 //     __attribute__((always_inline));
 int check_per_bound(t_moldyn *moldyn,t_3dvec *a);
+int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a);
 //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a)
 //     __attribute__((always_inline));
 

diff --git a/potentials/albe.c b/potentials/albe.c
index 08ce4fa..1ec3938 100644 (file)
--- a/potentials/albe.c
+++ b/potentials/albe.c
@@ -270,6 +270,11 @@ int albe_mult_3bp_k1(t_moldyn *moldyn,
                exchange->zeta_ij+=f_c_ik*g;
        }
 
+#ifdef DEBUG
+       if(ai==&(moldyn->atom[DATOM])) 
+               printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik);
+#endif
+
        /* store even more data for second k loop */
        exchange->g[kcount]=g;
        exchange->dg[kcount]=dg;
@@ -367,19 +372,17 @@ int albe_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        virial_calc(ai,&force,&(exchange->dist_ij));
 
 #ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
        if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
                printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
                printf("  adding %f %f %f\n",force.x,force.y,force.z);
-               if(ai==&(moldyn->atom[0]))
+               if(ai==&(moldyn->atom[DATOM]))
                        printf("  total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
-               if(aj==&(moldyn->atom[0]))
+               if(aj==&(moldyn->atom[DATOM]))
                        printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
                printf("  energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
                                                    f_c,b,f_a,f_r);
                printf("          %f %f %f\n",exchange->zeta_ij,.0,.0);
        }
-}
 #endif
 
        /* dzeta prefactor = - f_c f_a db, (* -0.5 due to force calc) */
@@ -477,7 +480,6 @@ int albe_mult_3bp_k2(t_moldyn *moldyn,
        v3_add(&(aj->f),&(aj->f),&force);
 
 #ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
        if(aj==&(moldyn->atom[DATOM])) {
                printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
                printf("  adding %f %f %f\n",force.x,force.y,force.z);
@@ -485,7 +487,6 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
                printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
                printf("    d ij ik = %f %f\n",d_ij,d_ik);
        }
-}
 #endif
 
        /* virial */
@@ -505,7 +506,6 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
        v3_add(&(ak->f),&(ak->f),&force);
 
 #ifdef DEBUG
-if(moldyn->time>DSTART&&moldyn->time<DEND) {
        if(ak==&(moldyn->atom[DATOM])) {
                printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
                printf("  adding %f %f %f\n",force.x,force.y,force.z);
@@ -513,7 +513,6 @@ if(moldyn->time>DSTART&&moldyn->time<DEND) {
                printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
                printf("    d ij ik = %f %f\n",d_ij,d_ik);
        }
-}
 #endif
 
        /* virial */

diff --git a/potentials/tersoff.c b/potentials/tersoff.c
index 5039ff7..4ca7141 100644 (file)
--- a/potentials/tersoff.c
+++ b/potentials/tersoff.c
@@ -85,20 +85,28 @@ int tersoff_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
                        p->h[1]=TM_H_C;
                        break;
                default:
-                       printf("[tersoff] WARNING: element1\n");
+                       printf("[tersoff] WARNING: element2\n");
                        return -1;
        }
 
        printf("[tersoff] parameter completion\n");
        p->S2[0]=p->S[0]*p->S[0];
        p->S2[1]=p->S[1]*p->S[1];
-       p->Smixed=sqrt(p->S[0]*p->S[1]);
-       p->S2mixed=p->Smixed*p->Smixed;
+       p->S2mixed=p->S[0]*p->S[1];
+       p->Smixed=sqrt(p->S2mixed);
        p->Rmixed=sqrt(p->R[0]*p->R[1]);
        p->Amixed=sqrt(p->A[0]*p->A[1]);
        p->Bmixed=sqrt(p->B[0]*p->B[1]);
        p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]);
        p->mu_m=0.5*(p->mu[0]+p->mu[1]);
+       p->betaini[0]=pow(p->beta[0],p->n[0]);
+       p->betaini[1]=pow(p->beta[1],p->n[1]);
+       p->ci2[0]=p->c[0]*p->c[0];
+       p->ci2[1]=p->c[1]*p->c[1];
+       p->di2[0]=p->d[0]*p->d[0];
+       p->di2[1]=p->d[1]*p->d[1];
+       p->ci2di2[0]=p->ci2[0]/p->di2[0];
+       p->ci2di2[1]=p->ci2[1]/p->di2[1];
 
        printf("[tersoff] mult parameter info:\n");
        printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
@@ -118,36 +126,6 @@ int tersoff_mult_set_params(t_moldyn *moldyn,int element1,int element2) {
        return 0;
 }
 
-/* tersoff 1 body part */
-int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
-
-       int brand;
-       t_tersoff_mult_params *params;
-       t_tersoff_exchange *exchange;
-       
-       brand=ai->brand;
-       params=moldyn->pot_params;
-       exchange=&(params->exchange);
-
-       /*
-        * simple: point constant parameters only depending on atom i to
-        *         their right values
-        */
-
-       exchange->beta_i=&(params->beta[brand]);
-       exchange->n_i=&(params->n[brand]);
-       exchange->c_i=&(params->c[brand]);
-       exchange->d_i=&(params->d[brand]);
-       exchange->h_i=&(params->h[brand]);
-
-       exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
-       exchange->ci2=params->c[brand]*params->c[brand];
-       exchange->di2=params->d[brand]*params->d[brand];
-       exchange->ci2di2=exchange->ci2/exchange->di2;
-
-       return 0;
-}
-       
 /* tersoff 2 body part */
 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
@@ -239,7 +217,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 #endif
 
        /* virial */
-       virial_calc(aj,&force,&dist_ij);
+       virial_calc(ai,&force,&dist_ij);
 
        /* energy 2bp contribution */
        energy=f_r*f_c;
@@ -364,10 +342,10 @@ int tersoff_mult_3bp_k1(t_moldyn *moldyn,
        cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
 
        /* g_ijk */
-       h_cos=*(exchange->h_i)-cos_theta;
-       d2_h_cos2=exchange->di2+(h_cos*h_cos);
-       frac=exchange->ci2/d2_h_cos2;
-       g=1.0+exchange->ci2di2-frac;
+       h_cos=params->h[brand]-cos_theta;
+       d2_h_cos2=params->di2[brand]+(h_cos*h_cos);
+       frac=params->ci2[brand]/d2_h_cos2;
+       g=1.0+params->ci2di2[brand]-frac;
        dg=-2.0*frac*h_cos/d2_h_cos2;
 
        /* zeta sum += f_c_ik * g_ijk */
@@ -385,11 +363,8 @@ int tersoff_mult_3bp_k1(t_moldyn *moldyn,
        }
 
 #ifdef DEBUG
-       if((ai==&(moldyn->atom[0]))|
-          (aj==&(moldyn->atom[864]))|
-          (ak==&(moldyn->atom[1003]))) {
-               printf(" -> %f %f %f\n",exchange->ci2di2,frac,h_cos);
-       }
+       if(ai==&(moldyn->atom[DATOM]))
+               printf("zeta_ij: %f %f %f %f\n",f_c_ik*g,f_c_ik,g,d_ik);
 #endif
 
        /* store even more data for second k loop */
@@ -425,8 +400,8 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        params=moldyn->pot_params;
        exchange=&(params->exchange);
 
-       brand=aj->brand;
-       if(brand==ai->brand) {
+       brand=ai->brand;
+       if(brand==aj->brand) {
                S=params->S[brand];
                R=params->R[brand];
                B=params->B[brand];
@@ -473,8 +448,8 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                db=0.0;
        }
        else {
-               ni=*(exchange->n_i);
-               tmp=exchange->betaini*pow(exchange->zeta_ij,ni-1.0);
+               ni=params->n[brand];
+               tmp=params->betaini[brand]*pow(exchange->zeta_ij,ni-1.0);
                b=(1.0+exchange->zeta_ij*tmp);
                db=chi*pow(b,-1.0/(2.0*ni)-1.0);
                b=db*b;
@@ -489,12 +464,12 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        v3_add(&(aj->f),&(aj->f),&force);
 
 #ifdef DEBUG
-       if((ai==&(moldyn->atom[0]))|(aj==&(moldyn->atom[0]))) {
+       if((ai==&(moldyn->atom[DATOM]))|(aj==&(moldyn->atom[DATOM]))) {
                printf("force 3bp (j2): [%d %d sum]\n",ai->tag,aj->tag);
                printf("adding %f %f %f\n",force.x,force.y,force.z);
-               if(ai==&(moldyn->atom[0]))
+               if(ai==&(moldyn->atom[DATOM]))
                        printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
-               if(aj==&(moldyn->atom[0]))
+               if(aj==&(moldyn->atom[DATOM]))
                        printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
                printf("energy: %f = %f %f %f %f\n",0.5*f_c*(b*f_a+f_r),
                                                    f_c,b,f_a,f_r);
@@ -503,7 +478,7 @@ int tersoff_mult_3bp_j2(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 #endif
 
        /* virial */
-       virial_calc(aj,&force,&(exchange->dist_ij));
+       virial_calc(ai,&force,&(exchange->dist_ij));
 
        /* dzeta prefactor = - 0.5 f_c f_a db */
        exchange->pre_dzeta=-0.5*f_a*f_c*db;
@@ -534,7 +509,7 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
        double pre_dzeta;
        double f_c_ik,df_c_ik;
        double dijdik_inv,fcdg,dfcg;
-       t_3dvec dcosdri,dcosdrj,dcosdrk;
+       t_3dvec dcosdrj,dcosdrk;
        t_3dvec force,tmp;
 
        params=moldyn->pot_params;
@@ -588,30 +563,11 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
        v3_scale(&dcosdrk,&dist_ij,dijdik_inv);
        v3_scale(&tmp,&dist_ik,-cos_theta/d_ik2);
        v3_add(&dcosdrk,&dcosdrk,&tmp);
-       v3_add(&dcosdri,&dcosdrj,&dcosdrk);
-       v3_scale(&dcosdri,&dcosdri,-1.0);
 
        /* f_c_ik * dg, df_c_ik * g */
        fcdg=f_c_ik*dg;
        dfcg=df_c_ik*g;
 
-       /* derivative wrt i */
-       v3_scale(&force,&dist_ik,dfcg);
-       v3_scale(&tmp,&dcosdri,fcdg);
-       v3_add(&force,&force,&tmp);
-       v3_scale(&force,&force,pre_dzeta);
-
-       /* force contribution */
-       v3_add(&(ai->f),&(ai->f),&force);
-       
-#ifdef DEBUG
-       if(ai==&(moldyn->atom[0])) {
-               printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
-               printf("adding %f %f %f\n",force.x,force.y,force.z);
-               printf("total i: %f %f %f\n",ai->f.x,ai->f.y,ai->f.z);
-       }
-#endif
-
        /* derivative wrt j */
        v3_scale(&force,&dcosdrj,fcdg*pre_dzeta);
 
@@ -619,17 +575,22 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
        v3_add(&(aj->f),&(aj->f),&force);
 
 #ifdef DEBUG
-       if(aj==&(moldyn->atom[0])) {
+       if(aj==&(moldyn->atom[DATOM])) {
                printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
-               printf("adding %f %f %f\n",force.x,force.y,force.z);
-               printf("total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+               printf("  adding %f %f %f\n",force.x,force.y,force.z);
+               printf("  total j: %f %f %f\n",aj->f.x,aj->f.y,aj->f.z);
+               printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
+               printf("    d ij ik = %f %f\n",d_ij,d_ik);
        }
 #endif
 
        /* virial */
-       v3_scale(&force,&force,-1.0);
        virial_calc(ai,&force,&dist_ij);
 
+       /* force contribution to atom i */
+       v3_scale(&force,&force,-1.0);
+       v3_add(&(ai->f),&(ai->f),&force);
+
        /* derivative wrt k */
        v3_scale(&force,&dist_ik,-1.0*dfcg); // dri rik = - drk rik
        v3_scale(&tmp,&dcosdrk,fcdg);
@@ -640,19 +601,24 @@ int tersoff_mult_3bp_k2(t_moldyn *moldyn,
        v3_add(&(ak->f),&(ak->f),&force);
 
 #ifdef DEBUG
-       if(ak==&(moldyn->atom[0])) {
+       if(ak==&(moldyn->atom[DATOM])) {
                printf("force 3bp (k2): [%d %d %d]\n",ai->tag,aj->tag,ak->tag);
-               printf("adding %f %f %f\n",force.x,force.y,force.z);
-               printf("total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
+               printf("  adding %f %f %f\n",force.x,force.y,force.z);
+               printf("  total k: %f %f %f\n",ak->f.x,ak->f.y,ak->f.z);
+               printf("  angle: %f\n",acos(cos_theta)*360.0/(2*M_PI));
+               printf("    d ij ik = %f %f\n",d_ij,d_ik);
        }
 #endif
 
        /* virial */
-       v3_scale(&force,&force,-1.0);
        virial_calc(ai,&force,&dist_ik);
+
+       /* force contribution to atom i */
+       v3_scale(&force,&force,-1.0);
+       v3_add(&(ai->f),&(ai->f),&force);
        
        /* increase k counter */
-       exchange->kcount++;     
+       exchange->kcount++;
 
        return 0;
 

diff --git a/potentials/tersoff.h b/potentials/tersoff.h
index a23503a..2e4f242 100644 (file)
--- a/potentials/tersoff.h
+++ b/potentials/tersoff.h
@@ -28,17 +28,6 @@ typedef struct s_tersoff_echange {
        double dg[TERSOFF_MAXN];
        double cos_theta[TERSOFF_MAXN];
 
-       double *beta_i;
-       double *n_i;
-       double *c_i;
-       double *d_i;
-       double *h_i;
-
-       double ci2;
-       double di2;
-       double ci2di2;
-       double betaini;
-
        double zeta_ij;
        double pre_dzeta;
 
@@ -70,12 +59,17 @@ typedef struct s_tersoff_mult_params {
        double d[2];
        double h[2];
 
+       double ci2[2];
+       double di2[2];
+       double ci2di2[2];
+       double betaini[2];
+
        t_tersoff_exchange exchange;    /* exchange between 2bp and 3bp calc */
 } t_tersoff_mult_params;
 
 /* function prototypes */
 int tersoff_mult_set_params(t_moldyn *moldyn,int element1,int element2);
-int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai);
+//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_3bp_j1(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
 int tersoff_mult_3bp_k1(t_moldyn *moldyn,