into the tersoff potential

diff --git a/moldyn.c b/moldyn.c
index d52fcd3868c2481e70cbc804b5a9cdf92542708d..ed831aed93ec4861c71a383033af863ef1378074 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -231,13 +231,15 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        printf("[moldyn] created lattice with %d atoms\n",count);
 
        while(count) {
-               moldyn->atom[count-1].element=element;
-               moldyn->atom[count-1].mass=mass;
-               moldyn->atom[count-1].attr=attr;
-               moldyn->atom[count-1].bnum=bnum;
                count-=1;
+               moldyn->atom[count].element=element;
+               moldyn->atom[count].mass=mass;
+               moldyn->atom[count].attr=attr;
+               moldyn->atom[count].bnum=bnum;
+               check_per_bound(moldyn,&(moldyn->atom[count].r));
        }
 
+
        return ret;
 }
 
@@ -529,7 +531,7 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
                }
        }
 
-       lc->dnlc=count2;
+       lc->dnlc=count1;
        lc->countn=27;
 
        return count2;
@@ -618,7 +620,6 @@ int moldyn_integrate(t_moldyn *moldyn) {
        /* sqaure of some variables */
        moldyn->tau_square=moldyn->tau*moldyn->tau;
        moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
-
        /* calculate initial forces */
        potential_force_calc(moldyn);
 
@@ -774,7 +775,7 @@ int potential_force_calc(t_moldyn *moldyn) {
        moldyn->energy=0.0;
 
        for(i=0;i<count;i++) {
-       
+
                /* reset force */
                v3_zero(&(atom[i].f));
 
@@ -823,11 +824,13 @@ int potential_force_calc(t_moldyn *moldyn) {
                                           !(btom->attr&ATOM_ATTR_3BP))
                                                continue;
 
+printf("DEBUG: problem exists here ...\n");
                                        link_cell_neighbour_index(moldyn,
                                           (btom->r.x+moldyn->dim.x/2)/lc->x,
                                           (btom->r.y+moldyn->dim.y/2)/lc->y,
                                           (btom->r.z+moldyn->dim.z/2)/lc->z,
                                           neighbourk);
+printf("DEBUG: as you won't see that!\n");
 
                                        for(k=0;k<lc->countn;k++) {
 
@@ -971,6 +974,20 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
  * tersoff potential & force for 2 sorts of atoms
  */
 
+/* create mixed terms from parameters and set them */
+int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
+
+       printf("[moldyn] tersoff parameter completion\n");
+       p->Smixed=sqrt(p->S[0]*p->S[1]);
+       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]);
+
+       return 0;
+}
+
 /* tersoff 1 body part */
 int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
 

diff --git a/moldyn.h b/moldyn.h
index 63dd14f0d45af67c6b8f9782079626f7ca3cfa63..2d4b5634b60c07fc942d1ee9e0eeb9424588c680 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -254,6 +254,7 @@ typedef struct s_tersoff_mult_params {
 
 #define K_BOLTZMANN            1.3807e-27                      /* Nm/K */
 #define AMU                    1.660540e-27                    /* kg */
+#define EV                     1.60217733e-19                  /* Nm */
 
 #define FCC                    0x01
 #define DIAMOND                        0x02
@@ -267,6 +268,32 @@ typedef struct s_tersoff_mult_params {
 #define LJ_SIGMA_SI            ((0.25*sqrt(3.0)*LC_SI)/1.122462)       /* m */
 #define LJ_EPSILON_SI          (2.1678*1.60e-19)                       /* Nm */
 
+#define TM_R_SI                        2.7e-10                         /* m */
+#define TM_S_SI                        3.0e-10                         /* m */
+#define TM_A_SI                        (1830.8*EV)                     /* Nm */
+#define TM_B_SI                        (471.18*EV)                     /* Nm */
+#define TM_LAMBDA_SI           2.4799e10                       /* 1/m */
+#define TM_MU_SI               1.7322e10                       /* 1/m */
+#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_H_SI                        (-0.59825)
+
+#define TM_R_C                 1.8e-10                         /* m */
+#define TM_S_C                 2.1e-10                         /* m */
+#define TM_A_C                 (1393.6*EV)                     /* Nm */
+#define TM_B_C                 (346.7*EV)                      /* Nm */
+#define TM_LAMBDA_C            3.4879e10                       /* 1/m */
+#define TM_MU_C                        2.2119e10                       /* 1/m */
+#define TM_BETA_C              1.5724e-7
+#define TM_N_C                 0.72751
+#define TM_C_C                 3.8049e4
+#define TM_D_C                 4.384
+#define TM_H_C                 (-0.57058)
+
+#define TM_CHI_SIC             0.9776
+
 
 /*
  *
@@ -325,6 +352,7 @@ int potential_force_calc(t_moldyn *moldyn);
 int check_per_bound(t_moldyn *moldyn,t_3dvec *a);
 int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
 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_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);

diff --git a/sic.c b/sic.c
index 7c40a7e6c8f603f8ff375b08b40c30e8f37ce9ea..2cd32eb04853a4c14518561999f9d6a429d9e5a1 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -15,7 +15,7 @@
 #include "posic.h"
 
 int main(int argc,char **argv) {
-
+printf("%d\n",sizeof(t_atom));
        /* main moldyn structure */
        t_moldyn md;
 
@@ -24,12 +24,6 @@ int main(int argc,char **argv) {
        t_ho_params ho;
        t_tersoff_mult_params tp;
 
-       /* misc variables, mainly to initialize stuff */
-       t_3dvec r,v;
-
-       /* temperature */
-       double t;
-
        /* initialize moldyn */
        printf("[sic] moldyn init\n");
        moldyn_init(&md,argc,argv);
@@ -40,10 +34,10 @@ int main(int argc,char **argv) {
 
        /* choose potential */
        printf("[sic] selecting potential\n");
-       //set_potential1b(&md,tersoff_mult_1bp,&tp);
-       //set_potential2b(&md,tersoff_mult_2bp,&tp);
-       //set_potential3b(&md,tersoff_mult_3bp,&tp);
-       set_potential2b(&md,lennard_jones,&lj);
+       set_potential1b(&md,tersoff_mult_1bp,&tp);
+       set_potential2b(&md,tersoff_mult_2bp,&tp);
+       set_potential3b(&md,tersoff_mult_3bp,&tp);
+       //set_potential2b(&md,lennard_jones,&lj);
 
        /*
         * potential parameters
@@ -59,9 +53,39 @@ int main(int argc,char **argv) {
        ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
        ho.spring_constant=1;
 
+       /*
+         * tersoff mult potential parameters for SiC
+        */
+       tp.S[0]=TM_S_SI;
+       tp.R[0]=TM_R_SI;
+       tp.A[0]=TM_A_SI;
+       tp.B[0]=TM_B_SI;
+       tp.lambda[0]=TM_LAMBDA_SI;
+       tp.mu[0]=TM_MU_SI;
+       tp.beta[0]=TM_BETA_SI;
+       tp.n[0]=TM_N_SI;
+       tp.c[0]=TM_C_SI;
+       tp.d[0]=TM_D_SI;
+       tp.h[0]=TM_H_SI;
+
+       tp.S[1]=TM_S_C;
+       tp.R[1]=TM_R_C;
+       tp.A[1]=TM_A_C;
+       tp.B[1]=TM_B_C;
+       tp.lambda[1]=TM_LAMBDA_C;
+       tp.mu[1]=TM_MU_C;
+       tp.beta[1]=TM_BETA_C;
+       tp.n[1]=TM_N_C;
+       tp.c[1]=TM_C_C;
+       tp.d[1]=TM_D_C;
+
+       tp.chi=TM_CHI_SIC;
+
+       tersoff_mult_complete_params(&tp);
+
        /* cutoff radius */
        printf("[sic] setting cutoff radius\n");
-       set_cutoff(&md,1*LC_SI);
+       set_cutoff(&md,TM_S_SI);
 
        /* set (initial) dimensions of simulation volume */
        printf("[sic] setting dimensions\n");
@@ -73,21 +97,16 @@ int main(int argc,char **argv) {
 
        /* create the lattice / place atoms */
        printf("[sic] creating atoms\n");
-       //memset(&v,0,sizeof(t_3dvec));
-       //r.y=0;
-       //r.z=0;
-       //r.x=0.23*sqrt(3.0)*LC_SI/2.0;
-       //add_atom(&md,SI,M_SI,0,ATOM_ATTR_2BP,&r,&v);
-       //r.x=-r.x;
-       //add_atom(&md,SI,M_SI,0,ATOM_ATTR_2BP,&r,&v);
-       create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,ATOM_ATTR_2BP,0,4,4,4);
+       create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP,
+                      0,4,4,4);
 
        /* setting a nearest neighbour distance for the moldyn checks */
        set_nn_dist(&md,sqrt(3.0)*LC_SI/4.0); /* diamond ! */
 
        /* set temperature */
        printf("[sic] setting temperature\n");
-       set_temperature(&md,273.0);
+       set_temperature(&md,273.0+450.0);
        
        /* initial thermal fluctuations of particles */
        printf("[sic] thermal init\n");