diff --git a/src/EXTRA-FIX/fix_deform_pressure.cpp b/src/EXTRA-FIX/fix_deform_pressure.cpp index c76550b40f..7f3a627b66 100644 --- a/src/EXTRA-FIX/fix_deform_pressure.cpp +++ b/src/EXTRA-FIX/fix_deform_pressure.cpp @@ -50,9 +50,160 @@ enum{NOCOUPLE=0,XYZ,XY,YZ,XZ}; FixDeformPressure::FixDeformPressure(LAMMPS *lmp, int narg, char **arg) : FixDeform(lmp, narg, arg), id_temp(nullptr), id_press(nullptr) { + // set defaults + + set_extra = new SetExtra[7]; + memset(set_extra, 0, 7 * sizeof(SetExtra)); + memset(&set_box, 0, sizeof(Set)); + + // parse child-specific arguments + + int index; + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg], "x") == 0 || + strcmp(arg[iarg], "y") == 0 || + strcmp(arg[iarg], "z") == 0) { + + if (strcmp(arg[iarg], "x") == 0) index = 0; + else if (strcmp(arg[iarg], "y") == 0) index = 1; + else if (strcmp(arg[iarg], "z") == 0) index = 2; + + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure", error); + if (strcmp(arg[iarg + 1], "pressure") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure", error); + set[index].style = PRESSURE; + if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) { + set_extra[index].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + } else { + set_extra[index].pstr = utils::strdup(&arg[iarg + 2][2]); + set_extra[index].pvar_flag = 1; + } + set_extra[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 4; + } else if (strcmp(arg[iarg + 1], "pressure/mean") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure/mean", error); + set[index].style = PMEAN; + if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) { + set_extra[index].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + } else { + set_extra[index].pstr = utils::strdup(&arg[iarg + 2][2]); + set_extra[index].pvar_flag = 1; + } + set_extra[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 4; + } else error->all(FLERR, "Illegal fix deform/pressure command argument: {}", arg[iarg + 1]); + + } else if (strcmp(arg[iarg], "xy") == 0 || + strcmp(arg[iarg], "xz") == 0 || + strcmp(arg[iarg], "yz") == 0) { + + if (strcmp(arg[iarg], "xy") == 0) index = 5; + else if (strcmp(arg[iarg], "xz") == 0) index = 4; + else if (strcmp(arg[iarg], "yz") == 0) index = 3; + + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure", error); + if (strcmp(arg[iarg + 1], "pressure") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure", error); + set[index].style = PRESSURE; + if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) { + set_extra[index].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + } else { + set_extra[index].pstr = utils::strdup(&arg[iarg + 2][2]); + set_extra[index].pvar_flag = 1; + } + set_extra[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 4; + } else error->all(FLERR, "Illegal fix deform/pressure command: {}", arg[iarg + 1]); + } else if (strcmp(arg[iarg], "box") == 0) { + if (strcmp(arg[iarg + 1], "volume") == 0) { + set_box.style = VOLUME; + iarg += 2; + } else if (strcmp(arg[iarg + 1], "pressure") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure", error); + set_box.style = PRESSURE; + if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) { + set_extra[6].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + } else { + set_extra[6].pstr = utils::strdup(&arg[iarg + 2][2]); + set_extra[6].pvar_flag = 1; + } + set_extra[6].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 4; + } else error->all(FLERR, "Illegal fix deform/pressure command argument: {}", arg[iarg + 1]); + } else break; + } + + // read options from end of input line + + options(narg - iarg, &arg[iarg]); + + // repeat: check triclinic + + dimension = domain->dimension; + if (triclinic == 0 && (set[3].style || set[4].style || set[5].style)) + error->all(FLERR, "Fix deform tilt factors require triclinic box"); + + // repeat: setup dimflags used by other classes to check for volume-change conflicts + + for (int i = 0; i < 6; i++) + if (set[i].style == NONE) dimflag[i] = 0; + else dimflag[i] = 1; + + if (set_box.style != NONE) { + dimflag[0] = 1; + dimflag[1] = 1; + dimflag[2] = 1; + } + + if (dimflag[0]) box_change |= BOX_CHANGE_X; + if (dimflag[1]) box_change |= BOX_CHANGE_Y; + if (dimflag[2]) box_change |= BOX_CHANGE_Z; + if (dimflag[3]) box_change |= BOX_CHANGE_YZ; + if (dimflag[4]) box_change |= BOX_CHANGE_XZ; + if (dimflag[5]) box_change |= BOX_CHANGE_XY; + + // repeat: no tensile deformation on shrink-wrapped dims + // b/c shrink wrap will change box-length + + for (int i = 0; i < 3; i++) + if (set_box.style && (domain->boundary[i][0] >= 2 || domain->boundary[i][1] >= 2)) + error->all(FLERR, "Cannot use fix deform on a shrink-wrapped boundary"); + + // repeat: no tilt deformation on shrink-wrapped 2nd dim + // b/c shrink wrap will change tilt factor in domain::reset_box() + + if (set[3].style && (domain->boundary[2][0] >= 2 || domain->boundary[2][1] >= 2)) + error->all(FLERR, "Cannot use fix deform tilt on a shrink-wrapped 2nd dim"); + if (set[4].style && (domain->boundary[2][0] >= 2 || domain->boundary[2][1] >= 2)) + error->all(FLERR, "Cannot use fix deform tilt on a shrink-wrapped 2nd dim"); + if (set[5].style && (domain->boundary[1][0] >= 2 || domain->boundary[1][1] >= 2)) + error->all(FLERR, "Cannot use fix deform tilt on a shrink-wrapped 2nd dim"); + + // repeat: set varflag + + for (int i = 0; i < 7; i++) + if (set_extra[i].pvar_flag) varflag = 1; + + // repeat: reneighboring only forced if flips can occur due to shape changes + + if ((!force_reneighbor) && flipflag && (set[3].style || set[4].style || set[5].style)) { + force_reneighbor = 1; + irregular = new Irregular(lmp); + } + + // set initial values at time fix deform is issued + + set_box.vol_initial = domain->xprd * domain->yprd * domain->zprd; + // populate coupled pressure controls if (pcouple != NOCOUPLE) { + + if (dimension == 2) + if (pcouple == XYZ || pcouple == XZ || pcouple == YZ) + error->all(FLERR, "Cannot couple Z dimension in fix deform/pressure in 2D"); + int coupled_indices[3] = {0}; int j = -1; double couple_gain, coupled_pressure; @@ -68,7 +219,7 @@ id_temp(nullptr), id_press(nullptr) // Check coupled styles and find reference for (int i = 0; i < 3; i++) { if (coupled_indices[i]) { - set[i].coupled_flag = 1; + set_extra[i].coupled_flag = 1; if (set[i].style != PRESSURE && set[i].style != NONE) error->all(FLERR, "Cannot couple non-pressure-controlled dimensions"); if (set[i].style == PRESSURE) @@ -86,12 +237,12 @@ id_temp(nullptr), id_press(nullptr) if (set[i].style == NONE) { set[i].style = PRESSURE; dimflag[i] = 1; - set[i].pgain = set[j].pgain; - if (set[j].pvar_flag) { - set[i].pstr = set[j].pstr; - set[i].pvar_flag = 1; + set_extra[i].pgain = set_extra[j].pgain; + if (set_extra[j].pvar_flag) { + set_extra[i].pstr = set_extra[j].pstr; + set_extra[i].pvar_flag = 1; } else { - set[i].ptarget = set[j].ptarget; + set_extra[i].ptarget = set_extra[j].ptarget; } } else { // Check for incompatibilities in style @@ -100,31 +251,18 @@ id_temp(nullptr), id_press(nullptr) if (set[j].style != PRESSURE) continue; // If pressure controlled, check for incompatibilities in parameters - if (set[i].pgain != set[j].pgain || set[i].pvar_flag != set[j].pvar_flag || - set[i].ptarget != set[j].ptarget) + if (set_extra[i].pgain != set_extra[j].pgain || set_extra[i].pvar_flag != set_extra[j].pvar_flag || + set_extra[i].ptarget != set_extra[j].ptarget) error->all(FLERR, "Coupled dimensions must have identical gain parameters"); - if (set[j].pvar_flag) - if (strcmp(set[i].pstr, set[j].pstr) != 0) + if (set_extra[j].pvar_flag) + if (strcmp(set_extra[i].pstr, set_extra[j].pstr) != 0) error->all(FLERR, "Coupled dimensions must have the same target pressure"); } } } } - // repeat some checks in child class to catch changes to pcouple - - if (dimflag[0]) box_change |= BOX_CHANGE_X; - if (dimflag[1]) box_change |= BOX_CHANGE_Y; - if (dimflag[2]) box_change |= BOX_CHANGE_Z; - - // no tensile deformation on shrink-wrapped dims - // b/c shrink wrap will change box-length - - for (int i = 0; i < 3; i++) - if ((set[i].style || set[6].style) && (domain->boundary[i][0] >= 2 || domain->boundary[i][1] >= 2)) - error->all(FLERR, "Cannot use fix deform/pressure on a shrink-wrapped boundary"); - // if vol/balance/p used, must have 2 free dimensions if (vol_balance_flag) { @@ -146,10 +284,15 @@ id_temp(nullptr), id_press(nullptr) // set pressure_flag pressure_flag = 0; - for (int i = 0; i < 7; i++) { - if (set[i].style == PRESSURE || set[i].style == PMEAN) pressure_flag = 1; - if (set[i].coupled_flag) pressure_flag = 1; + for (int i = 0; i < 6; i++) { + if (set[i].style == PRESSURE || set[i].style == PMEAN) { + pressure_flag = 1; + if (set_extra[i].pgain <= 0.0) + error->all(FLERR, "Illegal fix deform/pressure gain constant, must be positive"); + } + if (set_extra[i].coupled_flag) pressure_flag = 1; } + if (set_box.style == PRESSURE) pressure_flag = 1; if (vol_balance_flag) pressure_flag = 1; // check conflict between constant volume/pressure @@ -164,12 +307,12 @@ id_temp(nullptr), id_press(nullptr) if (set[i].style == PMEAN) error->all(FLERR, "Cannot use fix deform/pressure to assign constant volume and pressure"); - // check conflicts between x,y,z styles and iso + // check conflicts between x,y,z styles and box - if (set[6].style) + if (set_box.style) for (int i = 0; i < 3; i++) if (set[i].style == FINAL || set[i].style == DELTA || set[i].style == SCALE || set[i].style == PMEAN || set[i].style == VARIABLE) - error->all(FLERR, "Cannot use fix deform/pressure iso parameter with x, y, or z styles other than vel, erate, trate, pressure, and wiggle"); + error->all(FLERR, "Cannot use fix deform/pressure box parameter with x, y, or z styles other than vel, erate, trate, pressure, and wiggle"); // check pressure used for max rate and normalize error flag @@ -207,6 +350,14 @@ id_temp(nullptr), id_press(nullptr) FixDeformPressure::~FixDeformPressure() { + if (set_extra) + for (int i = 0; i < 7; i++) + delete[] set_extra[i].pstr; + delete[] set_extra; + + delete[] set_box.hstr; + delete[] set_box.hratestr; + // delete temperature and pressure if fix created them if (tflag) modify->delete_compute(id_temp); @@ -221,15 +372,17 @@ void FixDeformPressure::init() { FixDeform::init(); + set_box.vol_start = domain->xprd * domain->yprd * domain->zprd; + // check optional variables for PRESSURE or PMEAN style for (int i = 0; i < 7; i++) { - if (!set[i].pvar_flag) continue; - set[i].pvar = input->variable->find(set[i].pstr); - if (set[i].pvar < 0) - error->all(FLERR, "Variable name {} for fix deform/pressure does not exist", set[i].pstr); - if (!input->variable->equalstyle(set[i].pvar)) - error->all(FLERR, "Variable {} for fix deform/pressure is invalid style", set[i].pstr); + if (!set_extra[i].pvar_flag) continue; + set_extra[i].pvar = input->variable->find(set_extra[i].pstr); + if (set_extra[i].pvar < 0) + error->all(FLERR, "Variable name {} for fix deform/pressure does not exist", set_extra[i].pstr); + if (!input->variable->equalstyle(set_extra[i].pvar)) + error->all(FLERR, "Variable {} for fix deform/pressure is invalid style", set_extra[i].pstr); } // Find pressure/temp computes if needed @@ -265,7 +418,7 @@ void FixDeformPressure::end_of_step() // set new box size for strain-based dims - if (strain_flag) FixDeform::set_strain(); + if (strain_flag) FixDeform::apply_strain(); // set new box size for pressure-based dims @@ -274,30 +427,30 @@ void FixDeformPressure::end_of_step() pressure->compute_vector(); pressure->compute_scalar(); for (int i = 0; i < 3; i++) { - if (!set[i].saved) { - set[i].saved = 1; - set[i].prior_rate = 0.0; - set[i].prior_pressure = pressure->vector[i]; + if (!set_extra[i].saved) { + set_extra[i].saved = 1; + set_extra[i].prior_rate = 0.0; + set_extra[i].prior_pressure = pressure->vector[i]; } } - set_pressure(); + apply_pressure(); } // set new box size for VOLUME dims that are linked to other dims // NOTE: still need to set h_rate for these dims - if (volume_flag) set_volume(); + if (volume_flag) apply_volume(); - // apply any final isotropic scalings + // apply any final box scalings - if (set[6].style) set_iso(); + if (set_box.style) apply_box(); // Save pressure/strain rate if required if (pressure_flag) { for (int i = 0; i < 3; i++) { - set[i].prior_pressure = pressure->vector[i]; - set[i].prior_rate = ((set[i].hi_target - set[i].lo_target) / + set_extra[i].prior_pressure = pressure->vector[i]; + set_extra[i].prior_rate = ((set[i].hi_target - set[i].lo_target) / (domain->boxhi[i] - domain->boxlo[i]) - 1.0) / update->dt; } } @@ -305,7 +458,7 @@ void FixDeformPressure::end_of_step() if (varflag) modify->addstep_compute(update->ntimestep + nevery); - FixDeform::apply_deformation(); + FixDeform::update_domain(); // trigger virial computation, if needed, on next timestep @@ -314,16 +467,16 @@ void FixDeformPressure::end_of_step() } /* ---------------------------------------------------------------------- - set box size for pressure-based dimensions + apply pressure controls ------------------------------------------------------------------------- */ -void FixDeformPressure::set_pressure() +void FixDeformPressure::apply_pressure() { // If variable pressure, calculate current target for (int i = 0; i < 6; i++) if (set[i].style == PRESSURE) - if (set[i].pvar_flag) - set[i].ptarget = input->variable->compute_equal(set[i].pvar); + if (set_extra[i].pvar_flag) + set_extra[i].ptarget = input->variable->compute_equal(set_extra[i].pvar); // Find current (possibly coupled/hydrostatic) pressure for X, Y, Z double *tensor = pressure->vector; @@ -359,14 +512,14 @@ void FixDeformPressure::set_pressure() for (int i = 0; i < 3; i++) { if (set[i].style != PRESSURE && set[i].style != PMEAN) continue; - h_rate[i] = set[i].pgain * (p_current[i] - set[i].ptarget); + h_rate[i] = set_extra[i].pgain * (p_current[i] - set_extra[i].ptarget); if (normalize_pressure_flag) { - if (set[i].ptarget == 0) { + if (set_extra[i].ptarget == 0) { if (max_h_rate == 0) { error->all(FLERR, "Cannot normalize error for zero pressure without defining a max rate"); } else h_rate[i] = max_h_rate * h_rate[i] / fabs(h_rate[i]); - } else h_rate[i] /= fabs(set[i].ptarget); + } else h_rate[i] /= fabs(set_extra[i].ptarget); } if (max_h_rate != 0) @@ -398,13 +551,13 @@ void FixDeformPressure::set_pressure() pcurrent = tensor[3]; } - h_rate[i] = L * set[i].pgain * (pcurrent - set[i].ptarget); + h_rate[i] = L * set_extra[i].pgain * (pcurrent - set_extra[i].ptarget); if (normalize_pressure_flag) { - if (set[i].ptarget == 0) { + if (set_extra[i].ptarget == 0) { if (max_h_rate == 0) { error->all(FLERR, "Cannot normalize error for zero pressure without defining a max rate"); } else h_rate[i] = max_h_rate * h_rate[i] / fabs(h_rate[i]); - } else h_rate[i] /= fabs(set[i].ptarget); + } else h_rate[i] /= fabs(set_extra[i].ptarget); } if (max_h_rate != 0) @@ -416,10 +569,10 @@ void FixDeformPressure::set_pressure() } /* ---------------------------------------------------------------------- - set box size for VOLUME dimensions + apply volume controls ------------------------------------------------------------------------- */ -void FixDeformPressure::set_volume() +void FixDeformPressure::apply_volume() { double e1, e2; int linked_pressure = 0; @@ -446,8 +599,8 @@ void FixDeformPressure::set_volume() (set[fixed].hi_start - set[fixed].lo_start)); } else { double dt = update->dt; - double e1i = set[i].prior_rate; - double e2i = set[fixed].prior_rate; + double e1i = set_extra[i].prior_rate; + double e2i = set_extra[fixed].prior_rate; double L1i = domain->boxhi[i] - domain->boxlo[i]; double L2i = domain->boxhi[fixed] - domain->boxlo[fixed]; double L3i = domain->boxhi[dynamic1] - domain->boxlo[dynamic1]; @@ -457,8 +610,8 @@ void FixDeformPressure::set_volume() double e3 = (L3 / L3i - 1.0) / dt; double p1 = pressure->vector[i]; double p2 = pressure->vector[fixed]; - double p1i = set[i].prior_pressure; - double p2i = set[fixed].prior_pressure; + double p1i = set_extra[i].prior_pressure; + double p2i = set_extra[fixed].prior_pressure; double denominator; if (e3 == 0) { @@ -539,17 +692,17 @@ void FixDeformPressure::adjust_linked_rates(double &e_larger, double &e_smaller, } /* ---------------------------------------------------------------------- - apply isotropic controls + apply box controls ------------------------------------------------------------------------- */ -void FixDeformPressure::set_iso() +void FixDeformPressure::apply_box() { int i; double scale, shift; double v_rate; - if (set[6].style == VOLUME) { - double v0 = set[6].vol_start; + if (set_box.style == VOLUME) { + double v0 = set_box.vol_start; double v = 1.0; for (i = 0; i < 3; i++) v *= (set[i].hi_target - set[i].lo_target); @@ -566,28 +719,28 @@ void FixDeformPressure::set_iso() h_ratelo[i] = -0.5 * h_rate[i]; } - } else if (set[6].style == PRESSURE) { + } else if (set_box.style == PRESSURE) { // If variable pressure, calculate current target - if (set[6].pvar_flag) - set[6].ptarget = input->variable->compute_equal(set[6].pvar); + if (set_extra[6].pvar_flag) + set_extra[6].ptarget = input->variable->compute_equal(set_extra[6].pvar); - v_rate = set[6].pgain * (pressure->scalar- set[6].ptarget); + v_rate = set_extra[6].pgain * (pressure->scalar - set_extra[6].ptarget); if (normalize_pressure_flag) { - if (set[6].ptarget == 0) { + if (set_extra[6].ptarget == 0) { if (max_h_rate == 0) { error->all(FLERR, "Cannot normalize error for zero pressure without defining a max rate"); } else v_rate = max_h_rate * v_rate / fabs(v_rate); - } else v_rate /= fabs(set[6].ptarget); + } else v_rate /= fabs(set_extra[6].ptarget); } if (max_h_rate != 0) if (fabs(v_rate) > max_h_rate) v_rate = max_h_rate * v_rate / fabs(v_rate); - set[6].cumulative_strain += update->dt * v_rate; - scale = (1.0 + set[6].cumulative_strain); + set_extra[6].cumulative_strain += update->dt * v_rate; + scale = (1.0 + set_extra[6].cumulative_strain); for (i = 0; i < 3; i++) { shift = 0.5 * (set[i].hi_target - set[i].lo_target) * scale; set[i].lo_target = 0.5 * (set[i].lo_start + set[i].hi_start) - shift; @@ -601,235 +754,95 @@ void FixDeformPressure::set_iso() } } +/* ---------------------------------------------------------------------- + write Set data to restart file +------------------------------------------------------------------------- */ + +void FixDeformPressure::write_restart(FILE *fp) +{ + if (comm->me == 0) { + int size = 9 * sizeof(double) + 7 * sizeof(Set) + 7 * sizeof(SetExtra); + fwrite(&size, sizeof(int), 1, fp); + fwrite(h_rate, sizeof(double), 6, fp); + fwrite(h_ratelo, sizeof(double), 3, fp); + fwrite(set, sizeof(Set), 6, fp); + fwrite(&set_box, sizeof(Set), 1, fp); + fwrite(set_extra, sizeof(SetExtra), 7, fp); + } +} + +/* ---------------------------------------------------------------------- + use selected state info from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixDeformPressure::restart(char *buf) +{ + int n = 0; + auto list = (double *) buf; + for (int i = 0; i < 6; i++) + h_rate[i] = list[n++]; + for (int i = 0; i < 3; i++) + h_ratelo[i] = list[n++]; + + n = n * sizeof(double); + int samestyle = 1; + Set *set_restart = (Set *) &buf[n]; + for (int i = 0; i < 6; ++i) { + // restore data from initial state + set[i].lo_initial = set_restart[i].lo_initial; + set[i].hi_initial = set_restart[i].hi_initial; + set[i].vol_initial = set_restart[i].vol_initial; + set[i].tilt_initial = set_restart[i].tilt_initial; + // check if style settings are consistent (should do the whole set?) + if (set[i].style != set_restart[i].style) + samestyle = 0; + if (set[i].substyle != set_restart[i].substyle) + samestyle = 0; + } + n += 6 * sizeof(Set); + + // Only restore relevant box variables & check consistency + Set set_box_restart; + memcpy(&set_box_restart, (Set *) &buf[n], sizeof(Set)); + set_box.vol_initial = set_box_restart.vol_initial; + if (set_box.style != set_box_restart.style) + samestyle = 0; + + if (!samestyle) + error->all(FLERR, "Fix deform/pressure settings not consistent with restart"); + + n += sizeof(Set); + SetExtra *set_extra_restart = (SetExtra *) &buf[n]; + for (int i = 0; i < 7; ++i) { + set_extra[i].saved = set_extra_restart[i].saved; + set_extra[i].prior_rate = set_extra_restart[i].prior_rate; + set_extra[i].prior_pressure = set_extra_restart[i].prior_pressure; + set_extra[i].cumulative_strain = set_extra_restart[i].cumulative_strain; + } +} + + /* ---------------------------------------------------------------------- */ void FixDeformPressure::options(int narg, char **arg) { if (narg < 0) error->all(FLERR, "Illegal fix deform/pressure command"); - remapflag = Domain::X_REMAP; - scaleflag = 1; - flipflag = 1; - pcouple = NOCOUPLE; - dimension = domain->dimension; max_h_rate = 0.0; vol_balance_flag = 0; normalize_pressure_flag = 0; - int index; - int iarg = 4; + int iarg = 0; while (iarg < narg) { - if (strcmp(arg[iarg], "x") == 0 || - strcmp(arg[iarg], "y") == 0 || - strcmp(arg[iarg], "z") == 0) { - - if (strcmp(arg[iarg], "x") == 0) index = 0; - else if (strcmp(arg[iarg], "y") == 0) index = 1; - else if (strcmp(arg[iarg], "z") == 0) index = 2; - - if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure", error); - if (strcmp(arg[iarg + 1], "final") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure final", error); - set[index].style = FINAL; - set[index].flo = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - set[index].fhi = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "delta") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure delta", error); - set[index].style = DELTA; - set[index].dlo = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - set[index].dhi = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "scale") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure scale", error); - set[index].style = SCALE; - set[index].scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "vel") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure vel", error); - set[index].style = VEL; - set[index].vel = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "erate") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure erate", error); - set[index].style = ERATE; - set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "trate") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure trate", error); - set[index].style = TRATE; - set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "volume") == 0) { - set[index].style = VOLUME; - iarg += 2; - } else if (strcmp(arg[iarg + 1], "wiggle") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure wiggle", error); - set[index].style = WIGGLE; - set[index].amplitude = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - set[index].tperiod = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - if (set[index].tperiod <= 0.0) - error->all(FLERR, "Illegal fix deform/pressure wiggle period, must be positive"); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "variable") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure variable", error); - set[index].style = VARIABLE; - if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) - error->all(FLERR, "Illegal fix deform/pressure variable name {}", arg[iarg + 2]); - if (strstr(arg[iarg + 3], "v_") != arg[iarg + 3]) - error->all(FLERR, "Illegal fix deform/pressure variable name {}", arg[iarg + 3]); - delete[] set[index].hstr; - delete[] set[index].hratestr; - set[index].hstr = utils::strdup(&arg[iarg + 2][2]); - set[index].hratestr = utils::strdup(&arg[iarg + 3][2]); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "pressure") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure", error); - set[index].style = PRESSURE; - if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) { - set[index].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - } else { - set[index].pstr = utils::strdup(&arg[iarg + 2][2]); - set[index].pvar_flag = 1; - } - set[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - if (set[index].pgain <= 0.0) - error->all(FLERR, "Illegal fix deform/pressure pressure gain, must be positive"); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "pressure/mean") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure/mean", error); - set[index].style = PMEAN; - if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) { - set[index].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - } else { - set[index].pstr = utils::strdup(&arg[iarg + 2][2]); - set[index].pvar_flag = 1; - } - set[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - if (set[index].pgain <= 0.0) - error->all(FLERR, "Illegal fix deform/pressure pressure gain, must be positive"); - iarg += 4; - } else error->all(FLERR, "Illegal fix deform/pressure command argument: {}", arg[iarg + 1]); - - } else if (strcmp(arg[iarg], "xy") == 0 || - strcmp(arg[iarg], "xz") == 0 || - strcmp(arg[iarg], "yz") == 0) { - - if (triclinic == 0) - error->all(FLERR, "fix deform/pressure tilt factors require triclinic box"); - if (strcmp(arg[iarg], "xy") == 0) index = 5; - else if (strcmp(arg[iarg], "xz") == 0) index = 4; - else if (strcmp(arg[iarg], "yz") == 0) index = 3; - - if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure", error); - if (strcmp(arg[iarg + 1], "final") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure final", error); - set[index].style = FINAL; - set[index].ftilt = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "delta") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure delta", error); - set[index].style = DELTA; - set[index].dtilt = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "vel") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure vel", error); - set[index].style = VEL; - set[index].vel = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "erate") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure erate", error); - set[index].style = ERATE; - set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "trate") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure trate", error); - set[index].style = TRATE; - set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "wiggle") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure wiggle", error); - set[index].style = WIGGLE; - set[index].amplitude = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - set[index].tperiod = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - if (set[index].tperiod <= 0.0) - error->all(FLERR, "Illegal fix deform/pressure wiggle period, must be positive"); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "variable") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure variable", error); - set[index].style = VARIABLE; - if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) - error->all(FLERR, "Illegal fix deform/pressure variable name {}", arg[iarg + 2]); - if (strstr(arg[iarg + 3], "v_") != arg[iarg + 3]) - error->all(FLERR, "Illegal fix deform/pressure variable name {}", arg[iarg + 3]); - delete[] set[index].hstr; - delete[] set[index].hratestr; - set[index].hstr = utils::strdup(&arg[iarg + 2][2]); - set[index].hratestr = utils::strdup(&arg[iarg + 3][2]); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "pressure") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure", error); - set[index].style = PRESSURE; - if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) { - set[index].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - } else { - set[index].pstr = utils::strdup(&arg[iarg + 2][2]); - set[index].pvar_flag = 1; - } - set[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - if (set[index].pgain <= 0.0) - error->all(FLERR, "Illegal fix deform/pressure pressure gain, must be positive"); - iarg += 4; - } else error->all(FLERR, "Illegal fix deform/pressure command: {}", arg[iarg + 1]); - } else if (strcmp(arg[iarg], "iso") == 0) { - index = 6; - if (strcmp(arg[iarg + 1], "volume") == 0) { - set[index].style = VOLUME; - iarg += 2; - } else if (strcmp(arg[iarg + 1], "pressure") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure pressure", error); - set[index].style = PRESSURE; - if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) { - set[index].ptarget = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - } else { - set[index].pstr = utils::strdup(&arg[iarg + 2][2]); - set[index].pvar_flag = 1; - } - set[index].pgain = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - if (set[index].pgain <= 0.0) - error->all(FLERR, "Illegal fix deform/pressure pressure gain, must be positive"); - iarg += 4; - } - } else break; - } - - while (iarg < narg) { - if (strcmp(arg[iarg], "remap") == 0) { - if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure remap", error); - if (strcmp(arg[iarg + 1], "x") == 0) remapflag = Domain::X_REMAP; - else if (strcmp(arg[iarg + 1], "v") == 0) remapflag = Domain::V_REMAP; - else if (strcmp(arg[iarg + 1], "none") == 0) remapflag = Domain::NO_REMAP; - else error->all(FLERR, "Illegal fix deform/pressure remap command: {}", arg[iarg + 1]); - iarg += 2; - } else if (strcmp(arg[iarg], "units") == 0) { - if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure units", error); - if (strcmp(arg[iarg + 1], "box") == 0) scaleflag = 0; - else if (strcmp(arg[iarg + 1], "lattice") == 0) scaleflag = 1; - else error->all(FLERR, "Illegal fix deform/pressure units command: {}", arg[iarg + 1]); - iarg += 2; - } else if (strcmp(arg[iarg], "flip") == 0) { - if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure flip", error); - flipflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); - iarg += 2; - } else if (strcmp(arg[iarg], "couple") == 0) { + if (strcmp(arg[iarg], "couple") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure couple", error); if (strcmp(arg[iarg + 1], "xyz") == 0) pcouple = XYZ; else if (strcmp(arg[iarg + 1], "xy") == 0) pcouple = XY; else if (strcmp(arg[iarg + 1], "yz") == 0) pcouple = YZ; else if (strcmp(arg[iarg + 1], "xz") == 0) pcouple = XZ; else if (strcmp(arg[iarg + 1], "none") == 0) pcouple = NOCOUPLE; - else error->all(FLERR, "Illegal fix fix deform/pressure couple command: {}", arg[iarg + 1]); + else error->all(FLERR, "Illegal fix deform/pressure couple command: {}", arg[iarg + 1]); iarg += 2; } else if (strcmp(arg[iarg], "max/rate") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform/pressure max/rate", error); @@ -847,10 +860,6 @@ void FixDeformPressure::options(int narg, char **arg) iarg += 2; } else error->all(FLERR, "Illegal fix deform/pressure command: {}", arg[iarg]); } - - if (dimension == 2) - if (pcouple == XYZ || pcouple == XZ || pcouple == YZ) - error->all(FLERR, "Cannot couple Z dimension in fix deform/pressure in 2D"); } /* ---------------------------------------------------------------------- */ diff --git a/src/EXTRA-FIX/fix_deform_pressure.h b/src/EXTRA-FIX/fix_deform_pressure.h index 7e4ad6e35a..41b93f6e98 100644 --- a/src/EXTRA-FIX/fix_deform_pressure.h +++ b/src/EXTRA-FIX/fix_deform_pressure.h @@ -31,6 +31,8 @@ class FixDeformPressure : public FixDeform { void init() override; void setup(int) override; void end_of_step() override; + void write_restart(FILE *) override; + void restart(char *buf) override; int modify_param(int, char **) override; protected: @@ -46,10 +48,22 @@ class FixDeformPressure : public FixDeform { class Compute *temperature, *pressure; int tflag, pflag; + struct SetExtra { + double ptarget, pgain; + double prior_pressure, prior_rate; + double cumulative_strain; + int saved; + char *pstr; + int pvar, pvar_flag; + int coupled_flag; + }; + SetExtra *set_extra; + Set set_box; + void options(int, char **); - void set_pressure(); - void set_volume(); - void set_iso(); + void apply_volume() override; + void apply_pressure(); + void apply_box(); void couple(); void adjust_linked_rates(double&, double&, double, double, double); }; diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index cf45a12f5d..ea68574bf2 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -34,6 +34,8 @@ #include #include +#include +#include using namespace LAMMPS_NS; using namespace FixConst; @@ -56,15 +58,174 @@ irregular(nullptr), set(nullptr) nevery = utils::inumeric(FLERR, arg[3], false, lmp); if (nevery <= 0) error->all(FLERR, "Illegal fix deform command"); + // arguments for child classes + + std::unordered_set child_parameters; + std::unordered_map child_styles; + int nskip; + if (utils::strmatch(style, "pressure$")) { + child_parameters.insert("box"); + child_styles.insert({{"pressure", 4}, {"pressure/mean", 4}, {"volume", 3}}); + } + // set defaults - set = new Set[7]; - memset(set, 0, 7 * sizeof(Set)); + set = new Set[6]; + memset(set, 0, 6 * sizeof(Set)); // parse arguments + int index; + int iarg = 4; + while (iarg < narg) { + if (strcmp(arg[iarg], "x") == 0 || + strcmp(arg[iarg], "y") == 0 || + strcmp(arg[iarg], "z") == 0) { + + if (strcmp(arg[iarg], "x") == 0) index = 0; + else if (strcmp(arg[iarg], "y") == 0) index = 1; + else if (strcmp(arg[iarg], "z") == 0) index = 2; + + if (strcmp(arg[iarg + 1], "final") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform final", error); + set[index].style = FINAL; + set[index].flo = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + set[index].fhi = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 4; + } else if (strcmp(arg[iarg + 1], "delta") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform delta", error); + set[index].style = DELTA; + set[index].dlo = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + set[index].dhi = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + iarg += 4; + } else if (strcmp(arg[iarg + 1], "scale") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform scale", error); + set[index].style = SCALE; + set[index].scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg + 1], "vel") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform vel", error); + set[index].style = VEL; + set[index].vel = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg + 1], "erate") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform erate", error); + set[index].style = ERATE; + set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg + 1], "trate") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform trate", error); + set[index].style = TRATE; + set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg + 1], "volume") == 0) { + set[index].style = VOLUME; + iarg += 2; + } else if (strcmp(arg[iarg + 1], "wiggle") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform wiggle", error); + set[index].style = WIGGLE; + set[index].amplitude = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + set[index].tperiod = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + if (set[index].tperiod <= 0.0) + error->all(FLERR, "Illegal fix deform wiggle period, must be positive"); + iarg += 4; + } else if (strcmp(arg[iarg + 1], "variable") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform variable", error); + set[index].style = VARIABLE; + if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) + error->all(FLERR, "Illegal fix deform variable name {}", arg[iarg + 2]); + if (strstr(arg[iarg + 3], "v_") != arg[iarg + 3]) + error->all(FLERR, "Illegal fix deform variable name {}", arg[iarg + 3]); + delete[] set[index].hstr; + delete[] set[index].hratestr; + set[index].hstr = utils::strdup(&arg[iarg + 2][2]); + set[index].hratestr = utils::strdup(&arg[iarg + 3][2]); + iarg += 4; + } else if (child_styles.find(arg[iarg + 1]) != child_styles.end()) { + nskip = child_styles[arg[iarg + 1]]; + if (iarg + nskip > narg) + utils::missing_cmd_args(FLERR, fmt::format("fix {} {}", style, arg[iarg + 1]), error); + iarg += nskip; + } error->all(FLERR, "Illegal fix deform command argument: {}", arg[iarg + 1]); + + } else if (strcmp(arg[iarg], "xy") == 0 || + strcmp(arg[iarg], "xz") == 0 || + strcmp(arg[iarg], "yz") == 0) { + + if (strcmp(arg[iarg], "xy") == 0) index = 5; + else if (strcmp(arg[iarg], "xz") == 0) index = 4; + else if (strcmp(arg[iarg], "yz") == 0) index = 3; + + if (strcmp(arg[iarg + 1], "final") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform final", error); + set[index].style = FINAL; + set[index].ftilt = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg + 1], "delta") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform delta", error); + set[index].style = DELTA; + set[index].dtilt = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg + 1], "vel") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform vel", error); + set[index].style = VEL; + set[index].vel = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg + 1], "erate") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform erate", error); + set[index].style = ERATE; + set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg + 1], "trate") == 0) { + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform trate", error); + set[index].style = TRATE; + set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + iarg += 3; + } else if (strcmp(arg[iarg + 1], "wiggle") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform wiggle", error); + set[index].style = WIGGLE; + set[index].amplitude = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + set[index].tperiod = utils::numeric(FLERR, arg[iarg + 3], false, lmp); + if (set[index].tperiod <= 0.0) + error->all(FLERR, "Illegal fix deform wiggle period, must be positive"); + iarg += 4; + } else if (strcmp(arg[iarg + 1], "variable") == 0) { + if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform variable", error); + set[index].style = VARIABLE; + if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) + error->all(FLERR, "Illegal fix deform variable name {}", arg[iarg + 2]); + if (strstr(arg[iarg + 3], "v_") != arg[iarg + 3]) + error->all(FLERR, "Illegal fix deform variable name {}", arg[iarg + 3]); + delete[] set[index].hstr; + delete[] set[index].hratestr; + set[index].hstr = utils::strdup(&arg[iarg + 2][2]); + set[index].hratestr = utils::strdup(&arg[iarg + 3][2]); + iarg += 4; + } else if (child_styles.find(arg[iarg + 1]) != child_styles.end()) { + nskip = child_styles[arg[iarg + 1]]; + if (iarg + nskip > narg) + utils::missing_cmd_args(FLERR, fmt::format("fix {} {}", style, arg[iarg + 1]), error); + iarg += nskip; + } error->all(FLERR, "Illegal fix deform command argument: {}", arg[iarg + 1]); + } else if (child_parameters.find(arg[iarg]) != child_parameters.end()) { + if (child_styles.find(arg[iarg + 1]) != child_styles.end()) { + nskip = child_styles[arg[iarg + 1]]; + if (iarg + nskip > narg) + utils::missing_cmd_args(FLERR, fmt::format("fix {} {}", style, arg[iarg + 1]), error); + iarg += nskip; + } error->all(FLERR, "Illegal fix {} command argument: {}", style, arg[iarg + 1]); + } else break; + } + + // reparse all arguments for optional keywords + + options(narg - 4, &arg[4]); + + // check triclinic + triclinic = domain->triclinic; - options(narg, arg); + if (triclinic == 0 && (set[3].style || set[4].style || set[5].style)) + error->all(FLERR, "Fix deform tilt factors require triclinic box"); // no x remap effectively moves atoms within box, so set restart_pbc @@ -76,12 +237,6 @@ irregular(nullptr), set(nullptr) if (set[i].style == NONE) dimflag[i] = 0; else dimflag[i] = 1; - if (set[6].style != NONE) { - dimflag[0] = 1; - dimflag[1] = 1; - dimflag[2] = 1; - } - if (dimflag[0]) box_change |= BOX_CHANGE_X; if (dimflag[1]) box_change |= BOX_CHANGE_Y; if (dimflag[2]) box_change |= BOX_CHANGE_Z; @@ -93,7 +248,7 @@ irregular(nullptr), set(nullptr) // b/c shrink wrap will change box-length for (int i = 0; i < 3; i++) - if ((set[i].style || set[6].style) && (domain->boundary[i][0] >= 2 || domain->boundary[i][1] >= 2)) + if (set[i].style && (domain->boundary[i][0] >= 2 || domain->boundary[i][1] >= 2)) error->all(FLERR, "Cannot use fix deform on a shrink-wrapped boundary"); // no tilt deformation on shrink-wrapped 2nd dim @@ -113,7 +268,7 @@ irregular(nullptr), set(nullptr) if (set[i].style == FINAL || set[i].style == DELTA || set[i].style == VEL || set[i].style == WIGGLE) flag = 1; - double xscale,yscale,zscale; + double xscale, yscale, zscale; if (flag && scaleflag) { xscale = domain->lattice->xlattice; yscale = domain->lattice->ylattice; @@ -187,10 +342,8 @@ irregular(nullptr), set(nullptr) // set varflag varflag = 0; - for (int i = 0; i < 7; i++) { + for (int i = 0; i < 6; i++) if (set[i].style == VARIABLE) varflag = 1; - if (set[i].pvar_flag) varflag = 1; - } // set initial values at time fix deform is issued @@ -202,7 +355,6 @@ irregular(nullptr), set(nullptr) set[3].tilt_initial = domain->yz; set[4].tilt_initial = domain->xz; set[5].tilt_initial = domain->xy; - set[6].vol_initial = domain->xprd * domain->yprd * domain->zprd; // reneighboring only forced if flips can occur due to shape changes @@ -231,10 +383,9 @@ irregular(nullptr), set(nullptr) FixDeform::~FixDeform() { if (set) { - for (int i = 0; i < 7; i++) { + for (int i = 0; i < 6; i++) { delete[] set[i].hstr; delete[] set[i].hratestr; - delete[] set[i].pstr; } } delete[] set; @@ -282,7 +433,7 @@ void FixDeform::init() // check variables for VARIABLE style - for (int i = 0; i < 7; i++) { + for (int i = 0; i < 6; i++) { if (set[i].style != VARIABLE) continue; set[i].hvar = input->variable->find(set[i].hstr); if (set[i].hvar < 0) @@ -403,8 +554,6 @@ void FixDeform::init() } } - set[6].vol_start = domain->xprd * domain->yprd * domain->zprd; - // if using tilt TRATE, then initial tilt must be non-zero for (int i = 3; i < 6; i++) @@ -520,16 +669,16 @@ void FixDeform::end_of_step() // set new box size for strain-based dims - set_strain(); + apply_strain(); // set new box size for VOLUME dims that are linked to other dims // NOTE: still need to set h_rate for these dims - set_volume(); + apply_volume(); if (varflag) modify->addstep_compute(update->ntimestep + nevery); - apply_deformation(); + update_domain(); // redo KSpace coeffs since box has changed @@ -537,10 +686,10 @@ void FixDeform::end_of_step() } /* ---------------------------------------------------------------------- - set box size for strain-based dimensions + apply strain controls ------------------------------------------------------------------------- */ -void FixDeform::set_strain() +void FixDeform::apply_strain() { // for NONE, target is current box size // for TRATE, set target directly based on current time, also set h_rate @@ -620,10 +769,10 @@ void FixDeform::set_strain() } /* ---------------------------------------------------------------------- - set box size for VOLUME dimensions + apply volume controls ------------------------------------------------------------------------- */ -void FixDeform::set_volume() +void FixDeform::apply_volume() { double e1, e2; @@ -657,10 +806,10 @@ void FixDeform::set_volume() } /* ---------------------------------------------------------------------- - Apply calculated deformation + Update box domain ------------------------------------------------------------------------- */ -void FixDeform::apply_deformation() +void FixDeform::update_domain() { // tilt_target can be large positive or large negative value // add/subtract box lengths until tilt_target is closest to current value @@ -766,22 +915,22 @@ void FixDeform::apply_deformation() // reset global and local box to new size/shape // only if deform fix is controlling the dimension - if (set[0].style || set[6].style) { + if (dimflag[0]) { domain->boxlo[0] = set[0].lo_target; domain->boxhi[0] = set[0].hi_target; } - if (set[1].style || set[6].style) { + if (dimflag[1]) { domain->boxlo[1] = set[1].lo_target; domain->boxhi[1] = set[1].hi_target; } - if (set[2].style || set[6].style) { + if (dimflag[2]) { domain->boxlo[2] = set[2].lo_target; domain->boxhi[2] = set[2].hi_target; } if (triclinic) { - if (set[3].style) domain->yz = set[3].tilt_target; - if (set[4].style) domain->xz = set[4].tilt_target; - if (set[5].style) domain->xy = set[5].tilt_target; + if (dimflag[3]) domain->yz = set[3].tilt_target; + if (dimflag[4]) domain->xz = set[4].tilt_target; + if (dimflag[5]) domain->xy = set[5].tilt_target; } domain->set_global_box(); @@ -810,11 +959,11 @@ void FixDeform::apply_deformation() void FixDeform::write_restart(FILE *fp) { if (comm->me == 0) { - int size = 9 * sizeof(double) + 7 * sizeof(Set); + int size = 9 * sizeof(double) + 6 * sizeof(Set); fwrite(&size, sizeof(int), 1, fp); fwrite(h_rate, sizeof(double), 6, fp); fwrite(h_ratelo, sizeof(double), 3, fp); - fwrite(set, sizeof(Set), 7, fp); + fwrite(set, sizeof(Set), 6, fp); } } @@ -833,16 +982,12 @@ void FixDeform::restart(char *buf) int samestyle = 1; Set *set_restart = (Set *) &buf[n * sizeof(double)]; - for (int i = 0; i < 7; ++i) { + for (int i = 0; i < 6; ++i) { // restore data from initial state set[i].lo_initial = set_restart[i].lo_initial; set[i].hi_initial = set_restart[i].hi_initial; set[i].vol_initial = set_restart[i].vol_initial; set[i].tilt_initial = set_restart[i].tilt_initial; - set[i].saved = set_restart[i].saved; - set[i].prior_rate = set_restart[i].prior_rate; - set[i].prior_pressure = set_restart[i].prior_pressure; - set[i].cumulative_strain = set_restart[i].cumulative_strain; // check if style settings are consistent (should do the whole set?) if (set[i].style != set_restart[i].style) samestyle = 0; @@ -863,135 +1008,14 @@ void FixDeform::options(int narg, char **arg) scaleflag = 1; flipflag = 1; - int index; - int iarg = 4; - while (iarg < narg) { - if (strcmp(arg[iarg], "x") == 0 || - strcmp(arg[iarg], "y") == 0 || - strcmp(arg[iarg], "z") == 0) { + // arguments for child classes - if (strcmp(arg[iarg], "x") == 0) index = 0; - else if (strcmp(arg[iarg], "y") == 0) index = 1; - else if (strcmp(arg[iarg], "z") == 0) index = 2; - - if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform", error); - if (strcmp(arg[iarg + 1], "final") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform final", error); - set[index].style = FINAL; - set[index].flo = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - set[index].fhi = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "delta") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform delta", error); - set[index].style = DELTA; - set[index].dlo = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - set[index].dhi = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "scale") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform scale", error); - set[index].style = SCALE; - set[index].scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "vel") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform vel", error); - set[index].style = VEL; - set[index].vel = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "erate") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform erate", error); - set[index].style = ERATE; - set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "trate") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform trate", error); - set[index].style = TRATE; - set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "volume") == 0) { - set[index].style = VOLUME; - iarg += 2; - } else if (strcmp(arg[iarg + 1], "wiggle") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform wiggle", error); - set[index].style = WIGGLE; - set[index].amplitude = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - set[index].tperiod = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - if (set[index].tperiod <= 0.0) - error->all(FLERR, "Illegal fix deform wiggle period, must be positive"); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "variable") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform variable", error); - set[index].style = VARIABLE; - if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) - error->all(FLERR, "Illegal fix deform variable name {}", arg[iarg + 2]); - if (strstr(arg[iarg + 3], "v_") != arg[iarg + 3]) - error->all(FLERR, "Illegal fix deform variable name {}", arg[iarg + 3]); - delete[] set[index].hstr; - delete[] set[index].hratestr; - set[index].hstr = utils::strdup(&arg[iarg + 2][2]); - set[index].hratestr = utils::strdup(&arg[iarg + 3][2]); - iarg += 4; - } else error->all(FLERR, "Illegal fix deform command argument: {}", arg[iarg + 1]); - - } else if (strcmp(arg[iarg], "xy") == 0 || - strcmp(arg[iarg], "xz") == 0 || - strcmp(arg[iarg], "yz") == 0) { - - if (triclinic == 0) - error->all(FLERR, "Fix deform tilt factors require triclinic box"); - if (strcmp(arg[iarg], "xy") == 0) index = 5; - else if (strcmp(arg[iarg], "xz") == 0) index = 4; - else if (strcmp(arg[iarg], "yz") == 0) index = 3; - - if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform", error); - if (strcmp(arg[iarg + 1], "final") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform final", error); - set[index].style = FINAL; - set[index].ftilt = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "delta") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform delta", error); - set[index].style = DELTA; - set[index].dtilt = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "vel") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform vel", error); - set[index].style = VEL; - set[index].vel = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "erate") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform erate", error); - set[index].style = ERATE; - set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "trate") == 0) { - if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "fix deform trate", error); - set[index].style = TRATE; - set[index].rate = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - iarg += 3; - } else if (strcmp(arg[iarg + 1], "wiggle") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform wiggle", error); - set[index].style = WIGGLE; - set[index].amplitude = utils::numeric(FLERR, arg[iarg + 2], false, lmp); - set[index].tperiod = utils::numeric(FLERR, arg[iarg + 3], false, lmp); - if (set[index].tperiod <= 0.0) - error->all(FLERR, "Illegal fix deform wiggle period, must be positive"); - iarg += 4; - } else if (strcmp(arg[iarg + 1], "variable") == 0) { - if (iarg + 4 > narg) utils::missing_cmd_args(FLERR, "fix deform variable", error); - set[index].style = VARIABLE; - if (strstr(arg[iarg + 2], "v_") != arg[iarg + 2]) - error->all(FLERR, "Illegal fix deform variable name {}", arg[iarg + 2]); - if (strstr(arg[iarg + 3], "v_") != arg[iarg + 3]) - error->all(FLERR, "Illegal fix deform variable name {}", arg[iarg + 3]); - delete[] set[index].hstr; - delete[] set[index].hratestr; - set[index].hstr = utils::strdup(&arg[iarg + 2][2]); - set[index].hratestr = utils::strdup(&arg[iarg + 3][2]); - iarg += 4; - } else error->all(FLERR, "Illegal fix deform command: {}", arg[iarg + 1]); - } else break; - } + std::unordered_map child_options; + int nskip; + if (utils::strmatch(style, "pressure$")) + child_options.insert({{"couple", 2}, {"max/rate", 2}, {"normalize/pressure", 2}, {"vol/balance/p", 2}}); + int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg], "remap") == 0) { if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform remap", error); @@ -1010,6 +1034,11 @@ void FixDeform::options(int narg, char **arg) if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "fix deform flip", error); flipflag = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; + } else if (child_options.find(arg[iarg]) != child_options.end()) { + nskip = child_options[arg[iarg]]; + if (iarg + nskip > narg) + utils::missing_cmd_args(FLERR, fmt::format("fix {} {}", style, arg[iarg]), error); + iarg += nskip; } else error->all(FLERR, "Illegal fix deform command: {}", arg[iarg]); } } diff --git a/src/fix_deform.h b/src/fix_deform.h index 16a563cf4f..69f62fe558 100644 --- a/src/fix_deform.h +++ b/src/fix_deform.h @@ -35,8 +35,8 @@ class FixDeform : public Fix { void init() override; void pre_exchange() override; void end_of_step() override; - void write_restart(FILE *) override; - void restart(char *buf) override; + void virtual write_restart(FILE *) override; + void virtual restart(char *buf) override; double memory_usage() override; protected: @@ -59,22 +59,16 @@ class FixDeform : public Fix { double tilt_initial, tilt_start, tilt_stop, tilt_target, tilt_flip; double tilt_min, tilt_max; double vol_initial, vol_start; - double ptarget, pgain; - double prior_pressure, prior_rate; - double cumulative_strain; - int saved; int fixed, dynamic1, dynamic2; - char *hstr, *hratestr, *pstr; + char *hstr, *hratestr; int hvar, hratevar; - int pvar, pvar_flag; - int coupled_flag; }; Set *set; void options(int, char **); - void set_strain(); - void set_volume(); - void apply_deformation(); + void virtual apply_volume(); + void apply_strain(); + void update_domain(); }; } // namespace LAMMPS_NS