X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=e0881de5db49b1669c5a0b02641e66760dc957a0;hb=e2c5d8b0eb4ce6faeb48830634eef522dbdb52b0;hp=e554a2fc2cd62a7a33ee4228dc7c10a574584893;hpb=92ef07d77a4c879527180224acea73a3f6564497;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index e554a2f..e0881de 100644 --- a/moldyn.c +++ b/moldyn.c @@ -218,6 +218,14 @@ int set_potential_params(t_moldyn *moldyn,void *params) { return 0; } +int set_avg_skip(t_moldyn *moldyn,int skip) { + + printf("[moldyn] skip %d steps before starting average calc\n",skip); + moldyn->avg_skip=skip; + + return 0; +} + int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) { strncpy(moldyn->vlsdir,dir,127); @@ -404,11 +412,14 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { if(moldyn->rfd) { dprintf(moldyn->rfd,report_end); close(moldyn->rfd); - snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir); + snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1", + moldyn->vlsdir); system(sc); - snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir); + snprintf(sc,255,"cd %s && pdflatex report >/dev/null 2>&1", + moldyn->vlsdir); system(sc); - snprintf(sc,255,"cd %s && dvipdf report",moldyn->vlsdir); + snprintf(sc,255,"cd %s && dvipdf report >/dev/null 2>&1", + moldyn->vlsdir); system(sc); } if(&(moldyn->vis)) visual_tini(&(moldyn->vis)); @@ -501,6 +512,9 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, check_per_bound(moldyn,&(atom[ret].r)); } + /* update total system mass */ + total_mass_calc(moldyn); + return ret; } @@ -642,6 +656,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, atom[count].tag=count; atom[count].attr=attr; + /* update total system mass */ + total_mass_calc(moldyn); + return 0; } @@ -702,13 +719,23 @@ int thermal_init(t_moldyn *moldyn,u8 equi_init) { return 0; } +double total_mass_calc(t_moldyn *moldyn) { + + int i; + + moldyn->mass=0.0; + + for(i=0;icount;i++) + moldyn->mass+=moldyn->atom[i].mass; + + return moldyn->mass; +} + double temperature_calc(t_moldyn *moldyn) { /* assume up to date kinetic energy, which is 3/2 N k_B T */ moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count); - moldyn->t_sum+=moldyn->t; - moldyn->mean_t=moldyn->t_sum/moldyn->total_steps; return moldyn->t; } @@ -791,7 +818,7 @@ double pressure_calc(t_moldyn *moldyn) { /* * PV = NkT + - * W = 1/3 sum_i f_i r_i + * with W = 1/3 sum_i f_i r_i (- skipped!) * virial = sum_i f_i r_i * * => P = (2 Ekin + virial) / (3V) @@ -803,26 +830,98 @@ double pressure_calc(t_moldyn *moldyn) { v+=(virial->xx+virial->yy+virial->zz); } - /* virial sum and mean virial */ - moldyn->virial_sum+=v; - moldyn->mean_v=moldyn->virial_sum/moldyn->total_steps; - - /* assume up to date kinetic energy */ - moldyn->p=2.0*moldyn->ekin+moldyn->mean_v; - moldyn->p/=(3.0*moldyn->volume); - moldyn->p_sum+=moldyn->p; - moldyn->mean_p=moldyn->p_sum/moldyn->total_steps; + /* virial sum and average virial */ + if(moldyn->total_steps>=moldyn->avg_skip) { + moldyn->virial_sum+=v; + moldyn->virial_avg=moldyn->virial_sum/ + (moldyn->total_steps+1-moldyn->avg_skip); + moldyn->p=2.0*moldyn->k_avg+moldyn->virial_avg; + moldyn->p/=(3.0*moldyn->volume); + moldyn->p_sum+=moldyn->p; + moldyn->p_avg=moldyn->p_sum/ + (moldyn->total_steps+1-moldyn->avg_skip); + } /* pressure from 'absolute coordinates' virial */ virial=&(moldyn->virial); v=virial->xx+virial->yy+virial->zz; moldyn->gp=2.0*moldyn->ekin+v; moldyn->gp/=(3.0*moldyn->volume); - moldyn->gp_sum+=moldyn->gp; - moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps; + if(moldyn->total_steps>=moldyn->avg_skip) { + moldyn->gp_sum+=moldyn->gp; + moldyn->gp_avg=moldyn->gp_sum/ + (moldyn->total_steps+1-moldyn->avg_skip); + } return moldyn->p; -} +} + +int average_and_fluctuation_calc(t_moldyn *moldyn) { + + if(moldyn->total_stepsavg_skip) + return 0; + + /* assume up to date energies, temperature, pressure etc */ + + /* kinetic energy */ + moldyn->k_sum+=moldyn->ekin; + moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin); + moldyn->k_avg=moldyn->k_sum/(moldyn->total_steps+1-moldyn->avg_skip); + moldyn->k2_avg=moldyn->k2_sum/(moldyn->total_steps+1-moldyn->avg_skip); + moldyn->dk2_avg=moldyn->k2_avg-(moldyn->k_avg*moldyn->k_avg); + + /* potential energy */ + moldyn->v_sum+=moldyn->energy; + moldyn->v2_sum+=(moldyn->energy*moldyn->energy); + moldyn->v_avg=moldyn->v_sum/(moldyn->total_steps+1-moldyn->avg_skip); + moldyn->v2_avg=moldyn->v2_sum/(moldyn->total_steps+1-moldyn->avg_skip); + moldyn->dv2_avg=moldyn->v2_avg-(moldyn->v_avg*moldyn->v_avg); + + /* temperature */ + moldyn->t_sum+=moldyn->t; + moldyn->t_avg=moldyn->t_sum/(moldyn->total_steps+1-moldyn->avg_skip); + + /* virial */ + + + /* pressure */ + + + return 0; +} + +int get_heat_capacity(t_moldyn *moldyn) { + + double temp2,ighc; + + /* averages needed for heat capacity calc */ + if(moldyn->total_stepsavg_skip) + return 0; + + /* (temperature average)^2 */ + temp2=moldyn->t_avg*moldyn->t_avg; + printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n", + moldyn->t_avg); + + /* ideal gas contribution */ + ighc=3.0*moldyn->count*K_BOLTZMANN/2.0; + printf(" ideal gas contribution: %f\n", + ighc/moldyn->mass*KILOGRAM/JOULE); + + /* specific heat for nvt ensemble */ + moldyn->c_v_nvt=moldyn->dv2_avg/(K_BOLTZMANN*temp2)+ighc; + moldyn->c_v_nvt/=moldyn->mass; + + /* specific heat for nve ensemble */ + moldyn->c_v_nve=ighc/(1.0-(moldyn->dv2_avg/(ighc*K_BOLTZMANN*temp2))); + moldyn->c_v_nve/=moldyn->mass; + + printf(" NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE); + printf(" NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE); +printf(" --> sim: %f experimental: %f\n",moldyn->dv2_avg,1.5*moldyn->count*K_B2*moldyn->t_avg*moldyn->t_avg*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM))); + + return 0; +} double thermodynamic_pressure_calc(t_moldyn *moldyn) { @@ -1304,8 +1403,7 @@ return 0; e_kin_calc(moldyn); temperature_calc(moldyn); pressure_calc(moldyn); - //tp=thermodynamic_pressure_calc(moldyn); -//printf("thermodynamic p: %f\n",thermodynamic_pressure_calc(moldyn)/BAR); + average_and_fluctuation_calc(moldyn); /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) @@ -1335,15 +1433,15 @@ return 0; if(!(i%p)) { dprintf(moldyn->pfd, "%f %f %f %f %f\n",moldyn->time, - moldyn->p/BAR,moldyn->mean_p/BAR, - moldyn->gp/BAR,moldyn->mean_gp/BAR); + moldyn->p/BAR,moldyn->p_avg/BAR, + moldyn->gp/BAR,moldyn->gp_avg/BAR); } } if(t) { if(!(i%t)) { dprintf(moldyn->tfd, "%f %f %f\n", - moldyn->time,moldyn->t,moldyn->mean_t); + moldyn->time,moldyn->t,moldyn->t_avg); } } if(s) { @@ -1371,9 +1469,10 @@ return 0; if(!(i%10)) { printf("\rsched: %d, steps: %d, T: %f, P: %f %f V: %f", sched->count,i, - moldyn->mean_t, - moldyn->mean_p/BAR, - moldyn->mean_gp/BAR, + moldyn->t_avg, + //moldyn->p_avg/BAR, + moldyn->p/BAR, + moldyn->gp_avg/BAR, moldyn->volume); fflush(stdout); } @@ -1385,11 +1484,13 @@ return 0; } /* check for hooks */ - if(sched->hook) - sched->hook(moldyn,sched->hook_params); - - /* get a new info line */ - printf("\n"); + if(sched->count+1total_sched) + if(sched->hook) { + printf("\n ## schedule hook %d/%d start ##\n", + sched->count+1,sched->total_sched); + sched->hook(moldyn,sched->hook_params); + printf(" ## schedule hook end ##\n"); + } } @@ -1805,3 +1906,26 @@ int moldyn_bc_check(t_moldyn *moldyn) { return 0; } + +/* + * post processing functions + */ + +int get_line(int fd,char *line,int max) { + + int count,ret; + + count=0; + + while(1) { + if(count==max) return count; + ret=read(fd,line+count,1); + if(ret<=0) return ret; + if(line[count]=='\n') { + line[count]='\0'; + return count+1; + } + count+=1; + } +} +