return 0;
}
-int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params) {
+int set_potential1b(t_moldyn *moldyn,pf_func1b func) {
moldyn->func1b=func;
- moldyn->pot1b_params=params;
return 0;
}
-int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params) {
+int set_potential2b(t_moldyn *moldyn,pf_func2b func) {
moldyn->func2b=func;
- moldyn->pot2b_params=params;
return 0;
}
-int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params) {
+int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func) {
moldyn->func2b_post=func;
- moldyn->pot2b_params=params;
return 0;
}
-int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) {
+int set_potential3b(t_moldyn *moldyn,pf_func3b func) {
moldyn->func3b=func;
- moldyn->pot3b_params=params;
+
+ return 0;
+}
+
+int set_potential_params(t_moldyn *moldyn,void *params) {
+
+ moldyn->pot_params=params;
return 0;
}
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);
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);
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);
t_virial *virial;
/*
- * P = 1/(3V) sum_i ( p_i^2 / 2m + f_i r_i )
- *
- * virial = f_i r_i
+ * PV = NkT + <W>
+ * W = 1/3 sum_i f_i r_i
+ * virial = sum_i f_i r_i
+ *
+ * => P = (2 Ekin + virial) / (3V)
*/
v=0.0;
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) {
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 */
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;
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;
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;
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))
}
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) {
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;
* virial calculation
*/
-inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
+//inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
+int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
a->virial.xx+=f->x*d->x;
a->virial.yy+=f->y*d->y;
* periodic boundayr checking
*/
-inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
+//inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
+int check_per_bound(t_moldyn *moldyn,t_3dvec *a) {
double x,y,z;
t_3dvec *dim;
return 0;
}
-
-/*
- * example potentials
- */
-
-/* harmonic oscillator potential and force */
-
-int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
-
- t_ho_params *params;
- t_3dvec force,distance;
- double d,f;
- double sc,equi_dist;
-
- params=moldyn->pot2b_params;
- sc=params->spring_constant;
- equi_dist=params->equilibrium_distance;
-
- if(ai<aj) return 0;
-
- v3_sub(&distance,&(aj->r),&(ai->r));
-
- if(bc) check_per_bound(moldyn,&distance);
- d=v3_norm(&distance);
- if(d<=moldyn->cutoff) {
- moldyn->energy+=(0.5*sc*(d-equi_dist)*(d-equi_dist));
- /* f = -grad E; grad r_ij = -1 1/r_ij distance */
- f=sc*(1.0-equi_dist/d);
- v3_scale(&force,&distance,f);
- v3_add(&(ai->f),&(ai->f),&force);
- virial_calc(ai,&force,&distance);
- virial_calc(aj,&force,&distance); /* f and d signe switched */
- v3_scale(&force,&distance,-f);
- v3_add(&(aj->f),&(aj->f),&force);
- }
-
- return 0;
-}
-
/*
* debugging / critical check functions
*/
projects / physik / posic.git / blobdiff