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.0001*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;
}
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;
default:
- printf("[moldyn] unknown log mechanism: %02x\n",type);
+ printf("unknown log type: %02x\n",type);
return -1;
}
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:
+ ret=cubic_init(a,b,c,atom,&origin);
+ break;
case FCC:
ret=fcc_init(a,b,c,lc,atom,&origin);
break;
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
+
+HIER WEITER !
+}
+
/* 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++) {
atom=moldyn->atom;
+ double_ekin=0;
for(i=0;i<moldyn->count;i++)
double_ekin+=atom[i].mass*v3_absolute_square(&(atom[i].v));
return 0;
}
+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;
+}
+
double pressure_calc(t_moldyn *moldyn) {
int i;
- t_atom *atom;
- double p1,p2,p=0;
-
- for(i=0;i<moldyn->count;i++) {
-
+ double p1,p;
+ double v,h;
+ t_virial *virial;
+ v=0.0;
+ for(i=0;i<moldyn->count;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;
- p2=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt2);
- p2/=moldyn->volume;
+ printf("debug: vt1=%f v=%f nkt=%f\n",moldyn->vt1,v,h);
- printf("compare pressures: %f %f\n",p1/ATM,p2/ATM);
+ printf("compare pressures: %f %f %f\n",p1/ATM,p/ATM,h/moldyn->volume/ATM);
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;
+ }
+
+ /* save unscaled potential energy + atom/dim configuration */
+ u=moldyn->energy;
+ memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
+ 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);
+ potential_force_calc(moldyn);
+ tp->x=(moldyn->energy-u)/moldyn->dv;
+ p=tp->x*tp->x;
+
+ /* 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);
+ 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);
+ 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\n",(tp->x+tp->y+tp->z)/ATM);
+
+ scale=1.0+pow(moldyn->dv/moldyn->volume,ONE_THIRD);
+
+ scale_dim(moldyn,scale,1,1,1);
+ scale_dim(moldyn,scale,1,1,1);
+ potential_force_calc(moldyn);
+
+ 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;
+
+ return sqrt(p);
+}
+
double get_pressure(t_moldyn *moldyn) {
return moldyn->p;
}
+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;
/* 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);
+printf("-> %f\n",thermodynamic_pressure_calc(moldyn)/ATM);
+
/* 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);
}
if(m) {
if(!(i%m)) {
/* reset energy */
moldyn->energy=0.0;
- moldyn->vt2=0.0;
-
- /* get energy and force of every atom */
+ /* reset virial */
+ moldyn->vt1=0.0;
+
+ /* 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]));
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);
-
+temperature_calc(moldyn);
pressure_calc(moldyn);
+ideal_gas_law_pressure(moldyn);
return 0;
}