#include <math.h>
#include "moldyn.h"
+#include "report/report.h"
int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
moldyn->t_ref=t_ref;
- printf("[moldyn] temperature: %f\n",moldyn->t_ref);
+ printf("[moldyn] temperature [K]: %f\n",moldyn->t_ref);
return 0;
}
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;
}
moldyn->vis.dim.z=z;
}
+ moldyn->dv=0.000001*moldyn->volume;
+
printf("[moldyn] dimensions in A and A^3 respectively:\n");
printf(" x: %f\n",moldyn->dim.x);
printf(" y: %f\n",moldyn->dim.y);
printf(" z: %f\n",moldyn->dim.z);
printf(" volume: %f\n",moldyn->volume);
printf(" visualize simulation box: %s\n",visualize?"yes":"no");
+ printf(" delta volume (pressure calc): %f\n",moldyn->dv);
return 0;
}
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) {
}
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;
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;
count=moldyn->count;
/* how many atoms do we expect */
+ if(type==CUBIC) new*=1;
if(type==FCC) new*=4;
if(type==DIAMOND) new*=8;
v3_zero(&origin);
switch(type) {
+ case CUBIC:
+ 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:
return ret;
}
+/* cubic init */
+int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
+
+ int count;
+ 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;i<a;i++) {
+ r.y=o.y;
+ for(j=0;j<b;j++) {
+ r.z=o.z;
+ for(k=0;k<c;k++) {
+ v3_copy(&(atom[count].r),&r);
+ count+=1;
+ r.z+=lc;
+ }
+ r.y+=lc;
+ }
+ r.x+=lc;
+ }
+
+ for(i=0;i<count;i++) {
+ atom[i].r.x-=(a*lc)/2.0;
+ atom[i].r.y-=(b*lc)/2.0;
+ atom[i].r.z-=(c*lc)/2.0;
+ }
+
+ return count;
+}
+
/* fcc lattice init */
int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
double temperature_calc(t_moldyn *moldyn) {
- double double_ekin;
- int i;
- t_atom *atom;
+ /* assume up to date kinetic energy, which is 3/2 N k_B T */
- atom=moldyn->atom;
-
- for(i=0;i<moldyn->count;i++)
- double_ekin+=atom[i].mass*v3_absolute_square(&(atom[i].v));
-
- /* 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;
}
return 0;
}
+double ideal_gas_law_pressure(t_moldyn *moldyn) {
+
+ double p;
+
+ p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume;
+
+ return p;
+}
+
double pressure_calc(t_moldyn *moldyn) {
int i;
- t_atom *atom;
- double p1,p2,p=0;
-
+ 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;i<moldyn->count;i++) {
-
+ virial=&(moldyn->atom[i].virial);
+ v+=(virial->xx+virial->yy+virial->zz);
+ }
+
+ /* assume up to date kinetic energy */
+ moldyn->p=2.0*moldyn->ekin+v;
+ moldyn->p/=(3.0*moldyn->volume);
+
+ return moldyn->p;
+}
+
+double thermodynamic_pressure_calc(t_moldyn *moldyn) {
+
+ t_3dvec dim,*tp;
+ double u,p;
+ double scale;
+ t_atom *store;
+ tp=&(moldyn->tp);
+ store=malloc(moldyn->count*sizeof(t_atom));
+ if(store==NULL) {
+ printf("[moldyn] allocating store mem failed\n");
+ return -1;
}
- p1=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt1);
- p1/=moldyn->volume;
+ /* save unscaled potential energy + atom/dim configuration */
+ u=moldyn->energy;
+ memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
+ dim=moldyn->dim;
- p2=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt2);
- p2/=moldyn->volume;
+ /* 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);
+ 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;
- printf("compare pressures: %f %f\n",p1/ATM,p2/ATM);
+ /* restore atomic configuration + dim */
+ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
+ moldyn->dim=dim;
- return moldyn->p;
-}
+ /* 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);
+ 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;
+
+ /* restore atomic configuration + dim */
+ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
+ 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);
+ 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;
+
+ /* 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;
+
+ link_cell_shutdown(moldyn);
+ link_cell_init(moldyn,QUIET);
+
+ return sqrt(p);
+}
double get_pressure(t_moldyn *moldyn) {
}
+int scale_dim(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) {
+
+ t_3dvec *dim;
+
+ dim=&(moldyn->dim);
+
+ if(x) dim->x*=scale;
+ if(y) dim->y*=scale;
+ if(z) dim->z*=scale;
+
+ return 0;
+}
+
+int scale_atoms(t_moldyn *moldyn,double scale,u8 x,u8 y,u8 z) {
+
+ int i;
+ t_3dvec *r;
+
+ for(i=0;i<moldyn->count;i++) {
+ r=&(moldyn->atom[i].r);
+ if(x) r->x*=scale;
+ if(y) r->y*=scale;
+ if(z) r->z*=scale;
+ }
+
+ return 0;
+}
+
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;
- for(i=0;i<moldyn->count;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;
+ /* scale the atoms and dimensions */
+ scale_atoms(moldyn,scale,TRUE,TRUE,TRUE);
+ scale_dim(moldyn,scale,TRUE,TRUE,TRUE);
+
+ /* 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;
/* 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;
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;i<lc->cells;i++)
list_init_f(&(lc->subcell[i]));
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;
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);
/* 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);
scale_volume(moldyn);
/* check for log & visualization */
-//double ax;
-//double ao;
-//double av;
if(e) {
if(!(i%e))
-//ao=sqrt(0.1/M_SI);
-//ax=((0.28-0.25)*sqrt(3)*LC_SI/2)*cos(ao*i);
-//av=ao*(0.28-0.25)*sqrt(3)*LC_SI/2*sin(ao*i);
- update_e_kin(moldyn);
dprintf(moldyn->efd,
"%f %f %f %f\n",
- moldyn->time,moldyn->ekin,
- moldyn->energy,
- get_total_energy(moldyn));
-//moldyn->atom[0].r.x,ax,av*av*M_SI,0.1*ax*ax,av*av*M_SI+0.1*ax*ax);
+ moldyn->time,moldyn->ekin/energy_scale,
+ moldyn->energy/energy_scale,
+ get_total_energy(moldyn)/energy_scale);
}
if(m) {
if(!(i%m)) {
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);
}
}
/* reset energy */
moldyn->energy=0.0;
- moldyn->vt2=0.0;
-
- /* get energy and force of every atom */
+ /* reset force, site energy and virial of every atom */
for(i=0;i<count;i++) {
/* reset force */
v3_zero(&(itom[i].f));
- /* reset viral of atom i */
- virial=&(itom[i].virial);
+ /* reset virial */
+ virial=(&(itom[i].virial));
virial->xx=0.0;
virial->yy=0.0;
virial->zz=0.0;
virial->xy=0.0;
virial->xz=0.0;
virial->yz=0.0;
- moldyn->vt1=0.0;
-
+
/* reset site energy */
itom[i].e=0.0;
+ }
+
+ /* get energy,force and virial of every atom */
+ for(i=0;i<count;i++) {
+
/* single particle potential/force */
if(itom[i].attr&ATOM_ATTR_1BP)
moldyn->func1b(moldyn,&(itom[i]));
}
}
+
#ifdef DEBUG
printf("\n\n");
#endif
printf("\n\n");
#endif
- moldyn->vt2=0.0;
- for(i=0;i<count;i++)
- moldyn->vt2-=v3_scalar_product(&(itom[i].r),&(itom[i].f));
-
-//printf("compare: vt1: %f vt2: %f\n",moldyn->vt1,moldyn->vt2);
-
-//pressure_calc(moldyn);
-
return 0;
}
inline 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;
- 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;
}
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;
projects / physik / posic.git / blobdiff