From: hackbard Date: Fri, 26 Jan 2007 15:36:19 +0000 (+0000) Subject: pt scaling tests, though pressure estimation still a mess! X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=91e4a00285261865ffce4b9553d153b562c80ee6;p=physik%2Fposic.git pt scaling tests, though pressure estimation still a mess! --- diff --git a/moldyn.c b/moldyn.c index 881f0b6..0d5027c 100644 --- a/moldyn.c +++ b/moldyn.c @@ -73,7 +73,7 @@ int set_temperature(t_moldyn *moldyn,double t_ref) { moldyn->t_ref=t_ref; - printf("[moldyn] temperature: %f\n",moldyn->t_ref); + printf("[moldyn] temperature [K]: %f\n",moldyn->t_ref); return 0; } @@ -82,7 +82,7 @@ int set_pressure(t_moldyn *moldyn,double p_ref) { moldyn->p_ref=p_ref; - printf("[moldyn] pressure: %f\n",moldyn->p_ref); + printf("[moldyn] pressure [atm]: %f\n",moldyn->p_ref/ATM); return 0; } @@ -344,15 +344,18 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, switch(type) { case CUBIC: + set_nn_dist(moldyn,lc); origin.x=0.5*lc; origin.y=0.5*lc; origin.z=0.5*lc; ret=cubic_init(a,b,c,lc,atom,&origin); break; case FCC: + set_nn_dist(moldyn,0.5*sqrt(2.0)*lc); ret=fcc_init(a,b,c,lc,atom,NULL); break; case DIAMOND: + set_nn_dist(moldyn,0.25*sqrt(3.0)*lc); ret=diamond_init(a,b,c,lc,atom,&origin); break; default: @@ -720,7 +723,7 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { scale_dim(moldyn,scale,TRUE,0,0); scale_atoms(moldyn,scale,TRUE,0,0); link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); tp->x=(moldyn->energy-u)/moldyn->dv; p=tp->x*tp->x; @@ -734,7 +737,7 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { scale_dim(moldyn,scale,0,TRUE,0); scale_atoms(moldyn,scale,0,TRUE,0); link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); tp->y=(moldyn->energy-u)/moldyn->dv; p+=tp->y*tp->y; @@ -748,7 +751,7 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { scale_dim(moldyn,scale,0,0,TRUE); scale_atoms(moldyn,scale,0,0,TRUE); link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); tp->z=(moldyn->energy-u)/moldyn->dv; p+=tp->z*tp->z; @@ -765,7 +768,7 @@ printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x); scale_dim(moldyn,scale,1,1,1); scale_atoms(moldyn,scale,1,1,1); link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x); @@ -779,7 +782,7 @@ printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x); moldyn->energy=u; link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); return sqrt(p); } @@ -820,61 +823,48 @@ int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { int scale_volume(t_moldyn *moldyn) { - t_atom *atom; t_3dvec *dim,*vdim; - double scale,v; - t_virial virial; + double scale; t_linkcell *lc; - int i; - atom=moldyn->atom; - dim=&(moldyn->dim); vdim=&(moldyn->vis.dim); + dim=&(moldyn->dim); lc=&(moldyn->lc); - memset(&virial,0,sizeof(t_virial)); + /* scaling factor */ + if(moldyn->pt_scale&P_SCALE_BERENDSEN) { + scale=1.0-(moldyn->p_ref-moldyn->p)/moldyn->p_tc; + scale=pow(scale,ONE_THIRD); + } + else { + scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD); + } +moldyn->debug=scale; - for(i=0;icount;i++) { - virial.xx+=atom[i].virial.xx; - virial.yy+=atom[i].virial.yy; - virial.zz+=atom[i].virial.zz; - virial.xy+=atom[i].virial.xy; - virial.xz+=atom[i].virial.xz; - virial.yz+=atom[i].virial.yz; + /* scale the atoms and dimensions */ + scale_atoms(moldyn,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,scale,TRUE,TRUE,TRUE); + + /* visualize dimensions */ + if(vdim->x!=0) { + vdim->x=dim->x; + vdim->y=dim->y; + vdim->z=dim->z; } - /* just a guess so far ... */ - v=virial.xx+virial.yy+virial.zz; - -printf("%f\n",v); - /* get pressure from virial */ - moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+ONE_THIRD*v; - moldyn->p/=moldyn->volume; -printf("%f | %f\n",moldyn->p/(ATM),moldyn->p_ref/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)) { + /* recalculate scaled volume */ + moldyn->volume=dim->x*dim->y*dim->z; + + /* adjust/reinit linkcell */ + if(((int)(dim->x/moldyn->cutoff)!=lc->nx)|| + ((int)(dim->y/moldyn->cutoff)!=lc->ny)|| + ((int)(dim->z/moldyn->cutoff)!=lc->nx)) { link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); + } else { + lc->x*=scale; + lc->y*=scale; + lc->z*=scale; } return 0; @@ -939,7 +929,7 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) { /* linked list / cell method */ -int link_cell_init(t_moldyn *moldyn) { +int link_cell_init(t_moldyn *moldyn,u8 vol) { t_linkcell *lc; int i; @@ -960,7 +950,7 @@ int link_cell_init(t_moldyn *moldyn) { if(lc->cells<27) printf("[moldyn] FATAL: less then 27 subcells!\n"); - printf("[moldyn] initializing linked cells (%d)\n",lc->cells); + if(vol) printf("[moldyn] initializing linked cells (%d)\n",lc->cells); for(i=0;icells;i++) list_init_f(&(lc->subcell[i])); @@ -1137,7 +1127,7 @@ int moldyn_integrate(t_moldyn *moldyn) { atom=moldyn->atom; /* initialize linked cell method */ - link_cell_init(moldyn); + link_cell_init(moldyn,VERBOSE); /* logging & visualization */ e=moldyn->ewrite; @@ -1152,10 +1142,8 @@ int moldyn_integrate(t_moldyn *moldyn) { /* energy scaling factor */ energy_scale=moldyn->count*EV; -printf("debug: %f\n",moldyn->atom[0].f.x); /* calculate initial forces */ potential_force_calc(moldyn); -printf("debug: %f\n",moldyn->atom[0].f.x); /* some stupid checks before we actually start calculating bullshit */ if(moldyn->cutoff>0.5*moldyn->dim.x) @@ -1192,17 +1180,18 @@ printf("debug: %f\n",moldyn->atom[0].f.x); /* integration step */ moldyn->integrate(moldyn); + /* calculate kinetic energy, temperature and pressure */ + update_e_kin(moldyn); + temperature_calc(moldyn); + pressure_calc(moldyn); + //thermodynamic_pressure_calc(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)) scale_volume(moldyn); - update_e_kin(moldyn); - temperature_calc(moldyn); - pressure_calc(moldyn); - //thermodynamic_pressure_calc(moldyn); - /* check for log & visualization */ if(e) { if(!(i%e)) @@ -1237,8 +1226,9 @@ printf("debug: %f\n",moldyn->atom[0].f.x); if(!(i%v)) { visual_atoms(&(moldyn->vis),moldyn->time, moldyn->atom,moldyn->count); - printf("\rsched: %d, steps: %d, debug: %f", - sched->count,i,moldyn->p/ATM); + printf("\rsched: %d, steps: %d, T: %f, P: %f V: %f", + sched->count,i, + moldyn->t,moldyn->p/ATM,moldyn->volume); fflush(stdout); } } diff --git a/moldyn.h b/moldyn.h index 0dbc1f9..f2ccdb5 100644 --- a/moldyn.h +++ b/moldyn.h @@ -313,6 +313,9 @@ typedef struct s_tersoff_mult_params { #define TRUE 1 #define FALSE 0 +#define VERBOSE 1 +#define QUIET 0 + /* * * phsical values / constants @@ -428,7 +431,7 @@ t_3dvec get_total_p(t_moldyn *moldyn); double estimate_time_step(t_moldyn *moldyn,double nn_dist); -int link_cell_init(t_moldyn *moldyn); +int link_cell_init(t_moldyn *moldyn,u8 vol); int link_cell_update(t_moldyn *moldyn); int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell); int link_cell_shutdown(t_moldyn *moldyn); diff --git a/sic.c b/sic.c index 7e1249c..2a2c7bc 100644 --- a/sic.c +++ b/sic.c @@ -146,26 +146,15 @@ int main(int argc,char **argv) { // ATOM_ATTR_2BP|ATOM_ATTR_HB, // &r,&v); - /* setting a nearest neighbour distance for the moldyn checks */ - //set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */ - set_nn_dist(&md,LC_SI); - - /* set temperature */ - //set_temperature(&md,273.0+1410.0); - //set_temperature(&md,273.0+450.0); - //set_temperature(&md,273.0); - //set_temperature(&md,1.0); - //set_temperature(&md,0.0); + /* set temperature & pressure */ set_temperature(&md,atof(argv[2])+273.0); - - /* set pressure */ set_pressure(&md,ATM); /* set p/t scaling */ - //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0, - // T_SCALE_BERENDSEN,100.0); + set_pt_scale(&md,P_SCALE_BERENDSEN,0.001, + T_SCALE_BERENDSEN,100.0); //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0); - //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,0,0); + //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,0,0); /* initial thermal fluctuations of particles (in equilibrium) */ thermal_init(&md,TRUE); diff --git a/visual/visual.c b/visual/visual.c index a76c608..d6ce369 100644 --- a/visual/visual.c +++ b/visual/visual.c @@ -87,7 +87,8 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) { /* script file update */ dprintf(v->fd,"load xyz %s\n",file); - dprintf(v->fd,"spacefill 100\n"); + dprintf(v->fd,"spacefill\n"); + //dprintf(v->fd,"spacefill 100\n"); dprintf(v->fd,"rotate x 100\n"); dprintf(v->fd,"rotate y 10\n"); dprintf(v->fd,"set ambient 20\n");