fixed 2bp, todo: adjust create lattice, then go for 3bp

diff --git a/moldyn.c b/moldyn.c
index 7ea5d5080ab1ab185fa46b2137894239e9233902..aaab0deb5270c0157ebcbc4a07e58ddd12bf0943 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;
 }
 
@@ -802,6 +808,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 +892,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 */
 
@@ -1212,7 +1220,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 +1362,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;
 }

diff --git a/sic.c b/sic.c
index 1e29a85d9647341d316def20e9f1a097046d41f3..5edb699ef64d525c6b765ff64d0d9088f2090331 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -107,20 +107,24 @@ int main(int argc,char **argv) {
        set_pbc(&md,TRUE,TRUE,TRUE);
 
        /* create the lattice / place atoms */
-       //printf("[sic] creating atoms\n");
-       //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
-       //               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
-       //               0,5,5,5);
+       printf("[sic] creating atoms\n");
+       create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                      0,5,5,5);
 
        /* testing configuration */
-       r.x=2.7/2;      v.x=0;
-       r.y=0;          v.y=0;
-       r.z=0;          v.z=0;
-       add_atom(&md,SI,M_SI,0,ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,&r,&v);
-       r.x=-2.7/2;     v.x=0;
-       r.y=0;          v.y=0;
-       r.z=0;          v.z=0;
-       add_atom(&md,SI,M_SI,0,ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,&r,&v);
+       //r.x=2.45/2;   v.x=0;
+       //r.y=0;                v.y=0;
+       //r.z=0;                v.z=0;
+       //add_atom(&md,SI,M_SI,0,
+       //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+       //           &r,&v);
+       //r.x=-2.45/2;  v.x=0;
+       //r.y=0;                v.y=0;
+       //r.z=0;                v.z=0;
+       //add_atom(&md,SI,M_SI,0,
+       //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+       //           &r,&v);
 
        /* setting a nearest neighbour distance for the moldyn checks */
        set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */