From: hackbard Date: Thu, 30 Nov 2006 10:46:08 +0000 (+0000) Subject: vel scaling fixed/changed X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=180ff86ae35ca0ea6ee857e92e56120e65315d17;p=physik%2Fposic.git vel scaling fixed/changed --- diff --git a/moldyn.c b/moldyn.c index de7bbbf..d52fcd3 100644 --- a/moldyn.c +++ b/moldyn.c @@ -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) @@ -221,6 +228,7 @@ 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) { moldyn->atom[count-1].element=element; @@ -311,15 +319,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; @@ -328,17 +336,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; } @@ -594,9 +599,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; @@ -617,6 +622,17 @@ int moldyn_integrate(t_moldyn *moldyn) { /* 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; @@ -758,7 +774,6 @@ int potential_force_calc(t_moldyn *moldyn) { moldyn->energy=0.0; for(i=0;i