into the tersoff potential

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index d52fcd3..ed831ae 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) {