X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=f564075adece43f086ccd660cba312c4e6cd1ed3;hb=4871747c5c848e5881bea7949a41ceb589263841;hp=94d596e91878237cf41a127cd8aa947cbcd18997;hpb=fb951c04e522e4637618bf622fc67194c2a7b15f;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 94d596e..f564075 100644 --- a/moldyn.c +++ b/moldyn.c @@ -16,6 +16,7 @@ #include #include "moldyn.h" +#include "report/report.h" int moldyn_init(t_moldyn *moldyn,int argc,char **argv) { @@ -121,7 +122,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 +200,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 +256,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 +288,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 +344,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 +385,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); + + r.x=o.x; + for(i=0;iatom; - - double_ekin=0; - for(i=0;icount;i++) - double_ekin+=atom[i].mass*v3_absolute_square(&(atom[i].v)); + /* assume up to date kinetic energy, which is 3/2 N k_B T */ - /* kinetic energy = 3/2 N k_B T */ - moldyn->t=double_ekin/(3.0*K_BOLTZMANN*moldyn->count); + moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count); return moldyn->t; } @@ -599,7 +664,6 @@ double ideal_gas_law_pressure(t_moldyn *moldyn) { double p; p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume; -printf("temp = %f => ideal gas law pressure: %f\n",moldyn->t,p/ATM); return p; } @@ -607,27 +671,24 @@ 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; + /* + * P = 1/(3V) sum_i ( p_i^2 / 2m + f_i r_i ) + * + * virial = f_i r_i + */ + v=0.0; for(i=0;icount;i++) { virial=&(moldyn->atom[i].virial); v+=(virial->xx+virial->yy+virial->zz); } - 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); + /* assume up to date kinetic energy */ + moldyn->p=2.0*moldyn->ekin+v; + moldyn->p/=(3.0*moldyn->volume); return moldyn->p; } @@ -655,6 +716,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 +730,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 +744,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 +754,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 +775,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 +1128,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 +1146,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,17 +1193,19 @@ 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); + update_e_kin(moldyn); + temperature_calc(moldyn); + pressure_calc(moldyn); + //thermodynamic_pressure_calc(moldyn); /* check for log & visualization */ if(e) { if(!(i%e)) - 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 +1232,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->p/ATM); fflush(stdout); } } @@ -1246,9 +1326,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; } @@ -1503,11 +1577,10 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { d*=eps; v3_scale(&force,&distance,d); v3_add(&(aj->f),&(aj->f),&force); - v3_scale(&force,&distance,-1.0*d); /* f = - grad E */ + v3_scale(&force,&force,-1.0); /* f = - grad E */ 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;