#include <math.h>
#include "moldyn.h"
+#include "report/report.h"
int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
+ printf("[moldyn] init\n");
+
memset(moldyn,0,sizeof(t_moldyn));
rand_init(&(moldyn->random),NULL,1);
int moldyn_shutdown(t_moldyn *moldyn) {
printf("[moldyn] shutdown\n");
+
moldyn_log_shutdown(moldyn);
link_cell_shutdown(moldyn);
rand_close(&(moldyn->random));
int set_int_alg(t_moldyn *moldyn,u8 algo) {
+ printf("[moldyn] integration algorithm: ");
+
switch(algo) {
case MOLDYN_INTEGRATE_VERLET:
moldyn->integrate=velocity_verlet;
+ printf("velocity verlet\n");
break;
default:
printf("unknown integration algorithm: %02x\n",algo);
+ printf("unknown\n");
return -1;
}
moldyn->cutoff=cutoff;
+ printf("[moldyn] cutoff [A]: %f\n",moldyn->cutoff);
+
return 0;
}
moldyn->t_ref=t_ref;
+ printf("[moldyn] temperature: %f\n",moldyn->t_ref);
+
return 0;
}
moldyn->p_ref=p_ref;
+ printf("[moldyn] pressure: %f\n",moldyn->p_ref);
+
return 0;
}
moldyn->t_tc=ttc;
moldyn->p_tc=ptc;
+ printf("[moldyn] p/t scaling:\n");
+
+ printf(" p: %s",ptype?"yes":"no ");
+ if(ptype)
+ printf(" | type: %02x | factor: %f",ptype,ptc);
+ printf("\n");
+
+ printf(" t: %s",ttype?"yes":"no ");
+ if(ttype)
+ printf(" | type: %02x | factor: %f",ttype,ttc);
+ printf("\n");
+
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?"on":"off");
+ printf(" visualize simulation box: %s\n",visualize?"yes":"no");
+ printf(" delta volume (pressure calc): %f\n",moldyn->dv);
return 0;
}
int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
+ printf("[moldyn] periodic boundary conditions:\n");
+
if(x)
moldyn->status|=MOLDYN_STAT_PBX;
if(z)
moldyn->status|=MOLDYN_STAT_PBZ;
+ printf(" x: %s\n",x?"yes":"no");
+ printf(" y: %s\n",y?"yes":"no");
+ printf(" z: %s\n",z?"yes":"no");
+
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) {
char filename[128];
int ret;
+ printf("[moldyn] set log: ");
+
switch(type) {
case LOG_TOTAL_ENERGY:
moldyn->ewrite=timer;
return moldyn->efd;
}
dprintf(moldyn->efd,"# total energy log file\n");
+ printf("total energy (%d)\n",timer);
break;
case LOG_TOTAL_MOMENTUM:
moldyn->mwrite=timer;
return moldyn->mfd;
}
dprintf(moldyn->efd,"# total momentum log file\n");
+ printf("total momentum (%d)\n",timer);
break;
case SAVE_STEP:
moldyn->swrite=timer;
+ printf("save file (%d)\n",timer);
break;
case VISUAL_STEP:
moldyn->vwrite=timer;
printf("[moldyn] visual init failure\n");
return ret;
}
+ 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("[moldyn] unknown log mechanism: %02x\n",type);
+ 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:
+ 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);
+ ret=fcc_init(a,b,c,lc,atom,NULL);
break;
case DIAMOND:
ret=diamond_init(a,b,c,lc,atom,&origin);
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) {
atom=moldyn->atom;
random=&(moldyn->random);
+ printf("[moldyn] thermal init (equi init: %s)\n",equi_init?"yes":"no");
+
/* gaussian distribution of velocities */
v3_zero(&p_total);
for(i=0;i<moldyn->count;i++) {
double temperature_calc(t_moldyn *moldyn) {
- double double_ekin;
- int i;
- t_atom *atom;
-
- atom=moldyn->atom;
-
- for(i=0;i<moldyn->count;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;
}
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);
+ 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);
+ 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);
+ 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);
+ 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);
+
+ 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;
schedule->tau=ptr;
schedule->tau[count-1]=tau;
+ printf("[moldyn] schedule added:\n");
+ printf(" number: %d | runs: %d | tau: %f\n",count-1,runs,tau);
+
+
return 0;
}
-int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params) {
+int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
moldyn->schedule.hook=hook;
moldyn->schedule.hook_params=hook_params;
int fd;
char dir[128];
double ds;
+ double energy_scale;
sched=&(moldyn->schedule);
atom=moldyn->atom;
moldyn->tau_square=moldyn->tau*moldyn->tau;
moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
+ /* energy scaling factor */
+ energy_scale=moldyn->count*EV;
+
+printf("debug: %f\n",moldyn->atom[0].f.x);
/* calculate initial forces */
potential_force_calc(moldyn);
+printf("debug: %f\n",moldyn->atom[0].f.x);
/* some stupid checks before we actually start calculating bullshit */
if(moldyn->cutoff>0.5*moldyn->dim.x)
/* debugging, ignore */
moldyn->debug=0;
+ /* tell the world */
+ printf("[moldyn] integration start, go get a coffee ...\n");
+
/* executing the schedule */
for(sched->count=0;sched->count<sched->total_sched;sched->count++) {
if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
scale_volume(moldyn);
+ update_e_kin(moldyn);
+ temperature_calc(moldyn);
+ pressure_calc(moldyn);
+ //thermodynamic_pressure_calc(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, debug: %f",
+ sched->count,i,moldyn->p/ATM);
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