X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=73becdeed1a9b27220a714d55eed845d2ec71db2;hb=72df64eacc634e315f2205fc7bd2406223f92bf3;hp=aaab0deb5270c0157ebcbc4a07e58ddd12bf0943;hpb=ff1bf111f517d9fa69903cd56181c1e9e38665ff;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index aaab0de..73becde 100644 --- a/moldyn.c +++ b/moldyn.c @@ -236,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)"); @@ -262,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; } @@ -280,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; } @@ -801,14 +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); +//moldyn_bc_check(moldyn); /* neighbour list update */ link_cell_update(moldyn); @@ -840,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; @@ -971,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; @@ -1234,7 +1236,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { s_r=S-R; arg=M_PI*(d_ij-R)/s_r; f_c=0.5+0.5*cos(arg); - df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij)); + //df_c=-0.5*sin(arg)*(M_PI/(s_r*d_ij)); /* MARK! */ + df_c=0.5*sin(arg)*(M_PI/(s_r*d_ij)); /* two body contribution (ij, ji) */ v3_scale(&force,&dist_ij,-df_c*f_r-df_r*f_c); /* tell 3bp that S > r > R */ @@ -1505,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; @@ -1555,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); @@ -1577,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); } @@ -1606,19 +1613,19 @@ 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*0.5); - 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 ^ */ } @@ -1635,20 +1642,43 @@ 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) + 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*1e10,dim->x/2*1e10); + 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*1e10,dim->y/2*1e10); + 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*1e10,dim->z/2*1e10); + i,atom[i].r.z,dim->z/2); } return 0;