From ade81aa2afb15f22e98ed9595ff303d4fedfe122 Mon Sep 17 00:00:00 2001 From: hackbard Date: Wed, 6 Dec 2006 01:51:41 +0000 Subject: [PATCH] need sleep, continue debugging soon! --- clean | 2 +- moldyn.c | 98 +++++++++++++++++++++++++------------------------------- moldyn.h | 5 ++- 3 files changed, 48 insertions(+), 57 deletions(-) diff --git a/clean b/clean index fc5f86d..ff74c9b 100755 --- 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 --- 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;schedschedule.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;idim.z/2)/lc->z, neighbour_i); - countn=lc->countn; dnlc=lc->dnlc; - for(j=0;jattr&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;kcountn;k++) { -// -// that=&(neighbour_j[k]); -// list_reset(that); -// -// if(that->start==NULL) -// continue; -// -// bc_ijk=(kdnlc)?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;kstart==NULL) continue; - bc_ijk=(kfunc3b(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 --- 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 -- 2.20.1