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->atom=ptr;
atom=&(moldyn->atom[count]);
-
- v3_zero(&origin);
+
+ /* no atoms on the boundaries (only reason: it looks better!) */
+ origin.x=0.5*lc;
+ origin.y=0.5*lc;
+ origin.z=0.5*lc;
switch(type) {
case CUBIC:
- ret=cubic_init(a,b,c,lc,atom,NULL);
+ set_nn_dist(moldyn,lc);
+ ret=cubic_init(a,b,c,lc,atom,&origin);
break;
case FCC:
- ret=fcc_init(a,b,c,lc,atom,NULL);
+ v3_scale(&origin,&origin,0.5);
+ set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
+ ret=fcc_init(a,b,c,lc,atom,&origin);
break;
case DIAMOND:
+ v3_scale(&origin,&origin,0.25);
+ set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
ret=diamond_init(a,b,c,lc,atom,&origin);
break;
default:
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++) {
- r.z=o.z;
v3_copy(&(atom[count].r),&r);
count+=1;
r.z+=lc;
int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
int count;
- int i,j;
+ int i,j,k,l;
t_3dvec o,r,n;
t_3dvec basis[3];
- double help[3];
- double x,y,z;
-
- x=a*lc;
- y=b*lc;
- z=c*lc;
- if(origin) v3_copy(&o,origin);
- else v3_zero(&o);
+ count=0;
+ if(origin)
+ v3_copy(&o,origin);
+ else
+ v3_zero(&o);
/* construct the basis */
- for(i=0;i<3;i++) {
- for(j=0;j<3;j++) {
- if(i!=j) help[j]=0.5*lc;
- else help[j]=.0;
- }
- v3_set(&basis[i],help);
- }
+ memset(basis,0,3*sizeof(t_3dvec));
+ basis[0].x=0.5*lc;
+ basis[0].y=0.5*lc;
+ basis[1].x=0.5*lc;
+ basis[1].z=0.5*lc;
+ basis[2].y=0.5*lc;
+ basis[2].z=0.5*lc;
- v3_zero(&r);
- count=0;
-
/* fill up the room */
r.x=o.x;
- while(r.x<x) {
+ for(i=0;i<a;i++) {
r.y=o.y;
- while(r.y<y) {
+ for(j=0;j<b;j++) {
r.z=o.z;
- while(r.z<z) {
+ for(k=0;k<c;k++) {
+ /* first atom */
v3_copy(&(atom[count].r),&r);
- atom[count].element=1;
count+=1;
- for(i=0;i<3;i++) {
- v3_add(&n,&r,&basis[i]);
- if((n.x<x+o.x)&&
- (n.y<y+o.y)&&
- (n.z<z+o.z)) {
- v3_copy(&(atom[count].r),&n);
- count+=1;
- }
+ r.z+=lc;
+ /* the three face centered atoms */
+ for(l=0;l<3;l++) {
+ v3_add(&n,&r,&basis[l]);
+ v3_copy(&(atom[count].r),&n);
+ count+=1;
}
- r.z+=lc;
}
r.y+=lc;
}
r.x+=lc;
}
-
+
/* coordinate transformation */
- help[0]=x/2.0;
- help[1]=y/2.0;
- help[2]=z/2.0;
- v3_set(&n,help);
- for(i=0;i<count;i++)
- v3_sub(&(atom[i].r),&(atom[i].r),&n);
-
+ 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;
}
scale_dim(moldyn,scale,TRUE,0,0);
scale_atoms(moldyn,scale,TRUE,0,0);
link_cell_shutdown(moldyn);
- link_cell_init(moldyn);
+ link_cell_init(moldyn,QUIET);
potential_force_calc(moldyn);
tp->x=(moldyn->energy-u)/moldyn->dv;
p=tp->x*tp->x;
scale_dim(moldyn,scale,0,TRUE,0);
scale_atoms(moldyn,scale,0,TRUE,0);
link_cell_shutdown(moldyn);
- link_cell_init(moldyn);
+ link_cell_init(moldyn,QUIET);
potential_force_calc(moldyn);
tp->y=(moldyn->energy-u)/moldyn->dv;
p+=tp->y*tp->y;
scale_dim(moldyn,scale,0,0,TRUE);
scale_atoms(moldyn,scale,0,0,TRUE);
link_cell_shutdown(moldyn);
- link_cell_init(moldyn);
+ link_cell_init(moldyn,QUIET);
potential_force_calc(moldyn);
tp->z=(moldyn->energy-u)/moldyn->dv;
p+=tp->z*tp->z;
scale_dim(moldyn,scale,1,1,1);
scale_atoms(moldyn,scale,1,1,1);
link_cell_shutdown(moldyn);
- link_cell_init(moldyn);
+ link_cell_init(moldyn,QUIET);
potential_force_calc(moldyn);
printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x);
moldyn->energy=u;
link_cell_shutdown(moldyn);
- link_cell_init(moldyn);
+ link_cell_init(moldyn,QUIET);
return sqrt(p);
}
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]));
atom=moldyn->atom;
/* initialize linked cell method */
- link_cell_init(moldyn);
+ link_cell_init(moldyn,VERBOSE);
/* logging & visualization */
e=moldyn->ewrite;
/* 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);
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 */
if(e) {
if(!(i%e))
if(!(i%v)) {
visual_atoms(&(moldyn->vis),moldyn->time,
moldyn->atom,moldyn->count);
- printf("\rsched: %d, steps: %d, debug: %f | %f",
- sched->count,i,moldyn->p/ATM,moldyn->p/ATM);
+ printf("\rsched: %d, steps: %d, T: %f, P: %f V: %f",
+ sched->count,i,
+ moldyn->t,moldyn->p/ATM,moldyn->volume);
fflush(stdout);
}
}
v3_scale(&force,&force,-1.0); /* f = - grad E */
v3_add(&(ai->f),&(ai->f),&force);
virial_calc(ai,&force,&distance);
+if(force.x*distance.x<=0) printf("virial xx: %.15f -> %f %f %f\n",force.x*distance.x,distance.x,distance.y,distance.z);
virial_calc(aj,&force,&distance); /* f and d signe switched */
}
projects / physik / posic.git / blobdiff