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;
projects / physik / posic.git / blobdiff