+ /*
+ * dU = - p dV
+ *
+ * => p = - dU/dV
+ *
+ */
+
+ /* store atomic configuration + dimension */
+ store=malloc(moldyn->count*sizeof(t_atom));
+ if(store==NULL) {
+ printf("[moldyn] allocating store mem failed\n");
+ return -1;
+ }
+ memcpy(store,moldyn->atom,moldyn->count*sizeof(t_atom));
+ dim=moldyn->dim;
+
+ /* x1, y1 */
+ sd=0.00001;
+ h=(1.0-sd)*(1.0-sd)*(1.0-sd);
+ su=pow(2.0-h,ONE_THIRD)-1.0;
+ dv=(1.0-h)*moldyn->volume;
+
+ /* scale up dimension and atom positions */
+ scale_dim(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
+ scale_atoms(moldyn,SCALE_UP,su,TRUE,TRUE,TRUE);
+ link_cell_shutdown(moldyn);
+ link_cell_init(moldyn,QUIET);
+ potential_force_calc(moldyn);
+ y1=moldyn->energy;
+
+ /* restore atomic configuration + dim */
+ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
+ moldyn->dim=dim;
+
+ /* scale down dimension and atom positions */
+ scale_dim(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
+ scale_atoms(moldyn,SCALE_DOWN,sd,TRUE,TRUE,TRUE);
+ link_cell_shutdown(moldyn);
+ link_cell_init(moldyn,QUIET);
+ potential_force_calc(moldyn);
+ y0=moldyn->energy;
+
+ /* calculate pressure */
+ moldyn->tp=-(y1-y0)/(2.0*dv);
+
+ /* restore atomic configuration */
+ memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
+ moldyn->dim=dim;
+ link_cell_shutdown(moldyn);
+ link_cell_init(moldyn,QUIET);
+ //potential_force_calc(moldyn);
+
+ /* free store buffer */
+ if(store)
+ free(store);
+
+ return moldyn->tp;
+}
+
+double get_pressure(t_moldyn *moldyn) {
+
+ return moldyn->p;
+
+}
+
+int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
+
+ t_3dvec *dim;
+
+ dim=&(moldyn->dim);
+
+ if(dir==SCALE_UP)
+ scale=1.0+scale;
+
+ if(dir==SCALE_DOWN)
+ scale=1.0-scale;
+
+ if(x) dim->x*=scale;
+ if(y) dim->y*=scale;
+ if(z) dim->z*=scale;
+
+ return 0;
+}
+
+int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
+
+ int i;
+ t_3dvec *r;
+
+ if(dir==SCALE_UP)
+ scale=1.0+scale;
+
+ if(dir==SCALE_DOWN)
+ scale=1.0-scale;
+
+ 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_atoms_ind(t_moldyn *moldyn,double x,double y,double z) {
+
+ int i;
+ t_3dvec *r;
+
+ for(i=0;i<moldyn->count;i++) {
+ r=&(moldyn->atom[i].r);
+ r->x*=x;
+ r->y*=y;
+ r->z*=z;
+ }
+
+ return 0;
+}
+
+int scale_dim_ind(t_moldyn *moldyn,double x,double y,double z) {
+
+ t_3dvec *dim;
+
+ dim=&(moldyn->dim);
+
+ dim->x*=x;
+ dim->y*=y;
+ dim->z*=z;
+
+ return 0;
+}
+
+int scale_volume(t_moldyn *moldyn) {
+
+ t_3dvec *dim,*vdim;
+ double scale;
+ t_linkcell *lc;
+ //double sx,sy,sz;
+
+ vdim=&(moldyn->vis.dim);
+ dim=&(moldyn->dim);
+ lc=&(moldyn->lc);
+
+ /* scaling factor */
+ if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
+ scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau;
+ scale=pow(scale,ONE_THIRD);
+ }
+ else {
+ scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
+ }
+
+
+ /*
+ sx=1.0-(moldyn->p_ref-moldyn->px)*moldyn->p_tc*moldyn->tau;
+ sy=1.0-(moldyn->p_ref-moldyn->py)*moldyn->p_tc*moldyn->tau;
+ sz=1.0-(moldyn->p_ref-moldyn->pz)*moldyn->p_tc*moldyn->tau;
+ sx=pow(sx,ONE_THIRD);
+ sy=pow(sy,ONE_THIRD);
+ sz=pow(sz,ONE_THIRD);
+ */
+
+ /* scale the atoms and dimensions */
+ scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
+ scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
+ //scale_atoms_ind(moldyn,sx,sy,sz);
+ //scale_dim_ind(moldyn,sx,sy,sz);
+
+ /* visualize dimensions */
+ if(vdim->x!=0) {
+ vdim->x=dim->x;
+ vdim->y=dim->y;
+ vdim->z=dim->z;
+ }
+
+ /* 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,QUIET);
+ } else {
+ lc->x*=scale;
+ lc->y*=scale;
+ lc->z*=scale;
+ //lc->x*=sx;
+ //lc->y*=sx;
+ //lc->z*=sy;
+ }
+
+ return 0;
+
+}
+
+double e_kin_calc(t_moldyn *moldyn) {
+
+ int i;
+ t_atom *atom;
+
+ atom=moldyn->atom;
+ moldyn->ekin=0.0;
+ //moldyn->ekinx=0.0;
+ //moldyn->ekiny=0.0;
+ //moldyn->ekinz=0.0;
+
+ for(i=0;i<moldyn->count;i++) {
+ atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+ moldyn->ekin+=atom[i].ekin;
+ //moldyn->ekinx+=0.5*atom[i].mass*atom[i].v.x*atom[i].v.x;
+ //moldyn->ekiny+=0.5*atom[i].mass*atom[i].v.y*atom[i].v.y;
+ //moldyn->ekinz+=0.5*atom[i].mass*atom[i].v.z*atom[i].v.z;
+ }
+
+ return moldyn->ekin;
+}
+
+double get_total_energy(t_moldyn *moldyn) {
+
+ return(moldyn->ekin+moldyn->energy);
+}
+
+t_3dvec get_total_p(t_moldyn *moldyn) {
+
+ t_3dvec p,p_total;
+ int i;
+ t_atom *atom;
+
+ atom=moldyn->atom;
+
+ v3_zero(&p_total);
+ for(i=0;i<moldyn->count;i++) {
+ v3_scale(&p,&(atom[i].v),atom[i].mass);
+ v3_add(&p_total,&p_total,&p);
+ }
+
+ return p_total;
+}
+
+double estimate_time_step(t_moldyn *moldyn,double nn_dist) {