X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=0d5027c57df699447bcdf61a2c2b516235c2fa0d;hb=91e4a00285261865ffce4b9553d153b562c80ee6;hp=94d596e91878237cf41a127cd8aa947cbcd18997;hpb=fb951c04e522e4637618bf622fc67194c2a7b15f;p=physik%2Fposic.git diff --git a/moldyn.c b/moldyn.c index 94d596e..0d5027c 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) { @@ -72,7 +73,7 @@ int set_temperature(t_moldyn *moldyn,double t_ref) { moldyn->t_ref=t_ref; - printf("[moldyn] temperature: %f\n",moldyn->t_ref); + printf("[moldyn] temperature [K]: %f\n",moldyn->t_ref); return 0; } @@ -81,7 +82,7 @@ int set_pressure(t_moldyn *moldyn,double p_ref) { moldyn->p_ref=p_ref; - printf("[moldyn] pressure: %f\n",moldyn->p_ref); + printf("[moldyn] pressure [atm]: %f\n",moldyn->p_ref/ATM); return 0; } @@ -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,12 +344,18 @@ 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); + set_nn_dist(moldyn,lc); + origin.x=0.5*lc; + origin.y=0.5*lc; + origin.z=0.5*lc; + ret=cubic_init(a,b,c,lc,atom,&origin); break; case FCC: - ret=fcc_init(a,b,c,lc,atom,&origin); + set_nn_dist(moldyn,0.5*sqrt(2.0)*lc); + ret=fcc_init(a,b,c,lc,atom,NULL); break; case DIAMOND: + set_nn_dist(moldyn,0.25*sqrt(3.0)*lc); ret=diamond_init(a,b,c,lc,atom,&origin); break; default: @@ -340,9 +391,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;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 +670,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 +677,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 +722,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,QUIET); potential_force_calc(moldyn); tp->x=(moldyn->energy-u)/moldyn->dv; p=tp->x*tp->x; @@ -667,6 +736,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,QUIET); potential_force_calc(moldyn); tp->y=(moldyn->energy-u)/moldyn->dv; p+=tp->y*tp->y; @@ -679,6 +750,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,QUIET); potential_force_calc(moldyn); tp->z=(moldyn->energy-u)/moldyn->dv; p+=tp->z*tp->z; @@ -687,15 +760,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,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); + 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 +781,9 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) { /* restore energy */ moldyn->energy=u; + link_cell_shutdown(moldyn); + link_cell_init(moldyn,QUIET); + return sqrt(p); } @@ -743,61 +823,48 @@ int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) { int scale_volume(t_moldyn *moldyn) { - t_atom *atom; t_3dvec *dim,*vdim; - double scale,v; - t_virial virial; + double scale; t_linkcell *lc; - int i; - atom=moldyn->atom; - dim=&(moldyn->dim); vdim=&(moldyn->vis.dim); + dim=&(moldyn->dim); lc=&(moldyn->lc); - memset(&virial,0,sizeof(t_virial)); + /* scaling factor */ + if(moldyn->pt_scale&P_SCALE_BERENDSEN) { + scale=1.0-(moldyn->p_ref-moldyn->p)/moldyn->p_tc; + scale=pow(scale,ONE_THIRD); + } + else { + scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD); + } +moldyn->debug=scale; + + /* scale the atoms and dimensions */ + scale_atoms(moldyn,scale,TRUE,TRUE,TRUE); + scale_dim(moldyn,scale,TRUE,TRUE,TRUE); - for(i=0;icount;i++) { - virial.xx+=atom[i].virial.xx; - virial.yy+=atom[i].virial.yy; - virial.zz+=atom[i].virial.zz; - virial.xy+=atom[i].virial.xy; - virial.xz+=atom[i].virial.xz; - virial.yz+=atom[i].virial.yz; + /* visualize dimensions */ + if(vdim->x!=0) { + vdim->x=dim->x; + vdim->y=dim->y; + vdim->z=dim->z; } - /* just a guess so far ... */ - v=virial.xx+virial.yy+virial.zz; - -printf("%f\n",v); - /* get pressure from virial */ - moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+ONE_THIRD*v; - moldyn->p/=moldyn->volume; -printf("%f | %f\n",moldyn->p/(ATM),moldyn->p_ref/ATM); - - /* scale factor */ - if(moldyn->pt_scale&P_SCALE_BERENDSEN) - scale=3*sqrt(1-(moldyn->p_ref-moldyn->p)/moldyn->p_tc); - else - /* should actually never be used */ - scale=pow(moldyn->p/moldyn->p_ref,1.0/3.0); - -printf("scale = %f\n",scale); - /* actual scaling */ - dim->x*=scale; - dim->y*=scale; - dim->z*=scale; - if(vdim->x) vdim->x=dim->x; - if(vdim->y) vdim->y=dim->y; - if(vdim->z) vdim->z=dim->z; - moldyn->volume*=(scale*scale*scale); - - /* check whether we need a new linkcell init */ - if((dim->x/moldyn->cutoff!=lc->nx)|| - (dim->y/moldyn->cutoff!=lc->ny)|| - (dim->z/moldyn->cutoff!=lc->nx)) { + /* recalculate scaled volume */ + moldyn->volume=dim->x*dim->y*dim->z; + + /* adjust/reinit linkcell */ + if(((int)(dim->x/moldyn->cutoff)!=lc->nx)|| + ((int)(dim->y/moldyn->cutoff)!=lc->ny)|| + ((int)(dim->z/moldyn->cutoff)!=lc->nx)) { link_cell_shutdown(moldyn); - link_cell_init(moldyn); + link_cell_init(moldyn,QUIET); + } else { + lc->x*=scale; + lc->y*=scale; + lc->z*=scale; } return 0; @@ -862,7 +929,7 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) { /* linked list / cell method */ -int link_cell_init(t_moldyn *moldyn) { +int link_cell_init(t_moldyn *moldyn,u8 vol) { t_linkcell *lc; int i; @@ -883,7 +950,7 @@ int link_cell_init(t_moldyn *moldyn) { if(lc->cells<27) printf("[moldyn] FATAL: less then 27 subcells!\n"); - printf("[moldyn] initializing linked cells (%d)\n",lc->cells); + if(vol) printf("[moldyn] initializing linked cells (%d)\n",lc->cells); for(i=0;icells;i++) list_init_f(&(lc->subcell[i])); @@ -1054,12 +1121,13 @@ int moldyn_integrate(t_moldyn *moldyn) { int fd; char dir[128]; double ds; + double energy_scale; sched=&(moldyn->schedule); atom=moldyn->atom; /* initialize linked cell method */ - link_cell_init(moldyn); + link_cell_init(moldyn,VERBOSE); /* logging & visualization */ e=moldyn->ewrite; @@ -1071,6 +1139,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); @@ -1109,23 +1180,26 @@ int moldyn_integrate(t_moldyn *moldyn) { /* integration step */ moldyn->integrate(moldyn); + /* calculate kinetic energy, temperature and pressure */ + update_e_kin(moldyn); + temperature_calc(moldyn); + pressure_calc(moldyn); + //thermodynamic_pressure_calc(moldyn); + /* p/t scaling */ if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT)) scale_velocity(moldyn,FALSE); if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT)) scale_volume(moldyn); -printf("-> %f\n",thermodynamic_pressure_calc(moldyn)/ATM); - /* 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 +1226,9 @@ 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, T: %f, P: %f V: %f", + sched->count,i, + moldyn->t,moldyn->p/ATM,moldyn->volume); fflush(stdout); } } @@ -1246,9 +1321,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 +1572,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;