X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=ed831aed93ec4861c71a383033af863ef1378074;hb=dece53fc37f9ebb52b33c9743333c213be2d6f26;hp=eb2c6068caae95cd7efac77c5dff1e7e0eb285b8;hpb=3ce21f4433c04699bb6f47f0332a6f0e3c462f5e;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index eb2c606..ed831ae 100644 --- a/moldyn.c +++ b/moldyn.c @@ -41,9 +41,9 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { int moldyn_shutdown(t_moldyn *moldyn) { + printf("[moldyn] shutdown\n"); moldyn_log_shutdown(moldyn); link_cell_shutdown(moldyn); - moldyn_log_shutdown(moldyn); rand_close(&(moldyn->random)); free(moldyn->atom); @@ -93,6 +93,13 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) { return 0; } +int set_nn_dist(t_moldyn *moldyn,double dist) { + + moldyn->nnd=dist; + + return 0; +} + int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) { if(x) @@ -171,9 +178,10 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer) { int moldyn_log_shutdown(t_moldyn *moldyn) { + printf("[moldyn] log shutdown\n"); if(moldyn->efd) close(moldyn->efd); if(moldyn->mfd) close(moldyn->mfd); - if(moldyn->visual) visual_tini(moldyn->visual); + if(&(moldyn->vis)) visual_tini(&(moldyn->vis)); return 0; } @@ -184,17 +192,15 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, int count; int ret; t_3dvec origin; - t_atom *atom; count=a*b*c; - atom=moldyn->atom; if(type==FCC) count*=4; if(type==DIAMOND) count*=8; - atom=malloc(count*sizeof(t_atom)); - if(atom==NULL) { + moldyn->atom=malloc(count*sizeof(t_atom)); + if(moldyn->atom==NULL) { perror("malloc (atoms)"); return -1; } @@ -203,10 +209,10 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, switch(type) { case FCC: - ret=fcc_init(a,b,c,lc,atom,&origin); + ret=fcc_init(a,b,c,lc,moldyn->atom,&origin); break; case DIAMOND: - ret=diamond_init(a,b,c,lc,atom,&origin); + ret=diamond_init(a,b,c,lc,moldyn->atom,&origin); break; default: printf("unknown lattice type (%02x)\n",type); @@ -222,15 +228,18 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, } moldyn->count=count; + printf("[moldyn] created lattice with %d atoms\n",count); while(count) { - atom[count-1].element=element; - atom[count-1].mass=mass; - atom[count-1].attr=attr; - atom[count-1].bnum=bnum; count-=1; + moldyn->atom[count].element=element; + moldyn->atom[count].mass=mass; + moldyn->atom[count].attr=attr; + moldyn->atom[count].bnum=bnum; + check_per_bound(moldyn,&(moldyn->atom[count].r)); } + return ret; } @@ -312,15 +321,15 @@ int thermal_init(t_moldyn *moldyn) { } /* velocity scaling */ - scale_velocity(moldyn); + scale_velocity(moldyn,VSCALE_INIT_EQUI); return 0; } -int scale_velocity(t_moldyn *moldyn) { +int scale_velocity(t_moldyn *moldyn,u8 type) { int i; - double e,c; + double e,scale; t_atom *atom; atom=moldyn->atom; @@ -329,17 +338,14 @@ int scale_velocity(t_moldyn *moldyn) { * - velocity scaling (E = 3/2 N k T), E: kinetic energy */ - if(moldyn->t==0.0) { - printf("[moldyn] no velocity scaling for T = 0 K\n"); - return -1; - } - e=0.0; for(i=0;icount;i++) e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v)); - c=sqrt((2.0*e)/(3.0*moldyn->count*K_BOLTZMANN*moldyn->t)); + scale=(1.5*moldyn->count*K_BOLTZMANN*moldyn->t)/e; + if(type&VSCALE_INIT_EQUI) scale*=2.0; /* equipartition theorem */ + scale=sqrt(scale); for(i=0;icount;i++) - v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c)); + v3_scale(&(atom[i].v),&(atom[i].v),scale); return 0; } @@ -525,7 +531,7 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { } } - lc->dnlc=count2; + lc->dnlc=count1; lc->countn=27; return count2; @@ -595,9 +601,9 @@ int moldyn_integrate(t_moldyn *moldyn) { t_3dvec p; t_moldyn_schedule *schedule; t_atom *atom; - int fd; char fb[128]; + double ds; schedule=&(moldyn->schedule); atom=moldyn->atom; @@ -614,10 +620,20 @@ 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 */ + 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) + 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"); + 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"); + /* zero absolute time */ moldyn->time=0.0; @@ -759,7 +775,7 @@ int potential_force_calc(t_moldyn *moldyn) { moldyn->energy=0.0; for(i=0;iattr&ATOM_ATTR_3BP)) continue; +printf("DEBUG: problem exists here ...\n"); link_cell_neighbour_index(moldyn, (btom->r.x+moldyn->dim.x/2)/lc->x, (btom->r.y+moldyn->dim.y/2)/lc->y, (btom->r.z+moldyn->dim.z/2)/lc->z, neighbourk); +printf("DEBUG: as you won't see that!\n"); for(k=0;kcountn;k++) { @@ -956,6 +974,20 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { * tersoff potential & force for 2 sorts of atoms */ +/* create mixed terms from parameters and set them */ +int tersoff_mult_complete_params(t_tersoff_mult_params *p) { + + printf("[moldyn] tersoff parameter completion\n"); + p->Smixed=sqrt(p->S[0]*p->S[1]); + p->Rmixed=sqrt(p->R[0]*p->R[1]); + p->Amixed=sqrt(p->A[0]*p->A[1]); + p->Bmixed=sqrt(p->B[0]*p->B[1]); + p->lambda_m=0.5*(p->lambda[0]+p->lambda[1]); + p->mu_m=0.5*(p->mu[0]+p->mu[1]); + + return 0; +} + /* tersoff 1 body part */ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {