From: hackbard Date: Tue, 23 Jan 2007 16:36:18 +0000 (+0000) Subject: added report system, still messing around with virial pressure X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f3e4193447ac49a8953515910d4b4e6ce2c7608b;p=physik%2Fposic.git added report system, still messing around with virial pressure --- diff --git a/clean b/clean index 1f359db..e702e9e 100755 --- a/clean +++ b/clean @@ -1,2 +1,2 @@ -rm -rf $1 +rm -rf $1/* mkdir -p $1 diff --git a/moldyn.c b/moldyn.c index 94d596e..817a727 100644 --- a/moldyn.c +++ b/moldyn.c @@ -121,7 +121,7 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) { moldyn->vis.dim.z=z; } - moldyn->dv=0.0001*moldyn->volume; + moldyn->dv=0.000001*moldyn->volume; printf("[moldyn] dimensions in A and A^3 respectively:\n"); printf(" x: %f\n",moldyn->dim.x); @@ -199,6 +199,14 @@ int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) { return 0; } + +int moldyn_set_report(t_moldyn *moldyn,char *author,char *title) { + + strncpy(moldyn->rauthor,author,63); + strncpy(moldyn->rtitle,title,63); + + return 0; +} int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { @@ -247,6 +255,28 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { } printf("visual file (%d)\n",timer); break; + case CREATE_REPORT: + snprintf(filename,127,"%s/report.tex",moldyn->vlsdir); + moldyn->rfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->rfd<0) { + perror("[moldyn] report fd open"); + return moldyn->rfd; + } + snprintf(filename,127,"%s/plot.scr",moldyn->vlsdir); + moldyn->pfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->pfd<0) { + perror("[moldyn] plot fd open"); + return moldyn->pfd; + } + dprintf(moldyn->rfd,report_start, + moldyn->rauthor,moldyn->rtitle); + dprintf(moldyn->pfd,plot_script); + close(moldyn->pfd); + break; default: printf("unknown log type: %02x\n",type); return -1; @@ -257,9 +287,23 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { int moldyn_log_shutdown(t_moldyn *moldyn) { + char sc[256]; + printf("[moldyn] log shutdown\n"); if(moldyn->efd) close(moldyn->efd); if(moldyn->mfd) close(moldyn->mfd); + if(moldyn->rfd) { + dprintf(moldyn->rfd,report_end); + close(moldyn->rfd); + snprintf(sc,255,"cd %s && gnuplot plot.scr",moldyn->vlsdir); + system(sc); + snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir); + system(sc); + snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir); + system(sc); + snprintf(sc,255,"cd %s && dvipdf report",moldyn->vlsdir); + system(sc); + } if(&(moldyn->vis)) visual_tini(&(moldyn->vis)); return 0; @@ -299,10 +343,10 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, switch(type) { case CUBIC: - ret=cubic_init(a,b,c,atom,&origin); + ret=cubic_init(a,b,c,lc,atom,NULL); break; case FCC: - ret=fcc_init(a,b,c,lc,atom,&origin); + ret=fcc_init(a,b,c,lc,atom,NULL); break; case DIAMOND: ret=diamond_init(a,b,c,lc,atom,&origin); @@ -340,9 +384,38 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) { int count; - t_3dvec r + t_3dvec r; + int i,j,k; + t_3dvec o; + + count=0; + if(origin) + v3_copy(&o,origin); + else + v3_zero(&o); -HIER WEITER ! + r.x=o.x; + for(i=0;icount*moldyn->t*K_BOLTZMANN/moldyn->volume; -printf("temp = %f => ideal gas law pressure: %f\n",moldyn->t,p/ATM); return p; } @@ -607,8 +679,7 @@ printf("temp = %f => ideal gas law pressure: %f\n",moldyn->t,p/ATM); double pressure_calc(t_moldyn *moldyn) { int i; - double p1,p; - double v,h; + double v; t_virial *virial; v=0.0; @@ -616,18 +687,11 @@ double pressure_calc(t_moldyn *moldyn) { virial=&(moldyn->atom[i].virial); v+=(virial->xx+virial->yy+virial->zz); } + v*=ONE_THIRD; +printf("kieck mal: %f %f %f\n",v,moldyn->count*K_BOLTZMANN*moldyn->t,v/moldyn->count); - h=moldyn->count*K_BOLTZMANN*moldyn->t; - - p=(h-ONE_THIRD*v); - p/=moldyn->volume; - - p1=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt1); - p1/=moldyn->volume; - - printf("debug: vt1=%f v=%f nkt=%f\n",moldyn->vt1,v,h); - - printf("compare pressures: %f %f %f\n",p1/ATM,p/ATM,h/moldyn->volume/ATM); + moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+v; + moldyn->p/=moldyn->volume; return moldyn->p; } @@ -655,6 +719,8 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { scale=1.0+moldyn->dv/(moldyn->dim.y*moldyn->dim.z); scale_dim(moldyn,scale,TRUE,0,0); scale_atoms(moldyn,scale,TRUE,0,0); + link_cell_shutdown(moldyn); + link_cell_init(moldyn); potential_force_calc(moldyn); tp->x=(moldyn->energy-u)/moldyn->dv; p=tp->x*tp->x; @@ -667,6 +733,8 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.z); scale_dim(moldyn,scale,0,TRUE,0); scale_atoms(moldyn,scale,0,TRUE,0); + link_cell_shutdown(moldyn); + link_cell_init(moldyn); potential_force_calc(moldyn); tp->y=(moldyn->energy-u)/moldyn->dv; p+=tp->y*tp->y; @@ -679,6 +747,8 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.y); scale_dim(moldyn,scale,0,0,TRUE); scale_atoms(moldyn,scale,0,0,TRUE); + link_cell_shutdown(moldyn); + link_cell_init(moldyn); potential_force_calc(moldyn); tp->z=(moldyn->energy-u)/moldyn->dv; p+=tp->z*tp->z; @@ -687,15 +757,19 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; - printf("dU/dV komp addiert = %f\n",(tp->x+tp->y+tp->z)/ATM); + printf("dU/dV komp addiert = %f %f %f\n",tp->x,tp->y,tp->z); scale=1.0+pow(moldyn->dv/moldyn->volume,ONE_THIRD); +printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x); scale_dim(moldyn,scale,1,1,1); - scale_dim(moldyn,scale,1,1,1); + scale_atoms(moldyn,scale,1,1,1); + link_cell_shutdown(moldyn); + link_cell_init(moldyn); potential_force_calc(moldyn); +printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x); - printf("dU/dV einfach = %f\n",(moldyn->energy-u)/moldyn->dv/ATM); + printf("dU/dV einfach = %f\n",((moldyn->energy-u)/moldyn->dv)/ATM); /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); @@ -704,6 +778,9 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { /* restore energy */ moldyn->energy=u; + link_cell_shutdown(moldyn); + link_cell_init(moldyn); + return sqrt(p); } @@ -1054,6 +1131,7 @@ int moldyn_integrate(t_moldyn *moldyn) { int fd; char dir[128]; double ds; + double energy_scale; sched=&(moldyn->schedule); atom=moldyn->atom; @@ -1071,6 +1149,9 @@ int moldyn_integrate(t_moldyn *moldyn) { moldyn->tau_square=moldyn->tau*moldyn->tau; moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff; + /* energy scaling factor */ + energy_scale=moldyn->count*EV; + /* calculate initial forces */ potential_force_calc(moldyn); @@ -1115,7 +1196,9 @@ int moldyn_integrate(t_moldyn *moldyn) { if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT)) scale_volume(moldyn); -printf("-> %f\n",thermodynamic_pressure_calc(moldyn)/ATM); + temperature_calc(moldyn); + pressure_calc(moldyn); + //thermodynamic_pressure_calc(moldyn); /* check for log & visualization */ if(e) { @@ -1123,9 +1206,9 @@ printf("-> %f\n",thermodynamic_pressure_calc(moldyn)/ATM); update_e_kin(moldyn); dprintf(moldyn->efd, "%f %f %f %f\n", - moldyn->time,moldyn->ekin, - moldyn->energy, - get_total_energy(moldyn)); + moldyn->time,moldyn->ekin/energy_scale, + moldyn->energy/energy_scale, + get_total_energy(moldyn)/energy_scale); } if(m) { if(!(i%m)) { @@ -1152,8 +1235,8 @@ printf("-> %f\n",thermodynamic_pressure_calc(moldyn)/ATM); if(!(i%v)) { visual_atoms(&(moldyn->vis),moldyn->time, moldyn->atom,moldyn->count); - printf("\rsched: %d, steps: %d, debug: %d", - sched->count,i,moldyn->debug); + printf("\rsched: %d, steps: %d, debug: %f | %f", + sched->count,i,moldyn->p/ATM,moldyn->debug/ATM); fflush(stdout); } } @@ -1246,9 +1329,6 @@ int potential_force_calc(t_moldyn *moldyn) { /* reset energy */ moldyn->energy=0.0; - /* reset virial */ - moldyn->vt1=0.0; - /* reset force, site energy and virial of every atom */ for(i=0;ivirial.xx-=f->x*d->x; - a->virial.yy-=f->y*d->y; - a->virial.zz-=f->z*d->z; - a->virial.xy-=f->x*d->y; - a->virial.xz-=f->x*d->z; - a->virial.yz-=f->y*d->z; + a->virial.xx+=f->x*d->x; + a->virial.yy+=f->y*d->y; + a->virial.zz+=f->z*d->z; + a->virial.xy+=f->x*d->y; + a->virial.xz+=f->x*d->z; + a->virial.yz+=f->y*d->z; return 0; } @@ -1507,7 +1584,6 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { v3_add(&(ai->f),&(ai->f),&force); virial_calc(ai,&force,&distance); virial_calc(aj,&force,&distance); /* f and d signe switched */ - moldyn->vt1-=v3_scalar_product(&force,&distance); } return 0; diff --git a/moldyn.h b/moldyn.h index c910d7f..7329099 100644 --- a/moldyn.h +++ b/moldyn.h @@ -12,6 +12,7 @@ #include "random/random.h" #include "list/list.h" +#include "report/report.h" /* * @@ -82,7 +83,6 @@ typedef struct s_moldyn { t_3dvec dim; /* dimensions of the simulation volume */ double volume; /* volume of sim cell (dim.x*dim.y*dim.z) */ - double vt1,vt2; /* potential force function and parameter pointers */ int (*func1b)(struct s_moldyn *moldyn,t_atom *ai); @@ -139,6 +139,10 @@ typedef struct s_moldyn { int mfd; /* fd for momentum log */ unsigned int vwrite; /* how often to visualize atom information */ unsigned int swrite; /* how often to create a save file */ + int rfd; /* report file descriptor */ + char rtitle[64]; /* report title */ + char rauthor[64]; /* report author */ + int pfd; /* gnuplot script file descriptor */ u8 status; /* general moldyn properties */ @@ -162,7 +166,6 @@ typedef struct s_moldyn { #define P_SCALE_BERENDSEN 0x04 /* berendsen p control */ #define P_SCALE_DIRECT 0x08 /* direct p control */ - /* * * potential parameter structures @@ -307,6 +310,7 @@ typedef struct s_tersoff_mult_params { #define LOG_TOTAL_MOMENTUM 0x02 #define SAVE_STEP 0x04 #define VISUAL_STEP 0x08 +#define CREATE_REPORT 0x10 #define TRUE 1 #define FALSE 0 @@ -361,8 +365,9 @@ typedef struct s_tersoff_mult_params { * lattice constants */ -#define FCC 0x01 -#define DIAMOND 0x02 +#define CUBIC 0x01 +#define FCC 0x02 +#define DIAMOND 0x04 /* @@ -393,11 +398,13 @@ int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params); int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params); int moldyn_set_log_dir(t_moldyn *moldyn,char *dir); +int moldyn_set_report(t_moldyn *moldyn,char *author,char *title); int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer); int moldyn_log_shutdown(t_moldyn *moldyn); int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass, u8 attr,u8 brand,int a,int b,int c); +int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin); int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin); int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin); int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, diff --git a/report/report.h b/report/report.h new file mode 100644 index 0000000..7fd3180 --- /dev/null +++ b/report/report.h @@ -0,0 +1,64 @@ +/* + * report.h - header file of the report subsystem + * + * author: Frank Zirkelbach + * + */ + +#ifndef REPORT_H +#define REPORT_H + +static char report_start[]="\ +\\pdfoutput=0\n\ +\\documentclass[a4paper,11pt]{report}\n\ +\\usepackage[activate]{pdfcprot}\n\ +\\usepackage{verbatim}\n\ +\\usepackage{a4}\n\ +\\usepackage{a4wide}\n\ +\\usepackage[english,german]{babel}\n\ +\\usepackage[latin1]{inputenc}\n\ +\\usepackage[T1]{fontenc}\n\ +\\usepackage{amsmath}\n\ +\\usepackage{ae}\n\ +\\usepackage{aecompl}\n\ +\\usepackage[dvips]{graphicx}\n\ +\\graphicspath{{./}}\n\ +\\usepackage{color}\n\ +\\usepackage{pstricks}\n\ +\\usepackage{pst-node}\n\ +\\usepackage{rotating}\n\ +\\selectlanguage{english}\n\ +\\author{%s}\n\ +\\title{%s}\n\ +\\begin{document}\n\ +\\maketitle\n\ +"; + +static char report_end[]="\ +\\begin{figure}[!h]\n\ +\\begin{center}\n\ +\\includegraphics[width=10cm]{energy.eps}\n\ +\\caption{Kinetic, potential and total energy over time}\n\ +\\end{center}\n\ +\\end{figure}\n\ +\\end{document}\n\ +"; + +static char plot_script[]="\ +set autoscale \n\ +unset log \n\ +unset label \n\ +set xtic auto \n\ +set ytic auto \n\ +set title 'Energy vs. time' \n\ +set xlabel 'Time [fs]' \n\ +set ylabel 'Energy [eV]' \n\ +plot \"energy\" using 1:2 title 'Kinetic energy' with lines , \"energy\" using 1:3 title 'Potential energy' with lines , \"energy\" using 1:4 title 'Total energy' with lines \n\ +#set size 1.0, 0.6 \n\ +set terminal postscript eps enhanced color dashed lw 1 'Helvetica' 14 \n\ +set output 'energy.eps' \n\ +replot\ +"; + + +#endif diff --git a/sic.c b/sic.c index 5cfbb03..b5d4596 100644 --- a/sic.c +++ b/sic.c @@ -24,10 +24,8 @@ int hook(void *moldyn,void *hook_params) { if(md->schedule.count==1) set_pt_scale(md,0,0,0,0); - //set_temperature(md,md->t_ref-100.0); - return 0; } @@ -66,7 +64,7 @@ int main(int argc,char **argv) { /* cutoff radius */ //set_cutoff(&md,TM_S_SI); - set_cutoff(&md,2*LC_SI); + set_cutoff(&md,2*LC_SI*0.5*sqrt(1.5)); /* * potential parameters @@ -115,16 +113,17 @@ int main(int argc,char **argv) { tersoff_mult_complete_params(&tp); /* set (initial) dimensions of simulation volume */ - set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE); + set_dim(&md,8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),TRUE); /* set periodic boundary conditions in all directions */ set_pbc(&md,TRUE,TRUE,TRUE); /* create the lattice / place atoms */ - create_lattice(&md,DIAMOND,LC_SI,SI,M_SI, + create_lattice(&md,FCC,LC_SI*0.5*sqrt(1.5),SI,M_SI, + //create_lattice(&md,DIAMOND,LC_SI,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); + 0,8,8,8); moldyn_bc_check(&md); /* testing configuration */ @@ -161,24 +160,26 @@ int main(int argc,char **argv) { /* set p/t scaling */ //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0, // T_SCALE_BERENDSEN,100.0); - set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0); + //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0); //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,0,0); /* initial thermal fluctuations of particles (in equilibrium) */ thermal_init(&md,TRUE); /* create the simulation schedule */ - moldyn_add_schedule(&md,1001,3.0); - moldyn_add_schedule(&md,501,1.0); - moldyn_add_schedule(&md,501,1.0); + moldyn_add_schedule(&md,1001,1.0); + //moldyn_add_schedule(&md,501,1.0); + //moldyn_add_schedule(&md,501,1.0); /* schedule hook function */ - moldyn_set_schedule_hook(&md,&hook,NULL); + //moldyn_set_schedule_hook(&md,&hook,NULL); /* activate logging */ moldyn_set_log_dir(&md,argv[1]); - moldyn_set_log(&md,LOG_TOTAL_ENERGY,5); - moldyn_set_log(&md,VISUAL_STEP,50); + moldyn_set_report(&md,"Frank Zirkelbach","Test 1"); + moldyn_set_log(&md,LOG_TOTAL_ENERGY,1); + moldyn_set_log(&md,VISUAL_STEP,10); + moldyn_set_log(&md,CREATE_REPORT,0); /* * let's do the actual md algorithm now diff --git a/visual/visual.c b/visual/visual.c index a6fda4f..a76c608 100644 --- a/visual/visual.c +++ b/visual/visual.c @@ -87,7 +87,7 @@ 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 270\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");