From b3f72740276662c09ff43c294b04409a5bb2e013 Mon Sep 17 00:00:00 2001 From: hackbard Date: Fri, 8 Feb 2008 01:30:25 +0100 Subject: [PATCH] fixed missing init for new atoms (especially force!) + increased insert relaxation time + added some more critical debugging checks --- config.h | 2 +- moldyn.c | 48 ++++++++++++++++++++++++++++++++++-------------- sic.c | 6 ++++++ 3 files changed, 41 insertions(+), 15 deletions(-) diff --git a/config.h b/config.h index f0919b7..f43f006 100644 --- a/config.h +++ b/config.h @@ -41,7 +41,7 @@ #define INS_R_C 1.5 #define INS_DELTA_TC 1.0 -#define INS_RELAX 50 +#define INS_RELAX 100 #define INS_TAU 1.0 // postrun diff --git a/moldyn.c b/moldyn.c index ec96712..94beb04 100644 --- a/moldyn.c +++ b/moldyn.c @@ -595,7 +595,7 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, int count; atom=moldyn->atom; - count=(moldyn->count)++; + count=(moldyn->count)++; // asshole style! ptr=realloc(atom,(count+1)*sizeof(t_atom)); if(!ptr) { @@ -605,6 +605,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, moldyn->atom=ptr; atom=moldyn->atom; + + /* initialize new atom */ + memset(&(atom[count]),0,sizeof(t_atom)); atom[count].r=*r; atom[count].v=*v; atom[count].element=element; @@ -956,6 +959,8 @@ double pressure_calc(t_moldyn *moldyn) { int average_reset(t_moldyn *moldyn) { + printf("[moldyn] average reset\n"); + /* update skip value */ moldyn->avg_skip=moldyn->total_steps; @@ -1369,7 +1374,7 @@ int link_cell_update(t_moldyn *moldyn) { i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x); j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y); k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z); - + #ifdef STATIC_LISTS p=0; while(lc->subcell[i+j*nx+k*nx*ny][p]!=0) @@ -1419,6 +1424,10 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k, count2=27; a=nx*ny; + if(i>=nx||j>=ny||k>=nz) + printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n", + i,nx,j,ny,k,nz); + cell[0]=lc->subcell[i+j*nx+k*a]; for(ci=-1;ci<=1;ci++) { bx=0; @@ -1569,14 +1578,14 @@ int moldyn_integrate(t_moldyn *moldyn) { /* 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"); + printf("[moldyn] WARNING: cutoff > 0.5 x dim.x\n"); if(moldyn->cutoff>0.5*moldyn->dim.y) - printf("[moldyn] warning: cutoff > 0.5 x dim.y\n"); + printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n"); if(moldyn->cutoff>0.5*moldyn->dim.z) - printf("[moldyn] warning: cutoff > 0.5 x dim.z\n"); + printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n"); ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass; if(ds>0.05*moldyn->nnd) - printf("[moldyn] warning: forces too high / tau too small!\n"); + printf("[moldyn] WARNING: forces too high / tau too small!\n"); /* zero absolute time */ moldyn->time=0.0; @@ -2085,14 +2094,21 @@ int potential_force_calc(t_moldyn *moldyn) { } #endif - /* calculate global virial */ + /* some postprocessing */ for(i=0;igvir.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x; - moldyn->gvir.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y; - moldyn->gvir.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z; - moldyn->gvir.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x; - moldyn->gvir.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x; - moldyn->gvir.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y; + /* calculate global virial */ + moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x; + moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y; + moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z; + moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x; + moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x; + moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y; + + /* check forces regarding the given timestep */ + if(v3_norm(&(itom[i].f))>\ + 0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square) + printf("[moldyn] WARNING: pfc (high force: atom %d)\n", + i); } return 0; @@ -2308,6 +2324,9 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { slots=moldyn->cutoff/dr; o=2*slots; + if(slots*dr<=moldyn->cutoff) + printf("[moldyn] WARNING: pcc (low #slots)\n"); + printf("[moldyn] pair correlation calc info:\n"); printf(" time: %f\n",moldyn->time); printf(" count: %d\n",moldyn->count); @@ -2389,8 +2408,9 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { /* should never happen but it does 8) - * related to -ffloat-store problem! */ if(s>=slots) { - printf("[moldyn] WARNING pcc (%d/%d)\n", + printf("[moldyn] WARNING: pcc (%d/%d)", s,slots); + printf("\n"); s=slots-1; } diff --git a/sic.c b/sic.c index a7613e2..fb20a56 100644 --- a/sic.c +++ b/sic.c @@ -122,6 +122,9 @@ int sic_hook(void *moldyn,void *hook_params) { hp=hook_params; md=moldyn; + tau=1.0; + steps=0; + /* switch on t scaling */ if(md->schedule.count==0) set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0); @@ -438,6 +441,9 @@ int main(int argc,char **argv) { moldyn_set_log(&md,SAVE_STEP,LOG_S); moldyn_set_log(&md,CREATE_REPORT,0); + /* next neighbour distance for critical checking */ + set_nn_dist(&md,0.25*ALBE_LC_SI*sqrt(3.0)); + /* * let's do the actual md algorithm now * -- 2.20.1