diff --git a/src/SHOCK/fix_msst.cpp b/src/SHOCK/fix_msst.cpp index 6399b10fa20b6a0f408b35379982d0930f9413bf..b66f0ed4d35eb961ff3799ffd819c037e5cd6761 100644 --- a/src/SHOCK/fix_msst.cpp +++ b/src/SHOCK/fix_msst.cpp @@ -42,8 +42,8 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), old_velocity(NULL), rfix(NULL), - id_temp(NULL), id_press(NULL), id_pe(NULL), temperature(NULL), + Fix(lmp, narg, arg), old_velocity(NULL), rfix(NULL), + id_temp(NULL), id_press(NULL), id_pe(NULL), temperature(NULL), pressure(NULL), pe(NULL) { if (narg < 4) error->all(FLERR,"Illegal fix msst command"); @@ -131,7 +131,7 @@ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg],"beta") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command"); beta = force->numeric(FLERR,arg[iarg+1]); - if (beta < 0.0 || beta > 1.0) + if (beta < 0.0 || beta > 1.0) error->all(FLERR,"Illegal fix msst command"); iarg += 2; } else error->all(FLERR,"Illegal fix msst command"); @@ -348,7 +348,7 @@ void FixMSST::init() for (int i = 0; i < modify->nfix; i++) if (strcmp(modify->fix[i]->style,"external") == 0) fix_external = (FixExternal *) modify->fix[i]; - if (fix_external == NULL) + if (fix_external == NULL) error->all(FLERR,"Fix msst dftb cannot be used w/out fix external"); } } @@ -434,6 +434,7 @@ void FixMSST::setup(int vflag) // trigger virial computation on next timestep + pe->addstep(update->ntimestep+1); pressure->addstep(update->ntimestep+1); } @@ -468,7 +469,7 @@ void FixMSST::initial_integrate(int vflag) // must convert energy to mv^2 units if (dftb) { - TS_dftb = fix_external->compute_vector(0); + double TS_dftb = fix_external->compute_vector(0); TS = force->ftm2v*TS_dftb; if (update->ntimestep == 1) T0S0 = TS; } else { @@ -493,7 +494,7 @@ void FixMSST::initial_integrate(int vflag) TS_int += (update->dt*TS_dot); //TS_int += (update->dt*TS_dot)/total_mass; - // compute etot + extra terms for conserved quantity + // compute etot + extra terms for conserved quantity double e_scale = compute_etotal() + compute_scalar(); @@ -531,7 +532,7 @@ void FixMSST::initial_integrate(int vflag) for ( k = 0; k < 3; k++ ) { double C = f[i][k] * force->ftm2v / mass[type[i]]; TS_term = TS_dot/(mass[type[i]]*velocity_sum); - escale_term = force->ftm2v*beta*(e0-e_scale) / + escale_term = force->ftm2v*beta*(e0-e_scale) / (mass[type[i]]*velocity_sum); double D = mu * omega[sd] * omega[sd] / (velocity_sum * mass[type[i]] * vol ); @@ -591,7 +592,7 @@ void FixMSST::initial_integrate(int vflag) for ( k = 0; k < 3; k++ ) { double C = f[i][k] * force->ftm2v / mass[type[i]]; TS_term = TS_dot/(mass[type[i]]*velocity_sum); - escale_term = force->ftm2v*beta*(e0-e_scale) / + escale_term = force->ftm2v*beta*(e0-e_scale) / (mass[type[i]]*velocity_sum); double D = mu * omega[sd] * omega[sd] / (velocity_sum * mass[type[i]] * vol ); @@ -668,7 +669,7 @@ void FixMSST::final_integrate() { int i; double p_msst; // MSST driving pressure - double TS,TS_term,escale_term; + double TS_term,escale_term; // v update only for atoms in MSST group @@ -682,15 +683,7 @@ void FixMSST::final_integrate() int sd = direction; - // for DFTB, extract TS_dftb from fix external - // must convert energy to mv^2 units - - if (dftb) { - TS_dftb = fix_external->compute_vector(0); - TS = force->ftm2v*TS_dftb; - } else TS = 0.0; - - // compute etot + extra terms for conserved quantity + // compute etot + extra terms for conserved quantity double e_scale = compute_etotal() + compute_scalar(); @@ -702,7 +695,7 @@ void FixMSST::final_integrate() for ( int k = 0; k < 3; k++ ) { double C = f[i][k] * force->ftm2v / mass[type[i]]; TS_term = TS_dot/(mass[type[i]]*velocity_sum); - escale_term = force->ftm2v*beta*(e0-e_scale) / + escale_term = force->ftm2v*beta*(e0-e_scale) / (mass[type[i]]*velocity_sum); double D = mu * omega[sd] * omega[sd] / (velocity_sum * mass[type[i]] * vol ); @@ -875,7 +868,7 @@ void FixMSST::restart(char *buf) omega[direction] = list[n++]; e0 = list[n++]; v0 = list[n++]; - p0 = list[n++]; + p0 = list[n++]; TS_int = list[n++]; tscale = 0.0; // set tscale to zero for restart p0_set = 1; @@ -899,7 +892,7 @@ int FixMSST::modify_param(int narg, char **arg) strcpy(id_temp,arg[1]); int icompute = modify->find_compute(id_temp); - if (icompute < 0) + if (icompute < 0) error->all(FLERR,"Could not find fix_modify temperature ID"); temperature = modify->compute[icompute]; diff --git a/src/SHOCK/fix_msst.h b/src/SHOCK/fix_msst.h index 450256bcab58fd088dae7e6564d7d6de1bb2fe26..092017bfbb76f69185abbe4fc33e8f0e74c36de3 100644 --- a/src/SHOCK/fix_msst.h +++ b/src/SHOCK/fix_msst.h @@ -54,12 +54,12 @@ class FixMSST : public Fix { int dftb; // flag for use with DFTB+ double velocity_sum; // Sum of the velocities squared - double damping; // Damping function for TS force term at + double damping; // Damping function for TS force term at // small volume difference (v0 - vol) - double T0S0; // Initial TS term for DFTB+ simulations - double S_elec,S_elec_1,S_elec_2; // time history of electron entropy - // for DFTB+ simulaitons - double TS_dot; // time derivative of TS term for + double T0S0; // Initial TS term for DFTB+ simulations + double S_elec,S_elec_1,S_elec_2; // time history of electron entropy + // for DFTB+ simulaitons + double TS_dot; // time derivative of TS term for // DFTB+ simulations double **old_velocity; // Saved velocities @@ -88,12 +88,11 @@ class FixMSST : public Fix { int p0_set; // Is pressure set int v0_set; // Is volume set int e0_set; // Is energy set - double TS_int; // Needed for conserved quantity + double TS_int; // Needed for conserved quantity // with thermal electronic excitations double beta; // Energy conservation scaling factor int maxold; // allocated size of old_velocity - double TS_dftb; // value needed from DFTB+ via fix external class FixExternal *fix_external; // ptr to fix external // functions diff --git a/src/fix.cpp b/src/fix.cpp index ca9a69606c375e01032a05eaebd602f760ffe77e..e6e8a292b98f70bbc6365c7580ee78abd51b022a 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -161,7 +161,7 @@ void Fix::modify_params(int narg, char **arg) /* ---------------------------------------------------------------------- setup for energy, virial computation see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) - fixes call this if use ev_tally() + fixes call this if they use ev_tally() ------------------------------------------------------------------------- */ void Fix::ev_setup(int eflag, int vflag)