From ea612b88a0588b8f46fafaebf3b37fb46c83c0cf Mon Sep 17 00:00:00 2001 From: hackbard Date: Wed, 14 Mar 2007 15:12:05 +0000 Subject: [PATCH] improved log/report subsystem, playing around w/ pressure, sic hook testing --- moldyn.c | 192 ++++++++++++++++++++++++++++++++----------- moldyn.h | 16 +++- potentials/tersoff.c | 39 +++++---- report/report.h | 55 ++++++++++++- sic.c | 65 ++++++++++++--- 5 files changed, 282 insertions(+), 85 deletions(-) diff --git a/moldyn.c b/moldyn.c index ca6f84a..9a62b9f 100644 --- a/moldyn.c +++ b/moldyn.c @@ -246,6 +246,32 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { dprintf(moldyn->efd,"# total momentum log file\n"); printf("total momentum (%d)\n",timer); break; + case LOG_PRESSURE: + moldyn->pwrite=timer; + snprintf(filename,127,"%s/pressure",moldyn->vlsdir); + moldyn->pfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->pfd<0) { + perror("[moldyn] pressure log file\n"); + return moldyn->pfd; + } + dprintf(moldyn->pfd,"# pressure log file\n"); + printf("pressure (%d)\n",timer); + break; + case LOG_TEMPERATURE: + moldyn->twrite=timer; + snprintf(filename,127,"%s/temperature",moldyn->vlsdir); + moldyn->tfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->tfd<0) { + perror("[moldyn] temperature log file\n"); + return moldyn->tfd; + } + dprintf(moldyn->tfd,"# temperature log file\n"); + printf("temperature (%d)\n",timer); + break; case SAVE_STEP: moldyn->swrite=timer; printf("save file (%d)\n",timer); @@ -268,18 +294,47 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) { 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; + if(moldyn->efd) { + snprintf(filename,127,"%s/e_plot.scr", + moldyn->vlsdir); + moldyn->epfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->epfd<0) { + perror("[moldyn] energy plot fd open"); + return moldyn->epfd; + } + dprintf(moldyn->epfd,e_plot_script); + close(moldyn->epfd); + } + if(moldyn->pfd) { + snprintf(filename,127,"%s/pressure_plot.scr", + moldyn->vlsdir); + moldyn->ppfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->ppfd<0) { + perror("[moldyn] p plot fd open"); + return moldyn->ppfd; + } + dprintf(moldyn->ppfd,pressure_plot_script); + close(moldyn->ppfd); + } + if(moldyn->tfd) { + snprintf(filename,127,"%s/temperature_plot.scr", + moldyn->vlsdir); + moldyn->tpfd=open(filename, + O_WRONLY|O_CREAT|O_EXCL, + S_IRUSR|S_IWUSR); + if(moldyn->tpfd<0) { + perror("[moldyn] t plot fd open"); + return moldyn->tpfd; + } + dprintf(moldyn->tpfd,temperature_plot_script); + close(moldyn->tpfd); } 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); @@ -294,13 +349,35 @@ int moldyn_log_shutdown(t_moldyn *moldyn) { char sc[256]; printf("[moldyn] log shutdown\n"); - if(moldyn->efd) close(moldyn->efd); + if(moldyn->efd) { + close(moldyn->efd); + if(moldyn->rfd) { + dprintf(moldyn->rfd,report_energy); + snprintf(sc,255,"cd %s && gnuplot e_plot.scr", + moldyn->vlsdir); + system(sc); + } + } if(moldyn->mfd) close(moldyn->mfd); + if(moldyn->pfd) { + close(moldyn->pfd); + if(moldyn->rfd) + dprintf(moldyn->rfd,report_pressure); + snprintf(sc,255,"cd %s && gnuplot pressure_plot.scr", + moldyn->vlsdir); + system(sc); + } + if(moldyn->tfd) { + close(moldyn->tfd); + if(moldyn->rfd) + dprintf(moldyn->rfd,report_temperature); + snprintf(sc,255,"cd %s && gnuplot temperature_plot.scr", + moldyn->vlsdir); + system(sc); + } 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); @@ -676,9 +753,11 @@ double pressure_calc(t_moldyn *moldyn) { t_virial *virial; /* - * P = 1/(3V) sum_i ( p_i^2 / 2m + f_i r_i ) - * - * virial = f_i r_i + * PV = NkT + + * W = 1/3 sum_i f_i r_i + * virial = sum_i f_i r_i + * + * => P = (2 Ekin + virial) / (3V) */ v=0.0; @@ -698,9 +777,20 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { t_3dvec dim,*tp; double u,p; - double scale; + double scale,dv; t_atom *store; + /* + * dU = - p dV + * + * => p = - dU/dV + * + * dV: dx,y,z = 0.001 x,y,z + */ + + scale=1.00001; +printf("\n\nP-DEBUG:\n"); + tp=&(moldyn->tp); store=malloc(moldyn->count*sizeof(t_atom)); if(store==NULL) { @@ -714,27 +804,28 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { dim=moldyn->dim; /* derivative with respect to x direction */ - 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); + dv=0.00001*moldyn->dim.x*moldyn->dim.y*moldyn->dim.z; link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - tp->x=(moldyn->energy-u)/moldyn->dv; + tp->x=(moldyn->energy-u)/dv; p=tp->x*tp->x; +printf("e: %f eV de: %f eV dV: %f A^3\n",moldyn->energy/moldyn->count/EV,(moldyn->energy-u)/moldyn->count/EV,dv); /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; /* derivative with respect to y direction */ - 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); + dv=0.00001*moldyn->dim.y*moldyn->dim.x*moldyn->dim.z; link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - tp->y=(moldyn->energy-u)/moldyn->dv; + tp->y=(moldyn->energy-u)/dv; p+=tp->y*tp->y; /* restore atomic configuration + dim */ @@ -742,37 +833,19 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { moldyn->dim=dim; /* derivative with respect to z direction */ - 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); + dv=0.00001*moldyn->dim.z*moldyn->dim.x*moldyn->dim.y; link_cell_shutdown(moldyn); link_cell_init(moldyn,QUIET); potential_force_calc(moldyn); - tp->z=(moldyn->energy-u)/moldyn->dv; + tp->z=(moldyn->energy-u)/dv; p+=tp->z*tp->z; /* restore atomic configuration + dim */ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); moldyn->dim=dim; - 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_atoms(moldyn,scale,1,1,1); - link_cell_shutdown(moldyn); - link_cell_init(moldyn,QUIET); - 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); - - /* restore atomic configuration + dim */ - memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom)); - moldyn->dim=dim; - /* restore energy */ moldyn->energy=u; @@ -1109,14 +1182,15 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) { int moldyn_integrate(t_moldyn *moldyn) { int i; - unsigned int e,m,s,v; - t_3dvec p; + unsigned int e,m,s,v,p,t; + t_3dvec momentum; t_moldyn_schedule *sched; t_atom *atom; int fd; char dir[128]; double ds; double energy_scale; + //double tp; sched=&(moldyn->schedule); atom=moldyn->atom; @@ -1129,6 +1203,8 @@ int moldyn_integrate(t_moldyn *moldyn) { m=moldyn->mwrite; s=moldyn->swrite; v=moldyn->vwrite; + p=moldyn->pwrite; + t=moldyn->twrite; /* sqaure of some variables */ moldyn->tau_square=moldyn->tau*moldyn->tau; @@ -1179,7 +1255,7 @@ int moldyn_integrate(t_moldyn *moldyn) { update_e_kin(moldyn); temperature_calc(moldyn); pressure_calc(moldyn); - //thermodynamic_pressure_calc(moldyn); + //tp=thermodynamic_pressure_calc(moldyn); /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) @@ -1198,9 +1274,23 @@ int moldyn_integrate(t_moldyn *moldyn) { } if(m) { if(!(i%m)) { - p=get_total_p(moldyn); + momentum=get_total_p(moldyn); dprintf(moldyn->mfd, - "%f %f\n",moldyn->time,v3_norm(&p)); + "%f %f %f %f %f\n",moldyn->time, + momentum.x,momentum.y,momentum.z, + v3_norm(&momentum)); + } + } + if(p) { + if(!(i%p)) { + dprintf(moldyn->pfd, + "%f %f\n",moldyn->time,moldyn->p/ATM); + } + } + if(t) { + if(!(i%t)) { + dprintf(moldyn->tfd, + "%f %f\n",moldyn->time,moldyn->t); } } if(s) { @@ -1221,13 +1311,17 @@ int moldyn_integrate(t_moldyn *moldyn) { if(!(i%v)) { visual_atoms(&(moldyn->vis),moldyn->time, moldyn->atom,moldyn->count); - printf("\rsched: %d, steps: %d, T: %f, P: %f V: %f", - sched->count,i, - moldyn->t,moldyn->p/ATM,moldyn->volume); - fflush(stdout); } } + /* display progress */ + if(!(i%10)) { + printf("\rsched: %d, steps: %d, T: %f, P: %f V: %f", + sched->count,i, + moldyn->t,moldyn->p/ATM,moldyn->volume); + fflush(stdout); + } + /* increase absolute time */ moldyn->time+=moldyn->tau; diff --git a/moldyn.h b/moldyn.h index b140372..85a4950 100644 --- a/moldyn.h +++ b/moldyn.h @@ -133,12 +133,18 @@ typedef struct s_moldyn { int efd; /* fd for energy log */ unsigned int mwrite; /* how often to log momentum */ int mfd; /* fd for momentum log */ + unsigned int pwrite; /* how often to log pressure */ + int pfd; /* fd for pressure log */ + unsigned int twrite; /* how often to log temperature */ + int tfd; /* fd for temperature 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 */ + int epfd; /* energy gnuplot script file descriptor */ + int ppfd; /* pressure gnuplot script file descriptor */ + int tpfd; /* temperature gnuplot script file descriptor */ u8 status; /* general moldyn properties */ @@ -201,9 +207,11 @@ typedef struct s_moldyn { #define LOG_TOTAL_ENERGY 0x01 #define LOG_TOTAL_MOMENTUM 0x02 -#define SAVE_STEP 0x04 -#define VISUAL_STEP 0x08 -#define CREATE_REPORT 0x10 +#define LOG_PRESSURE 0x04 +#define LOG_TEMPERATURE 0x08 +#define SAVE_STEP 0x10 +#define VISUAL_STEP 0x20 +#define CREATE_REPORT 0x40 #define TRUE 1 #define FALSE 0 diff --git a/potentials/tersoff.c b/potentials/tersoff.c index 01364ea..16f771c 100644 --- a/potentials/tersoff.c +++ b/potentials/tersoff.c @@ -211,12 +211,13 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { v3_add(&(ai->f),&(ai->f),&force); /* virial */ - ai->virial.xx-=force.x*dist_ij.x; - ai->virial.yy-=force.y*dist_ij.y; - ai->virial.zz-=force.z*dist_ij.z; - ai->virial.xy-=force.x*dist_ij.y; - ai->virial.xz-=force.x*dist_ij.z; - ai->virial.yz-=force.y*dist_ij.z; + virial_calc(ai,&force,&dist_ij); + //ai->virial.xx-=force.x*dist_ij.x; + //ai->virial.yy-=force.y*dist_ij.y; + //ai->virial.zz-=force.z*dist_ij.z; + //ai->virial.xy-=force.x*dist_ij.y; + //ai->virial.xz-=force.x*dist_ij.z; + //ai->virial.yz-=force.y*dist_ij.z; #ifdef DEBUG if(ai==&(moldyn->atom[0])) { @@ -326,12 +327,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { v3_add(&(ai->f),&(ai->f),&force); /* virial */ - ai->virial.xx-=force.x*dist_ij->x; - ai->virial.yy-=force.y*dist_ij->y; - ai->virial.zz-=force.z*dist_ij->z; - ai->virial.xy-=force.x*dist_ij->y; - ai->virial.xz-=force.x*dist_ij->z; - ai->virial.yz-=force.y*dist_ij->z; + virial_calc(ai,&force,dist_ij); + //ai->virial.xx-=force.x*dist_ij->x; + //ai->virial.yy-=force.y*dist_ij->y; + //ai->virial.zz-=force.z*dist_ij->z; + //ai->virial.xy-=force.x*dist_ij->y; + //ai->virial.xz-=force.x*dist_ij->z; + //ai->virial.yz-=force.y*dist_ij->z; #ifdef DEBUG if(ai==&(moldyn->atom[0])) { @@ -380,12 +382,13 @@ if(ai==&(moldyn->atom[0])) { /* virial - plus sign, as dist_ij = - dist_ji - (really??) */ // TEST ... with a minus instead - ai->virial.xx-=force.x*dist_ij->x; - ai->virial.yy-=force.y*dist_ij->y; - ai->virial.zz-=force.z*dist_ij->z; - ai->virial.xy-=force.x*dist_ij->y; - ai->virial.xz-=force.x*dist_ij->z; - ai->virial.yz-=force.y*dist_ij->z; + virial_calc(ai,&force,dist_ij); + //ai->virial.xx-=force.x*dist_ij->x; + //ai->virial.yy-=force.y*dist_ij->y; + //ai->virial.zz-=force.z*dist_ij->z; + //ai->virial.xy-=force.x*dist_ij->y; + //ai->virial.xz-=force.x*dist_ij->z; + //ai->virial.yz-=force.y*dist_ij->z; #ifdef DEBUG if(ai==&(moldyn->atom[0])) { diff --git a/report/report.h b/report/report.h index e46ff2a..4fe42f7 100644 --- a/report/report.h +++ b/report/report.h @@ -34,17 +34,38 @@ static char report_start[]="\ \\maketitle\n\ "; -static char report_end[]="\ +static char report_energy[]="\ \\begin{figure}[!h]\n\ \\begin{center}\n\ \\includegraphics[width=16cm]{energy.eps}\n\ \\caption{Kinetic, potential and total energy over time}\n\ \\end{center}\n\ \\end{figure}\n\ +"; + +static char report_pressure[]="\ +\\begin{figure}[!h]\n\ +\\begin{center}\n\ +\\includegraphics[width=16cm]{pressure.eps}\n\ +\\caption{Pressure over time}\n\ +\\end{center}\n\ +\\end{figure}\n\ +"; + +static char report_temperature[]="\ +\\begin{figure}[!h]\n\ +\\begin{center}\n\ +\\includegraphics[width=16cm]{temperature.eps}\n\ +\\caption{Temperature over time}\n\ +\\end{center}\n\ +\\end{figure}\n\ +"; + +static char report_end[]="\ \\end{document}\n\ "; -static char plot_script[]="\ +static char e_plot_script[]="\ set autoscale \n\ unset log \n\ unset label \n\ @@ -55,7 +76,35 @@ set xlabel 'Time [fs]' \n\ set ylabel 'Energy per atom [eV]' \n\ set terminal postscript eps enhanced color solid lw 1 'Helvetica' 14 \n\ set output 'energy.eps' \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\ +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 \ +"; + +static char pressure_plot_script[]="\ +set autoscale \n\ +unset log \n\ +unset label \n\ +set xtic auto \n\ +set ytic auto \n\ +set title 'Pressure vs. time' \n\ +set xlabel 'Time [fs]' \n\ +set ylabel 'Pressure [atm]' \n\ +set terminal postscript eps enhanced color solid lw 1 'Helvetica' 14 \n\ +set output 'pressure.eps' \n\ +plot \"pressure\" using 1:2 title 'Pressure' with lines \ +"; + +static char temperature_plot_script[]="\ +set autoscale \n\ +unset log \n\ +unset label \n\ +set xtic auto \n\ +set ytic auto \n\ +set title 'Temperature vs. time' \n\ +set xlabel 'Time [fs]' \n\ +set ylabel 'Temperature [K]' \n\ +set terminal postscript eps enhanced color solid lw 1 'Helvetica' 14 \n\ +set output 'temperature.eps' \n\ +plot \"temperature\" using 1:2 title 'Temperature' with lines \ "; #endif diff --git a/sic.c b/sic.c index 167735f..682d89e 100644 --- a/sic.c +++ b/sic.c @@ -15,20 +15,51 @@ #include "potentials/lennard_jones.h" #include "potentials/tersoff.h" +#define INJECT 20 +#define NR_ATOMS 20 + int hook(void *moldyn,void *hook_params) { t_moldyn *md; + t_3dvec r,v,dist; + double d; + unsigned char run; + int i,j; + t_atom *atom; md=moldyn; - /* switch to direct scaling in first hook */ - if(md->schedule.count==0) + printf("\nschedule hook: "); + + if(!(md->schedule.count%2)) { + /* add carbon at random place, and enable t scaling */ + for(j=0;jrandom))*md->dim.x; + r.y=rand_get_double(&(md->random))*md->dim.y; + r.z=rand_get_double(&(md->random))*md->dim.z; + for(i=0;icount;i++) { + atom=&(md->atom[i]); + v3_sub(&dist,&(atom->r),&r); + d=v3_absolute_square(&dist); + if(d>TM_R_C) + run=0; + } + } + v.x=0; v.y=0; v.z=0; + add_atom(md,C,M_C,1, + ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB, + &r,&v); + } + printf("adding atoms & enable t scaling\n"); set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0); - /* switch off temp scaling in second hook */ - if(md->schedule.count==1) + } + else { + /* disable t scaling */ + printf("disabling t scaling\n"); set_pt_scale(md,0,0,0,0); - - //set_temperature(md,md->t_ref-100.0); + } return 0; } @@ -49,6 +80,9 @@ int main(int argc,char **argv) { t_ho_params ho; t_tersoff_mult_params tp; + /* atom injection counter */ + int inject; + /* testing location & velocity vector */ t_3dvec r,v; memset(&r,0,sizeof(t_3dvec)); @@ -176,18 +210,27 @@ int main(int argc,char **argv) { thermal_init(&md,TRUE); /* create the simulation schedule */ - moldyn_add_schedule(&md,10001,1.0); - //moldyn_add_schedule(&md,501,1.0); - //moldyn_add_schedule(&md,501,1.0); + /* initial configuration */ + moldyn_add_schedule(&md,500,1.0); + /* adding atoms */ + for(inject=0;inject