X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=b6fa8d3df41dbf12e89c0596463a6cba4207255b;hb=88ea16029cc07a837977f737f4536fa3e881dcee;hp=aaab0deb5270c0157ebcbc4a07e58ddd12bf0943;hpb=ff1bf111f517d9fa69903cd56181c1e9e38665ff;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index aaab0de..b6fa8d3 100644 --- a/moldyn.c +++ b/moldyn.c @@ -78,6 +78,13 @@ int set_temperature(t_moldyn *moldyn,double t_ref) { return 0; } +int set_pressure(t_moldyn *moldyn,double p_ref) { + + moldyn->p_ref=p_ref; + + return 0; +} + int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) { moldyn->pt_scale=(ptype|ttype); @@ -93,16 +100,19 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) { moldyn->dim.y=y; moldyn->dim.z=z; + moldyn->volume=x*y*z; + if(visualize) { moldyn->vis.dim.x=x; moldyn->vis.dim.y=y; moldyn->vis.dim.z=z; } - printf("[moldyn] dimensions in A:\n"); + printf("[moldyn] dimensions in A and A^2 respectively:\n"); printf(" x: %f\n",moldyn->dim.x); printf(" y: %f\n",moldyn->dim.y); printf(" z: %f\n",moldyn->dim.z); + printf(" volume: %f\n",moldyn->volume); printf(" visualize simulation box: %s\n",visualize?"on":"off"); return 0; @@ -236,10 +246,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 +273,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 +292,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; } @@ -389,7 +400,7 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { count+=1; } } - if(count!=0) moldyn->t=(2.0*e)/(3.0*count*K_BOLTZMANN); + if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN); else return 0; /* no atoms involved in scaling! */ /* (temporary) hack for e,t = 0 */ @@ -422,6 +433,57 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) { return 0; } +int scale_volume(t_moldyn *moldyn) { + + t_atom *atom; + t_3dvec *dim,*vdim; + double virial,scale; + t_linkcell *lc; + int i; + + atom=moldyn->atom; + dim=&(moldyn->dim); + vdim=&(moldyn->vis.dim); + lc=&(moldyn->lc); + + for(i=0;icount;i++) + virial+=v3_norm(&(atom[i].virial)); + +printf("%f\n",virial); + /* get pressure from virial */ + moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial; + moldyn->p/=moldyn->volume; +printf("%f\n",moldyn->p/(ATM)); + + /* scale factor */ + if(moldyn->pt_scale&P_SCALE_BERENDSEN) + scale=3*sqrt(1-(moldyn->p_ref-moldyn->p)/moldyn->p_tc); + else + /* should actually never be used */ + scale=pow(moldyn->p/moldyn->p_ref,1.0/3.0); + +printf("scale = %f\n",scale); + /* actual scaling */ + dim->x*=scale; + dim->y*=scale; + dim->z*=scale; + if(vdim->x) vdim->x=dim->x; + if(vdim->y) vdim->y=dim->y; + if(vdim->z) vdim->z=dim->z; + moldyn->volume*=(scale*scale*scale); + + /* check whether we need a new linkcell init */ + if((dim->x/moldyn->cutoff!=lc->nx)|| + (dim->y/moldyn->cutoff!=lc->ny)|| + (dim->z/moldyn->cutoff!=lc->nx)) { + link_cell_shutdown(moldyn); + link_cell_init(moldyn); + } + + return 0; + +} + double get_e_kin(t_moldyn *moldyn) { int i; @@ -723,6 +785,12 @@ int moldyn_integrate(t_moldyn *moldyn) { /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) scale_velocity(moldyn,FALSE); + if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT)) +{ +printf("going to do p scale ...\n"); + scale_volume(moldyn); +printf("done\n"); +} /* check for log & visualization */ if(e) { @@ -801,14 +869,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 +910,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; @@ -851,13 +920,16 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset energy */ moldyn->energy=0.0; - + /* get energy and force of every atom */ for(i=0;ifunc1b(moldyn,&(itom[i])); @@ -971,9 +1043,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 +1306,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 +1578,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 +1629,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 +1651,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,20 +1683,21 @@ 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 ^ */ + } return 0; @@ -1635,20 +1713,44 @@ 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"); + printf("-----------\natom.r.x:\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;