return 0;
}
+int set_pressure(t_moldyn *moldyn,double p_ref) {
+
+ moldyn->p_ref=p_ref;
+
+ return 0;
+}
+
int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc) {
moldyn->pt_scale=(ptype|ttype);
moldyn->dim.y=y;
moldyn->dim.z=z;
+ moldyn->volume=x*y*z;
+
if(visualize) {
moldyn->vis.dim.x=x;
moldyn->vis.dim.y=y;
moldyn->vis.dim.z=z;
}
- printf("[moldyn] dimensions in A:\n");
+ printf("[moldyn] dimensions in A and A^2 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");
return 0;
count+=1;
}
}
- if(count!=0) moldyn->t=(2.0*e)/(3.0*count*K_BOLTZMANN);
+ if(count!=0) moldyn->t=e/(1.5*count*K_BOLTZMANN);
else return 0; /* no atoms involved in scaling! */
/* (temporary) hack for e,t = 0 */
return 0;
}
+int scale_volume(t_moldyn *moldyn) {
+
+ t_atom *atom;
+ t_3dvec *dim,*vdim;
+ double virial,scale;
+ t_linkcell *lc;
+ int i;
+
+ atom=moldyn->atom;
+ dim=&(moldyn->dim);
+ vdim=&(moldyn->vis.dim);
+ lc=&(moldyn->lc);
+
+ for(i=0;i<moldyn->count;i++)
+ virial+=v3_norm(&(atom[i].virial));
+
+printf("%f\n",virial);
+ /* get pressure from virial */
+ moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*virial;
+ moldyn->p/=moldyn->volume;
+printf("%f\n",moldyn->p/(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)) {
+ link_cell_shutdown(moldyn);
+ link_cell_init(moldyn);
+ }
+
+ return 0;
+
+}
+
double get_e_kin(t_moldyn *moldyn) {
int i;
/* 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))
+{
+printf("going to do p scale ...\n");
+ scale_volume(moldyn);
+printf("done\n");
+}
/* check for log & visualization */
if(e) {
/* reset energy */
moldyn->energy=0.0;
-
+
/* get energy and force of every atom */
for(i=0;i<count;i++) {
/* reset force */
v3_zero(&(itom[i].f));
+ /* reset viral of atom i */
+ v3_zero(&(itom[i].virial));
+
/* single particle potential/force */
if(itom[i].attr&ATOM_ATTR_1BP)
moldyn->func1b(moldyn,&(itom[i]));
v3_scale(&temp2,&temp2,tmp*B*exp(-mu*d_jk)*f_c_jk*0.5);
v3_add(&(ai->f),&(ai->f),&temp2); /* -1 skipped in f_a calc ^ */
/* scaled with 0.5 ^ */
+
}
return 0;
printf("FATAL: atom %d: x: %.20f (%.20f)\n",
i,atom[i].r.x,dim->x/2);
printf("diagnostic:\n");
+ printf("-----------\natom.r.x:\n");
for(j=0;j<8;j++) {
memcpy(&byte,(u8 *)(&(atom[i].r.x))+j,1);
for(k=0;k<8;k++)
((byte)&(1<<k))?1:0,
(k==7)?'\n':'|');
}
- printf("---------------\n");
+ printf("---------------\nx=dim.x/2:\n");
for(j=0;j<8;j++) {
memcpy(&byte,(u8 *)(&x)+j,1);
for(k=0;k<8;k++)
((byte)&(1<<k))?1:0,
(k==7)?'\n':'|');
}
- if(atom[i].r.x==x) printf("hier gleich!\n");
- else printf("hier NICHT gleich!\n");
+ if(atom[i].r.x==x) printf("the same!\n");
+ else printf("different!\n");
}
if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
printf("FATAL: atom %d: y: %.20f (%.20f)\n",
t_3dvec r; /* position */
t_3dvec v; /* velocity */
t_3dvec f; /* force */
+ t_3dvec virial; /* virial */
int element; /* number of element in pse */
double mass; /* atom mass */
u8 bnum; /* brand number */
t_atom *atom; /* pointer to the atoms */
t_3dvec dim; /* dimensions of the simulation volume */
+ double volume; /* volume of sim cell (dim.x*dim.y*dim.z) */
/* potential force function and parameter pointers */
int (*func1b)(struct s_moldyn *moldyn,t_atom *ai);
double p; /* actual pressure */
/* pressure and temperature control (velocity/volume scaling) */
- /* (in units of tau!) */
+ /* (t_tc in units of tau, p_tc in units of tau * isoth. compressib.) */
unsigned char pt_scale; /* type of p and t scaling */
double t_tc; /* t berendsen control time constant */
double p_tc; /* p berendsen control time constant */
int debug; /* debugging stuff, ignore */
} t_moldyn;
-#define MOLDYN_STAT_PBX 0x08 /* periodic boudaries in x */
-#define MOLDYN_STAT_PBY 0x10 /* y */
-#define MOLDYN_STAT_PBZ 0x20 /* and z direction */
+#define MOLDYN_STAT_PBX 0x01 /* periodic boudaries in x */
+#define MOLDYN_STAT_PBY 0x02 /* y */
+#define MOLDYN_STAT_PBZ 0x04 /* and z direction */
-#define MOLDYN_1BP 0x00 /* care about single */
-#define MOLDYN_2BP 0x01 /* 2 body */
-#define MOLDYN_3BP 0x02 /* and 3 body particle pots */
+#define MOLDYN_PSCALE 0x08 /* size controlled by piston */
+
+#define MOLDYN_1BP 0x10 /* care about single */
+#define MOLDYN_2BP 0x20 /* 2 body */
+#define MOLDYN_3BP 0x40 /* and 3 body particle pots */
#define T_SCALE_BERENDSEN 0x01 /* berendsen t control */
#define T_SCALE_DIRECT 0x02 /* direct t control */
*
*/
+#define ONE_THIRD (1.0/3.0)
+
/*
* default values
*
* - length unit: 1 A (1 A = 1e-10 m)
* - time unit: 1 fs (1 fs = 1e-15 s)
* - mass unit: 1 amu (1 amu = 1.6605388628e-27 kg )
+ *
+ * fyi: in the following 1 N = (amu*A)/(fs*fs)
+ *
*/
#define METER 1e10 /* A */
#define AMU 1.6605388628e-27 /* kg */
#define KILOGRAM (1.0/AMU) /* amu */
#define NEWTON (METER*KILOGRAM/(SECOND*SECOND)) /* A amu / fs^2 */
+#define PASCAL (NEWTON/(METER*METER)) /* N / A^2 */
+#define ATM (1.0133e5*PASCAL) /* N / A^2 */
#define MOLDYN_TEMP 273.0
#define MOLDYN_TAU 1.0
*
* phsical values / constants
*
+ *
*/
#define K_BOLTZMANN (1.380650524e-23*METER*NEWTON) /* NA/K */
int set_int_alg(t_moldyn *moldyn,u8 algo);
int set_cutoff(t_moldyn *moldyn,double cutoff);
int set_temperature(t_moldyn *moldyn,double t_ref);
+int set_pressure(t_moldyn *moldyn,double p_ref);
int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc);
int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize);
int set_nn_dist(t_moldyn *moldyn,double dist);
int thermal_init(t_moldyn *moldyn,u8 equi_init);
int scale_velocity(t_moldyn *moldyn,u8 equi_init);
+int scale_volume(t_moldyn *moldyn);
double get_e_kin(t_moldyn *moldyn);
double get_e_pot(t_moldyn *moldyn);
create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
0,5,5,5);
+ moldyn_bc_check(&md);
/* testing configuration */
//r.x=2.45/2; v.x=0;
/* set temperature */
printf("[sic] setting temperature\n");
+ set_temperature(&md,273.0+1410.0);
//set_temperature(&md,273.0+450.0);
- set_temperature(&md,1.0);
//set_temperature(&md,273.0);
+ //set_temperature(&md,1.0);
+ //set_temperature(&md,0.0);
+
+ /* set pressure */
+ printf("[sic] setting pressure\n");
+ set_pressure(&md,ATM);
/* set p/t scaling */
printf("[sic] set p/t scaling\n");
- set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
+ //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,
+ // T_SCALE_BERENDSEN,100.0);
/* initial thermal fluctuations of particles (in equilibrium) */
printf("[sic] thermal init\n");
/* create the simulation schedule */
printf("[sic] adding schedule\n");
- moldyn_add_schedule(&md,1000,1.0);
+ moldyn_add_schedule(&md,100001,1.0);
/* activate logging */
printf("[sic] activate logging\n");
moldyn_set_log_dir(&md,"saves/test");
- moldyn_set_log(&md,LOG_TOTAL_ENERGY,10);
- moldyn_set_log(&md,VISUAL_STEP,10);
+ moldyn_set_log(&md,LOG_TOTAL_ENERGY,200);
+ moldyn_set_log(&md,VISUAL_STEP,200);
/*
* let's do the actual md algorithm now
projects / physik / posic.git / commitdiff