check_per_bound, comparing real numbers <- problem! think about it!

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 7ea5d50..44b53b6 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -99,6 +99,12 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
                moldyn->vis.dim.z=z;
        }
 
+       printf("[moldyn] dimensions in A:\n");
+       printf("  x: %f\n",moldyn->dim.x);
+       printf("  y: %f\n",moldyn->dim.y);
+       printf("  z: %f\n",moldyn->dim.z);
+       printf("  visualize simulation box: %s\n",visualize?"on":"off");
+
        return 0;
 }
 
@@ -230,10 +236,11 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 
        count=a*b*c;
 
+       /* how many atoms do we expect */
        if(type==FCC) count*=4;
-
        if(type==DIAMOND) count*=8;
 
+       /* allocate space for atoms */
        moldyn->atom=malloc(count*sizeof(t_atom));
        if(moldyn->atom==NULL) {
                perror("malloc (atoms)");
@@ -256,9 +263,10 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 
        /* debug */
        if(ret!=count) {
-               printf("ok, there is something wrong ...\n");
-               printf("calculated -> %d atoms\n",count);
-               printf("created -> %d atoms\n",ret);
+               printf("[moldyn] creating lattice failed\n");
+               printf("  amount of atoms\n");
+               printf("  - expected: %d\n",count);
+               printf("  - created: %d\n",ret);
                return -1;
        }
 
@@ -274,7 +282,6 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
                check_per_bound(moldyn,&(moldyn->atom[count].r));
        }
 
-
        return ret;
 }
 
@@ -802,6 +809,7 @@ int velocity_verlet(t_moldyn *moldyn) {
                v3_add(&(atom[i].v),&(atom[i].v),&delta);
        }
 
+moldyn_bc_check(moldyn);
        /* neighbour list update */
        link_cell_update(moldyn);
 
@@ -885,11 +893,12 @@ int potential_force_calc(t_moldyn *moldyn) {
                                        continue;
 
                                if((jtom->attr&ATOM_ATTR_2BP)&
-                                  (itom[i].attr&ATOM_ATTR_2BP))
+                                  (itom[i].attr&ATOM_ATTR_2BP)) {
                                        moldyn->func2b(moldyn,
                                                       &(itom[i]),
                                                       jtom,
                                                       bc_ij);
+                               }
 
                                /* 3 body potential/force */
 
@@ -963,9 +972,9 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
 
        dim=&(moldyn->dim);
 
-       x=0.5*dim->x;
-       y=0.5*dim->y;
-       z=0.5*dim->z;
+       x=dim->x/2;
+       y=dim->y/2;
+       z=dim->z/2;
 
        if(moldyn->status&MOLDYN_STAT_PBX) {
                if(a->x>=x) a->x-=dim->x;
@@ -979,6 +988,7 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
                if(a->z>=z) a->z-=dim->z;
                else if(-a->z>z) a->z+=dim->z;
        }
+printf("%f %f %f\n",a->x,x,a->x/x);
 
        return 0;
 }
@@ -1212,7 +1222,7 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 
        /* f_a, df_a calc (again, same for ij and ji) | save for later use! */
        exchange->f_a=-B*exp(-mu*d_ij);
-       exchange->df_a=-mu*exchange->f_a/d_ij;
+       exchange->df_a=mu*exchange->f_a/d_ij;
 
        /* f_c, df_c calc (again, same for ij and ji) */
        if(d_ij<R) {
@@ -1354,7 +1364,7 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        v3_scale(&force,&force,-0.5);
 
        /* add force */
-       v3_sub(&(ai->f),&(ai->f),&force);
+       v3_add(&(ai->f),&(ai->f),&force);
 
        return 0;
 }
@@ -1634,13 +1644,13 @@ int moldyn_bc_check(t_moldyn *moldyn) {
        for(i=0;i<moldyn->count;i++) {
                if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2)
                        printf("FATAL: atom %d: x: %.20f (%.20f)\n",
-                              i,atom[i].r.x*1e10,dim->x/2*1e10);
+                              i,atom[i].r.x,dim->x/2);
                if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
                        printf("FATAL: atom %d: y: %.20f (%.20f)\n",
-                              i,atom[i].r.y*1e10,dim->y/2*1e10);
+                              i,atom[i].r.y,dim->y/2);
                if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
                        printf("FATAL: atom %d: z: %.20f (%.20f)\n",
-                              i,atom[i].r.z*1e10,dim->z/2*1e10);
+                              i,atom[i].r.z,dim->z/2);
        }
 
        return 0;