need sleep, continue debugging soon!

diff --git a/clean b/clean
index fc5f86d..ff74c9b 100755 (executable)
--- a/clean
+++ b/clean
@@ -1 +1 @@
-rm -f saves/* video/*.ppm video/*.png
+rm -f saves/* video/*.ppm video/*.png > /dev/null 2>&1

diff --git a/moldyn.c b/moldyn.c
index 6af2dff..3608658 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -576,9 +576,8 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
        }
 
        lc->dnlc=count1;
-       lc->countn=27;
 
-       return count2;
+       return count1;
 }
 
 int link_cell_shutdown(t_moldyn *moldyn) {
@@ -664,10 +663,11 @@ 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);
 
-       /* do some checks before we actually start calculating bullshit */
+       /* some stupid checks before we actually start calculating bullshit */
        if(moldyn->cutoff>0.5*moldyn->dim.x)
                printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
        if(moldyn->cutoff>0.5*moldyn->dim.y)
@@ -680,6 +680,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
        /* zero absolute time */
        moldyn->time=0.0;
+
+       /* debugging, ignre */
+       moldyn->debug=0;
+
+       /* executing the schedule */
        for(sched=0;sched<moldyn->schedule.content_count;sched++) {
 
                /* setting amount of runs and finite time step size */
@@ -736,7 +741,8 @@ int moldyn_integrate(t_moldyn *moldyn) {
                        if(!(i%v)) {
                                visual_atoms(&(moldyn->vis),moldyn->time,
                                             moldyn->atom,moldyn->count);
-                               printf("\rsched: %d, steps: %d",sched,i);
+                               printf("\rsched: %d, steps: %d, theta: %d",
+                                      sched,i,moldyn->debug);
                                fflush(stdout);
                        }
                }
@@ -747,6 +753,9 @@ int moldyn_integrate(t_moldyn *moldyn) {
                if(schedule->hook)
                        schedule->hook(moldyn,schedule->hook_params);
 
+               /* get a new info line */
+               printf("\n");
+
        }
 
        return 0;
@@ -812,8 +821,8 @@ int potential_force_calc(t_moldyn *moldyn) {
        t_list neighbour_i2[27];
        //t_list neighbour_j[27];
        t_list *this,*that;
-       u8 bc_ij,bc_ijk;
-       int countn,dnlc;
+       u8 bc_ij,bc_ik;
+       int dnlc;
 
        count=moldyn->count;
        itom=moldyn->atom;
@@ -822,6 +831,7 @@ int potential_force_calc(t_moldyn *moldyn) {
        /* reset energy */
        moldyn->energy=0.0;
 
+       /* get energy and force of every atom */
        for(i=0;i<count;i++) {
 
                /* reset force */
@@ -840,10 +850,9 @@ int potential_force_calc(t_moldyn *moldyn) {
                                (itom[i].r.z+moldyn->dim.z/2)/lc->z,
                                neighbour_i);
 
-                       countn=lc->countn;
                        dnlc=lc->dnlc;
 
-                       for(j=0;j<countn;j++) {
+                       for(j=0;j<27;j++) {
 
                                this=&(neighbour_i[j]);
                                list_reset(this);
@@ -872,51 +881,12 @@ int potential_force_calc(t_moldyn *moldyn) {
                                           !(jtom->attr&ATOM_ATTR_3BP))
                                                continue;
 
-                       /* neighbourhood of atom j is not needed! */
-
-                       //              link_cell_neighbour_index(moldyn,
-                       //                 (jtom->r.x+moldyn->dim.x/2)/lc->x,
-                       //                 (jtom->r.y+moldyn->dim.y/2)/lc->y,
-                       //                 (jtom->r.z+moldyn->dim.z/2)/lc->z,
-                       //                 neighbour_j);
-
-//                                     /* neighbours of j */
-//                                     for(k=0;k<lc->countn;k++) {
-//
-//                                             that=&(neighbour_j[k]);
-//                                             list_reset(that);
-//                                     
-//                                             if(that->start==NULL)
-//                                                     continue;
-//
-//                                             bc_ijk=(k<lc->dnlc)?0:1;
-//
-//                                             do {
-//
-//                     ktom=that->current->data;
-//
-//                     if(!(ktom->attr&ATOM_ATTR_3BP))
-//                             continue;
-//
-//                     if(ktom==jtom)
-//                             continue;
-//
-//                     if(ktom==&(itom[i]))
-//                             continue;
-//
-//                     moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
-//
-/*                                             } while(list_next(that)!=\ */
-//                                                     L_NO_NEXT_ELEMENT);
-//
-//                                     }
-                       
                                        /* copy the neighbour lists */
                                        memcpy(neighbour_i2,neighbour_i,
                                               27*sizeof(t_list));
 
                                        /* get neighbours of i */
-                                       for(k=0;k<countn;k++) {
+                                       for(k=0;k<27;k++) {
 
                                                that=&(neighbour_i2[k]);
                                                list_reset(that);
@@ -924,7 +894,7 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                if(that->start==NULL)
                                                        continue;
 
-                                               bc_ijk=(k<dnlc)?0:1;
+                                               bc_ik=(k<dnlc)?0:1;
 
                                                do {
 
@@ -939,7 +909,7 @@ int potential_force_calc(t_moldyn *moldyn) {
                        if(ktom==&(itom[i]))
                                continue;
 
-                       moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ijk);
+                       moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ik|bc_ij);
 
                                                } while(list_next(that)!=\
                                                        L_NO_NEXT_ELEMENT);
@@ -950,11 +920,11 @@ int potential_force_calc(t_moldyn *moldyn) {
                
                                /* 2bp post function */
                                if(moldyn->func2b_post) {
-printf("before 2bp post: %.15f %.15f %.15f\n",itom[i].f.x,itom[i].r.x,itom[i].v.x);
+printf("pre(%d): %.15f %.15f %.15f\n",i,itom[i].f.x,itom[i].r.x,itom[i].v.x);
                                        moldyn->func2b_post(moldyn,
                                                            &(itom[i]),
                                                            jtom,bc_ij);
-printf("after 2bp post: %.15f %.15f %.15f\n",itom[i].f.x,itom[i].r.x,itom[i].v.x);
+printf("post(%d): %.15f %.15f %.15f\n",i,itom[i].f.x,itom[i].r.x,itom[i].v.x);
                                }
 
                        }
@@ -1252,7 +1222,7 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        params=moldyn->pot2b_params;
        exchange=&(params->exchange);
 
-       /* we do not run if f_c_ij was dtected to be 0! */
+       /* we do not run if f_c_ij was detected to be 0! */
        if(!(exchange->run2bp_post))
                return 0;
 
@@ -1271,6 +1241,7 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        help=pow(db_ij_scale1,-1.0/(2*n)-1);
        b_ij=chi*db_ij_scale1*help;
        db_ij_scale1=-chi/(2*n)*help;
+printf("debug: %.20f %.20f %.20f\n",db_ij->x,exchange->sum1_3bp,exchange->sum2_3bp);
 
        /* db_ij part */
        v3_scale(db_ij,db_ij,(db_ij_scale1*db_ij_scale2));
@@ -1395,11 +1366,26 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        d2=exchange->d2;
        c2d2=exchange->c2d2;
 
+       /* cosine of theta by scalaproduct,       *
+        * derivation of theta by law of cosines! */
        numer=d_ij2+d_ik*d_ik-d_jk*d_jk;
        denom=2*d_ij*d_ik;
        cos_theta=numer/denom;
-       /* prefere law of cosines, dot product -> nan (often) */
-       //cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
+       cos_theta=v3_scalar_product(&dist_ij,&dist_ik)/(d_ij*d_ik);
+printf("cos theta: %.25f\n",cos_theta);
+
+       /* hack - cos theta machine accuracy problems! */
+       if(cos_theta>1.0||cos_theta<-1.0) {
+               moldyn->debug++;
+               if(fabs(cos_theta)>1.0+ACCEPTABLE_ERROR)
+                       printf("[moldyn] WARNING: cos theta failure!\n");
+               if(cos_theta<0)
+                       cos_theta=-1.0;
+               else
+                       cos_theta=1.0;
+               printf("THETA CORRECTION\n");
+       }
+
        sin_theta=sqrt(1.0-(cos_theta*cos_theta));
        theta=acos(cos_theta);
        d_theta=(-1.0/sqrt(1.0-cos_theta*cos_theta))/(denom*denom);
@@ -1407,6 +1393,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        d_theta2=2*denom-numer*2*d_ij/d_ik;
        d_theta1*=d_theta;
        d_theta2*=d_theta;
+printf("FOO %.15f %.15f\n",sin_theta,cos_theta);
 
        h_cos=(h-cos_theta);
        d2_h_cos2=d2+(h_cos*h_cos);
@@ -1434,6 +1421,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        v3_scale(&temp,&dist_ik,d_theta2);
        v3_add(&force,&force,&temp);
 
+printf("DA:%.20f %.20f %.20f\n",d_theta1,force.x,temp.x);
        /* part 1 of db_ij */
        v3_scale(&force,&force,sin_theta*2*h_cos*f_c_ik*frac/d2_h_cos2);
 

diff --git a/moldyn.h b/moldyn.h
index b2f7522..c5fed60 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -48,7 +48,6 @@ typedef struct s_linkcell {
        double x,y,z;           /* the actual cell lengthes */
        t_list *subcell;        /* pointer to the cell lists */
        int dnlc;               /* direct neighbour lists counter */
-       int countn;             /* amount of neighbours */
 } t_linkcell;
 
 #include "visual/visual.h"
@@ -127,6 +126,8 @@ typedef struct s_moldyn {
        u8 status;              /* general moldyn properties */
 
        t_random random;        /* random interface */
+
+       int debug;              /* debugging stuff, ignore */
 } t_moldyn;
 
 #define MOLDYN_STAT_PBX                        0x08    /* periodic boudaries in x */
@@ -261,6 +262,8 @@ typedef struct s_tersoff_mult_params {
 #define TRUE                           1
 #define FALSE                          0
 
+#define ACCEPTABLE_ERROR               1e-15
+
 /*
  *
  * phsical values / constants