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);
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;
switch(type) {
case CUBIC:
- ret=cubic_init(a,b,c,atom,&origin);
+ ret=cubic_init(a,b,c,lc,atom,NULL);
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);
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;i<a;i++) {
+ r.y=o.y;
+ for(j=0;j<b;j++) {
+ for(k=0;k<c;k++) {
+ r.z=o.z;
+ 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 */
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;
- double p1,p;
- double v,h;
+ double v;
t_virial *virial;
v=0.0;
virial=&(moldyn->atom[i].virial);
v+=(virial->xx+virial->yy+virial->zz);
}
+ v*=ONE_THIRD;
+printf("kieck mal: %f %f %f\n",v,moldyn->count*K_BOLTZMANN*moldyn->t,v/moldyn->count);
- 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);
+ moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+v;
+ moldyn->p/=moldyn->volume;
return moldyn->p;
}
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;
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;
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;
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);
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));
/* restore energy */
moldyn->energy=u;
+ link_cell_shutdown(moldyn);
+ link_cell_init(moldyn);
+
return sqrt(p);
}
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;
+
/* calculate initial forces */
potential_force_calc(moldyn);
if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
scale_volume(moldyn);
-printf("-> %f\n",thermodynamic_pressure_calc(moldyn)/ATM);
+ temperature_calc(moldyn);
+ pressure_calc(moldyn);
+ //thermodynamic_pressure_calc(moldyn);
/* check for log & visualization */
if(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)) {
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 | %f",
+ sched->count,i,moldyn->p/ATM,moldyn->debug/ATM);
fflush(stdout);
}
}
/* reset energy */
moldyn->energy=0.0;
- /* reset virial */
- moldyn->vt1=0.0;
-
/* reset force, site energy and virial of every atom */
for(i=0;i<count;i++) {
}
}
+
#ifdef DEBUG
printf("\n\n");
#endif
printf("\n\n");
#endif
-temperature_calc(moldyn);
-pressure_calc(moldyn);
-ideal_gas_law_pressure(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;
}
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;
#include "random/random.h"
#include "list/list.h"
+#include "report/report.h"
/*
*
t_3dvec dim; /* dimensions of the simulation volume */
double volume; /* volume of sim cell (dim.x*dim.y*dim.z) */
- double vt1,vt2;
/* potential force function and parameter pointers */
int (*func1b)(struct s_moldyn *moldyn,t_atom *ai);
int mfd; /* fd for momentum log */
unsigned int vwrite; /* how often to visualize atom information */
unsigned int swrite; /* how often to create a save file */
+ int rfd; /* report file descriptor */
+ char rtitle[64]; /* report title */
+ char rauthor[64]; /* report author */
+ int pfd; /* gnuplot script file descriptor */
u8 status; /* general moldyn properties */
#define P_SCALE_BERENDSEN 0x04 /* berendsen p control */
#define P_SCALE_DIRECT 0x08 /* direct p control */
-
/*
*
* potential parameter structures
#define LOG_TOTAL_MOMENTUM 0x02
#define SAVE_STEP 0x04
#define VISUAL_STEP 0x08
+#define CREATE_REPORT 0x10
#define TRUE 1
#define FALSE 0
* lattice constants
*/
-#define FCC 0x01
-#define DIAMOND 0x02
+#define CUBIC 0x01
+#define FCC 0x02
+#define DIAMOND 0x04
/*
int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params);
int moldyn_set_log_dir(t_moldyn *moldyn,char *dir);
+int moldyn_set_report(t_moldyn *moldyn,char *author,char *title);
int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer);
int moldyn_log_shutdown(t_moldyn *moldyn);
int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
u8 attr,u8 brand,int a,int b,int c);
+int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);
int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);
int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);
int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
if(md->schedule.count==1)
set_pt_scale(md,0,0,0,0);
-
//set_temperature(md,md->t_ref-100.0);
-
return 0;
}
/* cutoff radius */
//set_cutoff(&md,TM_S_SI);
- set_cutoff(&md,2*LC_SI);
+ set_cutoff(&md,2*LC_SI*0.5*sqrt(1.5));
/*
* potential parameters
tersoff_mult_complete_params(&tp);
/* set (initial) dimensions of simulation volume */
- set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
+ set_dim(&md,8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),TRUE);
/* set periodic boundary conditions in all directions */
set_pbc(&md,TRUE,TRUE,TRUE);
/* create the lattice / place atoms */
- create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+ create_lattice(&md,FCC,LC_SI*0.5*sqrt(1.5),SI,M_SI,
+ //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
// ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
ATOM_ATTR_2BP|ATOM_ATTR_HB,
- 0,6,6,6);
+ 0,8,8,8);
moldyn_bc_check(&md);
/* testing configuration */
/* set p/t scaling */
//set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,
// T_SCALE_BERENDSEN,100.0);
- set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
+ //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
//set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,0,0);
/* initial thermal fluctuations of particles (in equilibrium) */
thermal_init(&md,TRUE);
/* create the simulation schedule */
- moldyn_add_schedule(&md,1001,3.0);
- moldyn_add_schedule(&md,501,1.0);
- moldyn_add_schedule(&md,501,1.0);
+ moldyn_add_schedule(&md,1001,1.0);
+ //moldyn_add_schedule(&md,501,1.0);
+ //moldyn_add_schedule(&md,501,1.0);
/* schedule hook function */
- moldyn_set_schedule_hook(&md,&hook,NULL);
+ //moldyn_set_schedule_hook(&md,&hook,NULL);
/* activate logging */
moldyn_set_log_dir(&md,argv[1]);
- moldyn_set_log(&md,LOG_TOTAL_ENERGY,5);
- moldyn_set_log(&md,VISUAL_STEP,50);
+ moldyn_set_report(&md,"Frank Zirkelbach","Test 1");
+ moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
+ moldyn_set_log(&md,VISUAL_STEP,10);
+ moldyn_set_log(&md,CREATE_REPORT,0);
/*
* let's do the actual md algorithm now