From 58fd691b276fbe87593036714f26dbfe7486cbeb Mon Sep 17 00:00:00 2001 From: hackbard Date: Wed, 27 Jun 2007 22:24:34 +0000 Subject: [PATCH] added post proc stuff --- moldyn.c | 119 +++++++++++++++++++++++++++++++++++-------------------- moldyn.h | 20 +++++++++- ppm2avi | 2 +- sic.c | 34 ++++++++++++---- 4 files changed, 123 insertions(+), 52 deletions(-) diff --git a/moldyn.c b/moldyn.c index aebd604..c130ef4 100644 --- a/moldyn.c +++ b/moldyn.c @@ -404,11 +404,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)); @@ -1385,8 +1388,9 @@ return 0; } /* check for hooks */ - if(sched->hook) - sched->hook(moldyn,sched->hook_params); + if(sched->count+1total_sched) + if(sched->hook) + sched->hook(moldyn,sched->hook_params); /* get a new info line */ printf("\n"); @@ -1809,17 +1813,17 @@ int moldyn_bc_check(t_moldyn *moldyn) { /* * postprocessing functions */ -#define LINE_MAX 128 -int read_line(int fd,char *line) { + +int get_line(int fd,char *line,int max) { int count,ret; count=0; while(1) { - if(count==LINE_MAX) return count; + if(count==max) return count; ret=read(fd,line+count,1); - if(ret<0) return ret; + if(ret<=0) return ret; if(line[count]=='\n') { line[count]='\0'; return count+1; @@ -1828,19 +1832,21 @@ int read_line(int fd,char *line) { } } -int calc_fluctuations(double start,double end,char *file) { +int calc_fluctuations(double start,double end,t_moldyn *moldyn) { int fd; int count,ret; double time,pot,kin,tot; - double p_m,k_m,t_m; - double p2_m,k2_m,t2_m; double p_sum,k_sum,t_sum; - char buf[LINE_MAX]; + char buf[64]; + char file[128+7]; + + printf("[moldyn] calculating energy fluctuations [eV]:\n"); + snprintf(file,128+7,"%s/energy",moldyn->vlsdir); fd=open(file,O_RDONLY); if(fd<0) { - perror("[moldyn] post proc open"); + perror("[moldyn] post proc energy open"); return fd; } @@ -1850,24 +1856,21 @@ int calc_fluctuations(double start,double end,char *file) { t_sum=0.0; count=0; while(1) { - ret=read_line(fd,buf); - if(ret<0) break; -printf("%d\n",ret); + ret=get_line(fd,buf,63); + if(ret<=0) break; if(buf[0]=='#') continue; sscanf(buf,"%lf %lf %lf %lf",&time,&kin,&pot,&tot); -printf("%f %f %f %f\n",time,pot,kin,tot); - //if(time>end) break; - if((time>=start)&(time<=end)) { - p_sum+=pot; - k_sum+=kin; - t_sum+=tot; - count+=1; - } + if(timeend) break; + p_sum+=pot; + k_sum+=kin; + t_sum+=tot; + count+=1; } - p_m=p_sum/count; - k_m=k_sum/count; - t_m=t_sum/count; + moldyn->p_m=p_sum/count; + moldyn->k_m=k_sum/count; + moldyn->t_m=t_sum/count; /* mean square fluctuations */ if(lseek(fd,SEEK_SET,0)<0) { @@ -1879,28 +1882,60 @@ printf("%f %f %f %f\n",time,pot,kin,tot); k_sum=0.0; t_sum=0.0; while(1) { - ret=read_line(fd,buf); - if(ret<0) break; + ret=get_line(fd,buf,63); + if(ret<=0) break; if(buf[0]=='#') continue; sscanf(buf,"%lf %lf %lf %lf",&time,&kin,&pot,&tot); + if(timeend) break; - if((time>=start)&(time<=end)) { - p_sum+=((pot-p_m)*(pot-p_m)); - k_sum+=((kin-k_m)*(kin-k_m)); - t_sum+=((tot-t_m)*(tot-t_m)); - count+=1; - } + k_sum+=((kin-moldyn->k_m)*(kin-moldyn->k_m)); + p_sum+=((pot-moldyn->p_m)*(pot-moldyn->p_m)); + t_sum+=((tot-moldyn->t_m)*(tot-moldyn->t_m)); + count+=1; } - p2_m=p_sum/count; - k2_m=k_sum/count; - t2_m=t_sum/count; + moldyn->dp2_m=p_sum/count; + moldyn->dk2_m=k_sum/count; + moldyn->dt2_m=t_sum/count; - printf("[moldyn] fluctuations (%f - %f):\n",start,end); - printf(" - averages : %f %f %f\n",k_m,p_m,t_m); - printf(" - mean square: %f %f %f\n",k2_m,p2_m,t2_m); + printf(" averages : %f %f %f\n",moldyn->k_m, + moldyn->p_m, + moldyn->t_m); + printf(" mean square: %f %f %f\n",moldyn->dk2_m, + moldyn->dp2_m, + moldyn->dt2_m); close(fd); return 0; } + +int get_heat_capacity(t_moldyn *moldyn) { + + double temp2,mass,ighc; + int i; + + /* (temperature average)^2 */ + temp2=2.0*moldyn->k_m*EV/(3.0*K_BOLTZMANN); + printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",temp2); + temp2*=temp2; + + /* total mass */ + mass=0.0; + for(i=0;icount;i++) + mass+=moldyn->atom[i].mass; + + /* ideal gas contribution */ + ighc=3.0*moldyn->count*K_BOLTZMANN/2.0; + printf(" ideal gas contribution: %f\n",ighc/mass*KILOGRAM/JOULE); + + moldyn->c_v_nvt=moldyn->dp2_m*moldyn->count*moldyn->count*EV/(K_BOLTZMANN*temp2)+ighc; + moldyn->c_v_nvt/=mass; + moldyn->c_v_nve=ighc/(1.0-(moldyn->dp2_m*moldyn->count*moldyn->count*EV/(ighc*K_BOLTZMANN*temp2))); + moldyn->c_v_nve/=mass; + + printf(" NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE); + printf(" NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE); + + return 0; +} diff --git a/moldyn.h b/moldyn.h index fc9bce7..484c456 100644 --- a/moldyn.h +++ b/moldyn.h @@ -139,9 +139,22 @@ typedef struct s_moldyn { double tau_square; /* delta t squared */ int total_steps; /* total steps */ + /* energy */ double energy; /* potential energy */ double ekin; /* kinetic energy */ + /* energy averages & fluctuations */ + double k_m; + double p_m; + double t_m; + double dk2_m; /* mean square kinetic energy fluctuations */ + double dp2_m; /* mean square potential energy fluctuations */ + double dt2_m; /* mean square fluctuations in total energy */ + + /* response functions */ + double c_v_nve; /* constant volume heat capacity (nve) */ + double c_v_nvt; /* constant volume heat capacity (nvt) */ + char vlsdir[128]; /* visualization/log/save directory */ t_visual vis; /* visualization interface structure */ u8 vlsprop; /* log/vis/save properties */ @@ -209,7 +222,9 @@ typedef struct s_moldyn { #define PASCAL (NEWTON/(METER*METER)) /* N / A^2 */ #define BAR ((1.0e5*PASCAL)) /* N / A^2 */ #define K_BOLTZMANN (1.380650524e-23*METER*NEWTON) /* NA/K */ +#define K_B2 (K_BOLTZMANN*K_BOLTZMANN) /* (NA)^2/K^2 */ #define EV (1.6021765314e-19*METER*NEWTON) /* NA */ +#define JOULE (NEWTON*METER) /* NA */ #define MOLDYN_TEMP 273.0 #define MOLDYN_TAU 1.0 @@ -427,7 +442,8 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a); int moldyn_bc_check(t_moldyn *moldyn); -int read_line(int fd,char *line); -int calc_fluctuations(double start,double end,char *file); +int get_line(int fd,char *line,int max); +int calc_fluctuations(double start,double end,t_moldyn *moldyn); +int get_heat_capacity(t_moldyn *moldyn); #endif diff --git a/ppm2avi b/ppm2avi index 74e65cf..fec403e 100755 --- a/ppm2avi +++ b/ppm2avi @@ -2,4 +2,4 @@ for i in $1/*.ppm; do convert $i $1/$(basename $i ppm)png done -mencoder mf://$1/*.png -ovc copy -o $1/md.avi +mencoder mf://$1/*.png -ovc copy -o $1/md.avi >/dev/null 2>&1 diff --git a/sic.c b/sic.c index cea2644..d50c35e 100644 --- a/sic.c +++ b/sic.c @@ -35,6 +35,12 @@ int hook(void *moldyn,void *hook_params) { md=moldyn; +// vortrag +set_temperature(moldyn,(4-md->schedule.count)*1000.0); +set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0); + +return 0; + printf("\nschedule hook: "); if(!(md->schedule.count%2)) { @@ -203,11 +209,11 @@ int main(int argc,char **argv) { albe_mult_complete_params(&ap); /* set (initial) dimensions of simulation volume */ - //set_dim(&md,6*LC_SI_ALBE,6*LC_SI_ALBE,6*LC_SI_ALBE,TRUE); + set_dim(&md,6*LC_SI_ALBE,6*LC_SI_ALBE,6*LC_SI_ALBE,TRUE); //set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE); //set_dim(&md,6*LC_C_ALBE,6*LC_C_ALBE,6*LC_C_ALBE,TRUE); //set_dim(&md,6*LC_C,6*LC_C,6*LC_C,TRUE); - set_dim(&md,6*LC_SIC_ALBE,6*LC_SIC_ALBE,6*LC_SIC_ALBE,TRUE); + //set_dim(&md,6*LC_SIC_ALBE,6*LC_SIC_ALBE,6*LC_SIC_ALBE,TRUE); /* set periodic boundary conditions in all directions */ set_pbc(&md,TRUE,TRUE,TRUE); @@ -215,14 +221,14 @@ int main(int argc,char **argv) { /* create the lattice / place atoms */ //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI, //create_lattice(&md,DIAMOND,LC_C_ALBE,C,M_C, - //create_lattice(&md,DIAMOND,LC_SI_ALBE,SI,M_SI, - // ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + create_lattice(&md,DIAMOND,LC_SI_ALBE,SI,M_SI, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, // ATOM_ATTR_2BP|ATOM_ATTR_HB, - // 0,6,6,6,NULL); + 0,6,6,6,NULL); // 1,6,6,6,NULL); /* create centered zinc blende lattice */ - /**/ + /* r.x=0.5*0.25*LC_SIC_ALBE; r.y=r.x; r.z=r.x; create_lattice(&md,FCC,LC_SIC_ALBE,SI,M_SI, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, @@ -231,7 +237,7 @@ int main(int argc,char **argv) { create_lattice(&md,FCC,LC_SIC_ALBE,C,M_C, ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, 1,6,6,6,&r); - /**/ + */ moldyn_bc_check(&md); @@ -283,6 +289,11 @@ int main(int argc,char **argv) { /* create the simulation schedule */ /* initial configuration */ moldyn_add_schedule(&md,10000,1.0); + //moldyn_add_schedule(&md,1000,1.0); + //moldyn_add_schedule(&md,1000,1.0); + //moldyn_add_schedule(&md,1000,1.0); + //moldyn_add_schedule(&md,1000,1.0); + //moldyn_add_schedule(&md,1000,1.0); /* adding atoms */ //for(inject=0;inject