X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=73becdeed1a9b27220a714d55eed845d2ec71db2;hb=72df64eacc634e315f2205fc7bd2406223f92bf3;hp=6736e457b28b9b97fe1b5dd89c9093a15d7e1c30;hpb=b30d3d558354102cf5d0c06392961032b21bf9a9;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 6736e45..73becde 100644 --- 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; } @@ -155,38 +161,56 @@ int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) { return 0; } -int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer) { +int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) { + + strncpy(moldyn->vlsdir,dir,127); + + return 0; +} + +int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { + + char filename[128]; + int ret; switch(type) { case LOG_TOTAL_ENERGY: moldyn->ewrite=timer; - moldyn->efd=open(fb,O_WRONLY|O_CREAT|O_TRUNC); + snprintf(filename,127,"%s/energy",moldyn->vlsdir); + moldyn->efd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); if(moldyn->efd<0) { - perror("[moldyn] efd open"); + perror("[moldyn] energy log fd open"); return moldyn->efd; } dprintf(moldyn->efd,"# total energy log file\n"); break; case LOG_TOTAL_MOMENTUM: moldyn->mwrite=timer; - moldyn->mfd=open(fb,O_WRONLY|O_CREAT|O_TRUNC); + snprintf(filename,127,"%s/momentum",moldyn->vlsdir); + moldyn->mfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); if(moldyn->mfd<0) { - perror("[moldyn] mfd open"); + perror("[moldyn] momentum log fd open"); return moldyn->mfd; } dprintf(moldyn->efd,"# total momentum log file\n"); break; case SAVE_STEP: moldyn->swrite=timer; - strncpy(moldyn->sfb,fb,63); break; case VISUAL_STEP: moldyn->vwrite=timer; - strncpy(moldyn->vfb,fb,63); - visual_init(&(moldyn->vis),fb); + ret=visual_init(&(moldyn->vis),moldyn->vlsdir); + if(ret<0) { + printf("[moldyn] visual init failure\n"); + return ret; + } break; default: - printf("unknown log mechanism: %02x\n",type); + printf("[moldyn] unknown log mechanism: %02x\n",type); return -1; } @@ -212,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)"); @@ -238,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; } @@ -256,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; } @@ -371,8 +396,10 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { /* (temporary) hack for e,t = 0 */ if(e==0.0) { moldyn->t=0.0; - if(moldyn->t_ref!=0.0) + if(moldyn->t_ref!=0.0) { thermal_init(moldyn,equi_init); + return 0; + } else return 0; /* no scaling needed */ } @@ -384,13 +411,14 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { scale*=2.0; else if(moldyn->pt_scale&T_SCALE_BERENDSEN) - scale=1.0+moldyn->tau*(scale-1.0)/moldyn->t_tc; + scale=1.0+(scale-1.0)/moldyn->t_tc; scale=sqrt(scale); /* velocity scaling */ - for(i=0;icount;i++) + for(i=0;icount;i++) { if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB)) v3_scale(&(atom[i].v),&(atom[i].v),scale); + } return 0; } @@ -447,12 +475,6 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) { /* nn_dist is the nearest neighbour distance */ - if(moldyn->t==5.0) { - printf("[moldyn] i do not estimate timesteps below %f K!\n", - MOLDYN_CRITICAL_EST_TEMP); - return 23.42; - } - tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t); return tau; @@ -645,7 +667,7 @@ int moldyn_integrate(t_moldyn *moldyn) { t_moldyn_schedule *schedule; t_atom *atom; int fd; - char fb[128]; + char dir[128]; double ds; schedule=&(moldyn->schedule); @@ -681,7 +703,7 @@ int moldyn_integrate(t_moldyn *moldyn) { /* zero absolute time */ moldyn->time=0.0; - /* debugging, ignre */ + /* debugging, ignore */ moldyn->debug=0; /* executing the schedule */ @@ -703,14 +725,11 @@ int moldyn_integrate(t_moldyn *moldyn) { if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) scale_velocity(moldyn,FALSE); - /* increase absolute time */ - moldyn->time+=moldyn->tau; - /* check for log & visualization */ if(e) { if(!(i%e)) dprintf(moldyn->efd, - "%.15f %.45f %.45f %.45f\n", + "%f %f %f %f\n", moldyn->time,update_e_kin(moldyn), moldyn->energy, get_total_energy(moldyn)); @@ -719,15 +738,14 @@ int moldyn_integrate(t_moldyn *moldyn) { if(!(i%m)) { p=get_total_p(moldyn); dprintf(moldyn->mfd, - "%.15f %.45f\n",moldyn->time, - v3_norm(&p)); + "%f %f\n",moldyn->time,v3_norm(&p)); } } if(s) { if(!(i%s)) { - snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb, - moldyn->t,i*moldyn->tau); - fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT); + snprintf(dir,128,"%s/s-%07.f.save", + moldyn->vlsdir,moldyn->time); + fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT); if(fd<0) perror("[moldyn] save fd open"); else { write(fd,moldyn,sizeof(t_moldyn)); @@ -741,12 +759,15 @@ int moldyn_integrate(t_moldyn *moldyn) { if(!(i%v)) { visual_atoms(&(moldyn->vis),moldyn->time, moldyn->atom,moldyn->count); - printf("\rsched: %d, steps: %d, theta: %d", + printf("\rsched: %d, steps: %d, debug: %d", sched,i,moldyn->debug); fflush(stdout); } } + /* increase absolute time */ + moldyn->time+=moldyn->tau; + } /* check for hooks */ @@ -781,13 +802,16 @@ int velocity_verlet(t_moldyn *moldyn) { v3_add(&(atom[i].r),&(atom[i].r),&delta); v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass); v3_add(&(atom[i].r),&(atom[i].r),&delta); +//if(i==5) printf("v: %f %f %f\n",atom[i].r.x,(atom[i].r.x+moldyn->dim.x/2)/moldyn->lc.x,2*atom[i].r.x/moldyn->dim.x); check_per_bound(moldyn,&(atom[i].r)); +//if(i==5) printf("n: %f %f %f\n",atom[i].r.x,(atom[i].r.x+moldyn->dim.x/2)/moldyn->lc.x,2*atom[i].r.x/moldyn->dim.x); /* velocities */ v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass); v3_add(&(atom[i].v),&(atom[i].v),&delta); } +//moldyn_bc_check(moldyn); /* neighbour list update */ link_cell_update(moldyn); @@ -819,7 +843,6 @@ int potential_force_calc(t_moldyn *moldyn) { t_linkcell *lc; t_list neighbour_i[27]; t_list neighbour_i2[27]; - //t_list neighbour_j[27]; t_list *this,*that; u8 bc_ij,bc_ik; int dnlc; @@ -841,94 +864,99 @@ int potential_force_calc(t_moldyn *moldyn) { if(itom[i].attr&ATOM_ATTR_1BP) moldyn->func1b(moldyn,&(itom[i])); + if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP))) + continue; + /* 2 body pair potential/force */ - if(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) { - link_cell_neighbour_index(moldyn, - (itom[i].r.x+moldyn->dim.x/2)/lc->x, - (itom[i].r.y+moldyn->dim.y/2)/lc->y, - (itom[i].r.z+moldyn->dim.z/2)/lc->z, - neighbour_i); + link_cell_neighbour_index(moldyn, + (itom[i].r.x+moldyn->dim.x/2)/lc->x, + (itom[i].r.y+moldyn->dim.y/2)/lc->y, + (itom[i].r.z+moldyn->dim.z/2)/lc->z, + neighbour_i); - dnlc=lc->dnlc; + dnlc=lc->dnlc; - for(j=0;j<27;j++) { + for(j=0;j<27;j++) { - this=&(neighbour_i[j]); - list_reset(this); + this=&(neighbour_i[j]); + list_reset(this); - if(this->start==NULL) - continue; + if(this->start==NULL) + continue; - bc_ij=(jcurrent->data; + do { + jtom=this->current->data; - if(jtom==&(itom[i])) - continue; + if(jtom==&(itom[i])) + continue; - if((jtom->attr&ATOM_ATTR_2BP)& - (itom[i].attr&ATOM_ATTR_2BP)) - moldyn->func2b(moldyn, - &(itom[i]), - jtom, - bc_ij); + if((jtom->attr&ATOM_ATTR_2BP)& + (itom[i].attr&ATOM_ATTR_2BP)) { + moldyn->func2b(moldyn, + &(itom[i]), + jtom, + bc_ij); + } - /* 3 body potential/force */ + /* 3 body potential/force */ - if(!(itom[i].attr&ATOM_ATTR_3BP)|| - !(jtom->attr&ATOM_ATTR_3BP)) - continue; + if(!(itom[i].attr&ATOM_ATTR_3BP)|| + !(jtom->attr&ATOM_ATTR_3BP)) + continue; - /* copy the neighbour lists */ - memcpy(neighbour_i2,neighbour_i, - 27*sizeof(t_list)); + /* copy the neighbour lists */ + memcpy(neighbour_i2,neighbour_i, + 27*sizeof(t_list)); - /* get neighbours of i */ - for(k=0;k<27;k++) { + /* get neighbours of i */ + for(k=0;k<27;k++) { - that=&(neighbour_i2[k]); - list_reset(that); + that=&(neighbour_i2[k]); + list_reset(that); - if(that->start==NULL) - continue; + if(that->start==NULL) + continue; - bc_ik=(kcurrent->data; + ktom=that->current->data; - if(!(ktom->attr&ATOM_ATTR_3BP)) - continue; + if(!(ktom->attr&ATOM_ATTR_3BP)) + continue; - if(ktom==jtom) - continue; + if(ktom==jtom) + continue; - if(ktom==&(itom[i])) - continue; + if(ktom==&(itom[i])) + continue; - moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ik|bc_ij); + moldyn->func3b(moldyn, + &(itom[i]), + jtom, + ktom, + bc_ik|bc_ij); - } while(list_next(that)!=\ - L_NO_NEXT_ELEMENT); + } while(list_next(that)!=\ + L_NO_NEXT_ELEMENT); + + } - } - - } while(list_next(this)!=L_NO_NEXT_ELEMENT); - /* 2bp post function */ if(moldyn->func2b_post) { -printf("DEBUG: vor 2bp post\n"); moldyn->func2b_post(moldyn, &(itom[i]), jtom,bc_ij); -printf("DEBUG: nach 2bp post\n"); } - - } + + } while(list_next(this)!=L_NO_NEXT_ELEMENT); + } + } return 0; @@ -945,9 +973,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; @@ -1051,8 +1079,8 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) { p->mu_m=0.5*(p->mu[0]+p->mu[1]); printf("[moldyn] tersoff mult parameter info:\n"); - printf(" S (m) | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed); - printf(" R (m) | %.12f | %.12f | %.12f\n",p->R[0],p->R[1],p->Rmixed); + printf(" S (A) | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed); + printf(" R (A) | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed); printf(" A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV); printf(" B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV); printf(" lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1], @@ -1119,6 +1147,9 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* clear 3bp and 2bp post run */ exchange->run3bp=0; exchange->run2bp_post=0; + + /* reset S > r > R mark */ + exchange->d_ij_between_rs=0; /* * calc of 2bp contribution of V_ij and dV_ij/ji @@ -1187,11 +1218,11 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* f_r_ij = f_r_ji, df_r_ij = df_r_ji */ f_r=A*exp(-lambda*d_ij); - df_r=-lambda*f_r/d_ij; + df_r=lambda*f_r/d_ij; /* 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 > R */ + exchange->d_ij_between_rs=1; } /* add forces of 2bp (ij, ji) contribution @@ -1281,40 +1315,57 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { /* Vij and dVij */ zeta=exchange->zeta_ij; - tmp=betaini*pow(zeta,ni-1.0); /* beta^n * zeta^n-1 */ - b=(1+zeta*tmp); /* 1 + beta^n * zeta^n */ - db=chi*pow(b,-1.0/(2*ni)-1); /* chi * (...)^(-1/2n - 1) */ - b=db*b; /* b_ij */ - db*=-0.5*tmp; /* db_ij */ - v3_scale(&force,&(exchange->dzeta_ij),f_a*db); - v3_scale(&temp,dist_ij,df_a*b); - v3_add(&force,&force,&temp); - v3_scale(&force,&force,f_c); + if(zeta==0.0) { + moldyn->debug++; /* just for debugging ... */ + db=0.0; + b=chi; + v3_scale(&force,dist_ij,df_a*b*f_c); + } + else { + tmp=betaini*pow(zeta,ni-1.0); /* beta^n * zeta^n-1 */ + b=(1+zeta*tmp); /* 1 + beta^n zeta^n */ + db=chi*pow(b,-1.0/(2*ni)-1); /* x(...)^(-1/2n - 1) */ + b=db*b; /* b_ij */ + db*=-0.5*tmp; /* db_ij */ + v3_scale(&force,&(exchange->dzeta_ij),f_a*db); + v3_scale(&temp,dist_ij,df_a*b); + v3_add(&force,&force,&temp); + v3_scale(&force,&force,f_c); + } v3_scale(&temp,dist_ij,df_c*b*f_a); v3_add(&force,&force,&temp); + v3_scale(&force,&force,-0.5); + + /* add force */ + v3_add(&(ai->f),&(ai->f),&force); /* add energy of 3bp sum */ moldyn->energy+=(0.5*f_c*b*f_a); - /* add force (sub, as F = - dVij) */ - v3_sub(&(ai->f),&(ai->f),&force); - /* dVji */ zeta=exchange->zeta_ji; - tmp=betajnj*pow(zeta,nj-1.0); /* beta^n * zeta^n-1 */ - b=(1+zeta*tmp); /* 1 + beta^n * zeta^n */ - db=chi*pow(b,-1.0/(2*nj)-1); /* chi * (...)^(-1/2n - 1) */ - b=db*b; /* b_ij */ - db*=-0.5*tmp; /* db_ij */ - v3_scale(&force,&(exchange->dzeta_ji),f_a*db); - v3_scale(&temp,dist_ij,df_a*b); - v3_add(&force,&force,&temp); - v3_scale(&force,&force,f_c); + if(zeta==0.0) { + moldyn->debug++; + b=chi; + v3_scale(&force,dist_ij,df_a*b*f_c); + } + else { + tmp=betajnj*pow(zeta,nj-1.0); /* beta^n * zeta^n-1 */ + b=(1+zeta*tmp); /* 1 + beta^n zeta^n */ + db=chi*pow(b,-1.0/(2*nj)-1); /* x(...)^(-1/2n - 1) */ + b=db*b; /* b_ij */ + db*=-0.5*tmp; /* db_ij */ + v3_scale(&force,&(exchange->dzeta_ji),f_a*db); + v3_scale(&temp,dist_ij,df_a*b); + v3_add(&force,&force,&temp); + v3_scale(&force,&force,f_c); + } v3_scale(&temp,dist_ij,df_c*b*f_a); v3_add(&force,&force,&temp); + v3_scale(&force,&force,-0.5); - /* add force (sub, as F = - dVji) */ - v3_sub(&(ai->f),&(ai->f),&force); + /* add force */ + v3_add(&(ai->f),&(ai->f),&force); return 0; } @@ -1457,7 +1508,8 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { s_r=S-R; arg=M_PI*(d_ik-R)/s_r; f_c_ik=0.5+0.5*cos(arg); - df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); + //df_c_ik=-0.5*sin(arg)*(M_PI/(s_r*d_ik)); /* MARK */ + df_c_ik=0.5*sin(arg)*(M_PI/(s_r*d_ik)); /* zeta_ij */ exchange->zeta_ij+=f_c_ik*g; @@ -1507,13 +1559,13 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { c2d2=exchange->cj2dj2; /* cosine of theta_jik by scalaproduct */ - rr=v3_scalar_product(&dist_ij,&dist_jk); /* times -1 */ + rr=-v3_scalar_product(&dist_ij,&dist_jk); /* -1, as ij -> ji */ dd=d_ij*d_jk; cos_theta=rr/dd; /* d_costheta */ - d_costheta1=1.0/(d_jk*d_ij); - d_costheta2=cos_theta/(d_ij*d_ij); /* in fact -cos(), but ^ */ + d_costheta1=1.0/dd; + d_costheta2=cos_theta/(d_ij*d_ij); /* some usefull values */ h_cos=(h-cos_theta); @@ -1529,6 +1581,9 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { v3_add(&temp1,&temp1,&temp2); v3_scale(&temp1,&temp1,-2.0*frac*h_cos/d2_h_cos2); /* dg */ + /* store dg in temp2 and use it for dVjk later */ + v3_copy(&temp2,&temp1); + /* f_c_jk + {d,}zeta contribution (df_c_jk = 0) */ dzeta=&(exchange->dzeta_ji); if(d_jkzeta_ji+=f_c_jk*g; - /* dzeta_ij */ + /* dzeta_ji */ v3_scale(&temp1,&temp1,f_c_jk); v3_add(dzeta,dzeta,&temp1); } @@ -1558,19 +1613,72 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) { /* dV_jk stuff | add force contribution on atom i immediately */ if(exchange->d_ij_between_rs) { zeta=f_c*g; - v3_scale(&temp1,&temp1,f_c); - v3_scale(&temp2,&dist_ij,df_c); - v3_add(&temp1,&temp1,&temp2); + v3_scale(&temp1,&temp2,f_c); + v3_scale(&temp2,&dist_ij,df_c*g); + v3_add(&temp2,&temp2,&temp1); /* -> dzeta_jk in temp2 */ } else { zeta=g; - // dzeta_jk is simply dg, which is temp1 + // dzeta_jk is simply dg, which is stored in temp2 } /* betajnj * zeta_jk ^ nj-1 */ tmp=exchange->betajnj*pow(zeta,(n-1.0)); - tmp=-chi/2.0*pow(1+tmp*zeta,-1.0/(2.0*n)-1)*tmp; - v3_scale(&temp1,&temp1,tmp*B*exp(-mu*d_jk)*f_c_jk); - v3_add(&(ai->f),&(ai->f),&temp1); /* -1 skipped in f_a calc ^ */ + tmp=-chi/2.0*pow((1+tmp*zeta),(-1.0/(2.0*n)-1))*tmp; + v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5); + v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */ + /* scaled with 0.5 ^ */ + } + + return 0; +} + + +/* + * debugging / critical check functions + */ + +int moldyn_bc_check(t_moldyn *moldyn) { + + t_atom *atom; + t_3dvec *dim; + int i; +double x; +u8 byte; +int j,k; + + atom=moldyn->atom; + dim=&(moldyn->dim); +x=dim->x/2; + + for(i=0;icount;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,dim->x/2); + printf("diagnostic:\n"); + for(j=0;j<8;j++) { + memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1); + for(k=0;k<8;k++) + printf("%d%c", + ((byte)&(1<=dim->y/2||-atom[i].r.y>dim->y/2) + printf("FATAL: atom %d: y: %.20f (%.20f)\n", + 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,dim->z/2); } return 0;