diff --git a/src/BOCS/fix_bocs.cpp b/src/BOCS/fix_bocs.cpp index d35facdc5a..4918f8d879 100644 --- a/src/BOCS/fix_bocs.cpp +++ b/src/BOCS/fix_bocs.cpp @@ -69,12 +69,9 @@ const int NUM_INPUT_DATA_COLUMNS = 2; // columns in the pressure correction ---------------------------------------------------------------------- */ FixBocs::FixBocs(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - rfix(nullptr), id_dilate(nullptr), irregular(nullptr), - id_temp(nullptr), id_press(nullptr), - eta(nullptr), eta_dot(nullptr), eta_dotdot(nullptr), - eta_mass(nullptr), etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr), - etap_mass(nullptr) + Fix(lmp, narg, arg), id_dilate(nullptr), irregular(nullptr), id_temp(nullptr), + id_press(nullptr), eta(nullptr), eta_dot(nullptr), eta_dotdot(nullptr), eta_mass(nullptr), + etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr), etap_mass(nullptr) { if (lmp->citeme) lmp->citeme->add(cite_user_bocs_package); @@ -379,9 +376,6 @@ FixBocs::FixBocs(LAMMPS *lmp, int narg, char **arg) : if (deviatoric_flag) size_vector += 1; } - nrigid = 0; - rfix = nullptr; - if (pre_exchange_flag) irregular = new Irregular(lmp); else irregular = nullptr; @@ -424,31 +418,29 @@ FixBocs::~FixBocs() { if (copymode) return; - delete [] id_dilate; - delete [] rfix; - + delete[] id_dilate; delete irregular; // delete temperature and pressure if fix created them if (tcomputeflag) modify->delete_compute(id_temp); - delete [] id_temp; + delete[] id_temp; if (tstat_flag) { - delete [] eta; - delete [] eta_dot; - delete [] eta_dotdot; - delete [] eta_mass; + delete[] eta; + delete[] eta_dot; + delete[] eta_dotdot; + delete[] eta_mass; } if (pstat_flag) { if (pcomputeflag) modify->delete_compute(id_press); - delete [] id_press; + delete[] id_press; if (mpchain) { - delete [] etap; - delete [] etap_dot; - delete [] etap_dotdot; - delete [] etap_mass; + delete[] etap; + delete[] etap_dot; + delete[] etap_dotdot; + delete[] etap_mass; } } if (p_match_coeffs) free(p_match_coeffs); @@ -596,20 +588,10 @@ void FixBocs::init() } // detect if any rigid fixes exist so rigid bodies move when box is remapped - // rfix[] = indices to each fix rigid - delete [] rfix; - nrigid = 0; - rfix = nullptr; - - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) nrigid++; - if (nrigid) { - rfix = new int[nrigid]; - nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; - } + rfix.clear(); + for (auto &ifix : modify->get_fix_list()) + if (ifix->rigid_flag) rfix.push_back(ifix); } // NJD MRD 2 functions @@ -1204,9 +1186,7 @@ void FixBocs::remap() domain->x2lamda(x[i],x[i]); } - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(0); + for (auto &ifix : rfix) ifix->deform(0); // reset global and local box to new size/shape @@ -1351,9 +1331,7 @@ void FixBocs::remap() domain->lamda2x(x[i],x[i]); } - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(1); + for (auto &ifix : rfix) ifix->deform(1); } /* ---------------------------------------------------------------------- @@ -1512,7 +1490,7 @@ int FixBocs::modify_param(int narg, char **arg) modify->delete_compute(id_temp); tcomputeflag = 0; } - delete [] id_temp; + delete[] id_temp; id_temp = utils::strdup(arg[1]); int icompute = modify->find_compute(arg[1]); @@ -1544,29 +1522,23 @@ int FixBocs::modify_param(int narg, char **arg) modify->delete_compute(id_press); pcomputeflag = 0; } - delete [] id_press; + delete[] id_press; id_press = utils::strdup(arg[1]); - int icompute = modify->find_compute(arg[1]); - if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID"); - pressure = modify->compute[icompute]; + pressure = modify->get_compute_by_id(id_press); + if (!pressure) error->all(FLERR,"Could not find fix_modify pressure ID {}", id_press); + if (!pressure->pressflag) + error->all(FLERR, "Fix_modify pressure ID {} does not compute pressure", id_press); - if (p_match_flag) // NJD MRD - { - if (p_basis_type == BASIS_ANALYTIC) - { - (dynamic_cast(pressure))->send_cg_info(p_basis_type, N_p_match, - p_match_coeffs, N_mol, vavg); + if (p_match_flag) { + auto bocspress = dynamic_cast(pressure); + if (bocspress) { + if (p_basis_type == BASIS_ANALYTIC) { + bocspress->send_cg_info(p_basis_type, N_p_match, p_match_coeffs, N_mol, vavg); + } else if (p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE) { + bocspress->send_cg_info(p_basis_type, splines, spline_length); + } } - else if (p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE ) - { - (dynamic_cast(pressure))->send_cg_info(p_basis_type, splines, spline_length ); - } - } - - if (pressure->pressflag == 0) - { - error->all(FLERR, "Fix_modify pressure ID does not compute pressure"); } return 2; } diff --git a/src/BOCS/fix_bocs.h b/src/BOCS/fix_bocs.h index fd47fda4d7..71fbc273d8 100644 --- a/src/BOCS/fix_bocs.h +++ b/src/BOCS/fix_bocs.h @@ -75,9 +75,8 @@ class FixBocs : public Fix { double drag, tdrag_factor; // drag factor on particle thermostat double pdrag_factor; // drag factor on barostat int kspace_flag; // 1 if KSpace invoked, 0 if not - int nrigid; // number of rigid fixes int dilate_group_bit; // mask for dilation group - int *rfix; // indices of rigid fixes + std::vector rfix; // list of rigid fixes char *id_dilate; // group name to dilate class Irregular *irregular; // for migrating atoms after box flips diff --git a/src/DRUDE/fix_tgnh_drude.cpp b/src/DRUDE/fix_tgnh_drude.cpp index 273f163303..987408fe63 100644 --- a/src/DRUDE/fix_tgnh_drude.cpp +++ b/src/DRUDE/fix_tgnh_drude.cpp @@ -52,14 +52,13 @@ enum{ISO,ANISO,TRICLINIC}; ---------------------------------------------------------------------- */ FixTGNHDrude::FixTGNHDrude(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - rfix(nullptr), irregular(nullptr), id_temp(nullptr), id_press(nullptr), - etamol(nullptr), etamol_dot(nullptr), etamol_dotdot(nullptr), etamol_mass(nullptr), - etaint(nullptr), etaint_dot(nullptr), etaint_dotdot(nullptr), etaint_mass(nullptr), - etadrude(nullptr), etadrude_dot(nullptr), etadrude_dotdot(nullptr), etadrude_mass(nullptr), - etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr), etap_mass(nullptr) + Fix(lmp, narg, arg), irregular(nullptr), id_temp(nullptr), id_press(nullptr), etamol(nullptr), + etamol_dot(nullptr), etamol_dotdot(nullptr), etamol_mass(nullptr), etaint(nullptr), + etaint_dot(nullptr), etaint_dotdot(nullptr), etaint_mass(nullptr), etadrude(nullptr), + etadrude_dot(nullptr), etadrude_dotdot(nullptr), etadrude_mass(nullptr), etap(nullptr), + etap_dot(nullptr), etap_dotdot(nullptr), etap_mass(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal fix nvt/npt/nph command"); + if (narg < 4) error->all(FLERR, "Illegal fix {} command", style); restart_global = 1; dynamic_group_allow = 0; @@ -507,9 +506,6 @@ FixTGNHDrude::FixTGNHDrude(LAMMPS *lmp, int narg, char **arg) : } } - nrigid = 0; - rfix = nullptr; - if (pre_exchange_flag) irregular = new Irregular(lmp); else irregular = nullptr; @@ -519,15 +515,15 @@ FixTGNHDrude::FixTGNHDrude(LAMMPS *lmp, int narg, char **arg) : vol0 = t0 = 0.0; // find fix drude - int ifix; - for (ifix = 0; ifix < modify->nfix; ifix++) - if (strcmp(modify->fix[ifix]->style,"drude") == 0) break; - if (ifix == modify->nfix) error->all(FLERR, "fix tgnh/drude requires fix drude"); - fix_drude = dynamic_cast(modify->fix[ifix]); + + auto fdrude = modify->get_fix_by_style("^drude"); + if (fdrude.size() < 1) error->all(FLERR, "Fix {} requires fix drude", style); + fix_drude = dynamic_cast(fdrude[0]); + if (!fix_drude) error->all(FLERR, "Fix {} requires fix drude", style); // make sure ghost atoms have velocity if (!comm->ghost_velocity) - error->all(FLERR,"fix tgnh/drude requires ghost velocities. Use comm_modify vel yes"); + error->all(FLERR,"Fix {} requires ghost velocities. Use comm_modify vel yes", style); } /* ---------------------------------------------------------------------- */ @@ -536,38 +532,36 @@ FixTGNHDrude::~FixTGNHDrude() { if (copymode) return; - delete [] rfix; - delete irregular; // delete temperature and pressure if fix created them if (tcomputeflag) modify->delete_compute(id_temp); - delete [] id_temp; + delete[] id_temp; if (tstat_flag) { - delete [] etaint; - delete [] etaint_dot; - delete [] etaint_dotdot; - delete [] etaint_mass; - delete [] etamol; - delete [] etamol_dot; - delete [] etamol_dotdot; - delete [] etamol_mass; - delete [] etadrude; - delete [] etadrude_dot; - delete [] etadrude_dotdot; - delete [] etadrude_mass; + delete[] etaint; + delete[] etaint_dot; + delete[] etaint_dotdot; + delete[] etaint_mass; + delete[] etamol; + delete[] etamol_dot; + delete[] etamol_dotdot; + delete[] etamol_mass; + delete[] etadrude; + delete[] etadrude_dot; + delete[] etadrude_dotdot; + delete[] etadrude_mass; } if (pstat_flag) { if (pcomputeflag) modify->delete_compute(id_press); - delete [] id_press; + delete[] id_press; if (mpchain) { - delete [] etap; - delete [] etap_dot; - delete [] etap_dotdot; - delete [] etap_mass; + delete[] etap; + delete[] etap_dot; + delete[] etap_dotdot; + delete[] etap_mass; } } } @@ -605,19 +599,15 @@ void FixTGNHDrude::init() // set temperature and pressure ptrs - int icompute = modify->find_compute(id_temp); - if (icompute < 0) - error->all(FLERR,"Temperature ID for fix nvt/npt does not exist"); - temperature = modify->compute[icompute]; + temperature = modify->get_compute_by_id(id_temp); + if (!temperature) error->all(FLERR,"Temperature ID for fix {} does not exist", style); if (temperature->tempbias) which = BIAS; else which = NOBIAS; if (pstat_flag) { - icompute = modify->find_compute(id_press); - if (icompute < 0) - error->all(FLERR,"Pressure ID for fix npt/nph does not exist"); - pressure = modify->compute[icompute]; + pressure = modify->get_compute_by_id(id_press); + if (!pressure) error->all(FLERR,"Pressure ID for fix {} does not exist", id_press); } // set timesteps and frequencies @@ -670,20 +660,10 @@ void FixTGNHDrude::init() } // detect if any rigid fixes exist so rigid bodies move when box is remapped - // rfix[] = indices to each fix rigid - delete [] rfix; - nrigid = 0; - rfix = nullptr; - - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) nrigid++; - if (nrigid) { - rfix = new int[nrigid]; - nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; - } + rfix.clear(); + for (auto &ifix : modify->get_fix_list()) + if (ifix->rigid_flag) rfix.push_back(ifix); } /* ---------------------------------------------------------------------- @@ -1111,9 +1091,7 @@ void FixTGNHDrude::remap() domain->x2lamda(nlocal); - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(0); + for (auto &ifix : rfix) ifix->deform(0); // reset global and local box to new size/shape @@ -1253,9 +1231,7 @@ void FixTGNHDrude::remap() domain->lamda2x(nlocal); - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(1); + for (auto &ifix : rfix) ifix->deform(1); } /* ---------------------------------------------------------------------- @@ -1426,27 +1402,23 @@ int FixTGNHDrude::modify_param(int narg, char **arg) modify->delete_compute(id_temp); tcomputeflag = 0; } - delete [] id_temp; + delete[] id_temp; id_temp = utils::strdup(arg[1]); - int icompute = modify->find_compute(arg[1]); - if (icompute < 0) - error->all(FLERR,"Could not find fix_modify temperature ID"); - temperature = modify->compute[icompute]; + temperature = modify->get_compute_by_id(id_temp); + if (!temperature) error->all(FLERR,"Could not find fix_modify temperature ID {}", id_temp); if (temperature->tempflag == 0) - error->all(FLERR, - "Fix_modify temperature ID does not compute temperature"); + error->all(FLERR, "Fix_modify temperature ID {} does not compute temperature", id_temp); if (temperature->igroup != 0 && comm->me == 0) error->warning(FLERR,"Temperature for fix modify is not for group all"); // reset id_temp of pressure to new temperature ID if (pstat_flag) { - icompute = modify->find_compute(id_press); - if (icompute < 0) - error->all(FLERR,"Pressure ID for fix modify does not exist"); - modify->compute[icompute]->reset_extra_compute_fix(id_temp); + pressure = modify->get_compute_by_id(id_press); + if (!pressure) error->all(FLERR,"Pressure ID {} for fix modify does not exist", id_press); + pressure->reset_extra_compute_fix(id_temp); } return 2; @@ -1458,15 +1430,14 @@ int FixTGNHDrude::modify_param(int narg, char **arg) modify->delete_compute(id_press); pcomputeflag = 0; } - delete [] id_press; + delete[] id_press; id_press = utils::strdup(arg[1]); - int icompute = modify->find_compute(arg[1]); - if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID"); - pressure = modify->compute[icompute]; + pressure = modify->get_compute_by_id(id_press); + if (!pressure) error->all(FLERR,"Could not find fix_modify pressure ID {}", id_press); if (pressure->pressflag == 0) - error->all(FLERR,"Fix_modify pressure ID does not compute pressure"); + error->all(FLERR,"Fix_modify pressure ID {} does not compute pressure", id_press); return 2; } diff --git a/src/DRUDE/fix_tgnh_drude.h b/src/DRUDE/fix_tgnh_drude.h index adfa69671a..b2724809b4 100644 --- a/src/DRUDE/fix_tgnh_drude.h +++ b/src/DRUDE/fix_tgnh_drude.h @@ -63,8 +63,7 @@ class FixTGNHDrude : public Fix { double omega_mass[6]; double p_current[6]; int kspace_flag; // 1 if KSpace invoked, 0 if not - int nrigid; // number of rigid fixes - int *rfix; // indices of rigid fixes + std::vector rfix; // indices of rigid fixes class Irregular *irregular; // for migrating atoms after box flips int nlevels_respa; diff --git a/src/EXTRA-FIX/fix_npt_cauchy.cpp b/src/EXTRA-FIX/fix_npt_cauchy.cpp index feb5a95c6f..f3dfd1af36 100644 --- a/src/EXTRA-FIX/fix_npt_cauchy.cpp +++ b/src/EXTRA-FIX/fix_npt_cauchy.cpp @@ -54,14 +54,12 @@ enum{ISO,ANISO,TRICLINIC}; ---------------------------------------------------------------------- */ FixNPTCauchy::FixNPTCauchy(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), - rfix(nullptr), id_dilate(nullptr), irregular(nullptr), - id_temp(nullptr), id_press(nullptr), - eta(nullptr), eta_dot(nullptr), eta_dotdot(nullptr), - eta_mass(nullptr), etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr), - etap_mass(nullptr), id_store(nullptr), init_store(nullptr) + Fix(lmp, narg, arg), id_dilate(nullptr), irregular(nullptr), id_temp(nullptr), + id_press(nullptr), eta(nullptr), eta_dot(nullptr), eta_dotdot(nullptr), eta_mass(nullptr), + etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr), etap_mass(nullptr), id_store(nullptr), + init_store(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal fix npt/cauchy command"); + if (narg < 4) error->all(FLERR, "Illegal fix npt/cauchy command"); dynamic_group_allow = 1; ecouple_flag = 1; @@ -571,9 +569,6 @@ FixNPTCauchy::FixNPTCauchy(LAMMPS *lmp, int narg, char **arg) : if (deviatoric_flag) size_vector += 1; } - nrigid = 0; - rfix = nullptr; - if (pre_exchange_flag) irregular = new Irregular(lmp); else irregular = nullptr; @@ -619,8 +614,6 @@ FixNPTCauchy::~FixNPTCauchy() if (copymode) return; delete[] id_dilate; - delete[] rfix; - delete[] id_store; delete irregular; @@ -690,19 +683,16 @@ void FixNPTCauchy::init() // set temperature and pressure ptrs - int icompute = modify->find_compute(id_temp); - if (icompute < 0) - error->all(FLERR,"Temperature ID for fix npt/cauchy does not exist"); - temperature = modify->compute[icompute]; + temperature = modify->get_compute_by_id(id_temp); + if (!temperature) + error->all(FLERR,"Temperature ID {} for fix npt/cauchy does not exist", id_temp); if (temperature->tempbias) which = BIAS; else which = NOBIAS; if (pstat_flag) { - icompute = modify->find_compute(id_press); - if (icompute < 0) - error->all(FLERR,"Pressure ID for fix npt/cauchy does not exist"); - pressure = modify->compute[icompute]; + pressure = modify->get_compute_by_id(id_press); + if (!pressure) error->all(FLERR,"Pressure ID {} for fix npt/cauchy does not exist", id_press); } // set timesteps and frequencies @@ -759,20 +749,10 @@ void FixNPTCauchy::init() } // detect if any rigid fixes exist so rigid bodies move when box is remapped - // rfix[] = indices to each fix rigid - delete[] rfix; - nrigid = 0; - rfix = nullptr; - - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) nrigid++; - if (nrigid) { - rfix = new int[nrigid]; - nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (modify->fix[i]->rigid_flag) rfix[nrigid++] = i; - } + rfix.clear(); + for (auto &ifix : modify->get_fix_list()) + if (ifix->rigid_flag) rfix.push_back(ifix); } /* ---------------------------------------------------------------------- @@ -1121,9 +1101,7 @@ void FixNPTCauchy::remap() domain->x2lamda(x[i],x[i]); } - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(0); + for (auto &ifix : rfix) ifix->deform(0); // reset global and local box to new size/shape @@ -1268,9 +1246,7 @@ void FixNPTCauchy::remap() domain->lamda2x(x[i],x[i]); } - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(1); + for (auto &ifix : rfix) ifix->deform(1); } /* ---------------------------------------------------------------------- @@ -1432,23 +1408,20 @@ int FixNPTCauchy::modify_param(int narg, char **arg) delete[] id_temp; id_temp = utils::strdup(arg[1]); - int icompute = modify->find_compute(arg[1]); - if (icompute < 0) - error->all(FLERR,"Could not find fix_modify temperature ID"); - temperature = modify->compute[icompute]; + temperature = modify->get_compute_by_id(id_temp); + if (!temperature) error->all(FLERR,"Could not find fix_modify temperature ID {}", id_temp); if (temperature->tempflag == 0) - error->all(FLERR,"Fix_modify temperature ID does not compute temperature"); + error->all(FLERR,"Fix_modify temperature ID {} does not compute temperature", id_temp); if (temperature->igroup != 0 && comm->me == 0) error->warning(FLERR,"Temperature for fix modify is not for group all"); // reset id_temp of pressure to new temperature ID if (pstat_flag) { - icompute = modify->find_compute(id_press); - if (icompute < 0) - error->all(FLERR,"Pressure ID for fix modify does not exist"); - modify->compute[icompute]->reset_extra_compute_fix(id_temp); + pressure = modify->get_compute_by_id(id_press); + if (!pressure) error->all(FLERR,"Pressure ID {} for fix modify does not exist", id_press); + pressure->reset_extra_compute_fix(id_temp); } return 2; @@ -1463,12 +1436,11 @@ int FixNPTCauchy::modify_param(int narg, char **arg) delete[] id_press; id_press = utils::strdup(arg[1]); - int icompute = modify->find_compute(arg[1]); - if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID"); - pressure = modify->compute[icompute]; + pressure = modify->get_compute_by_id(id_press); + if (!pressure) error->all(FLERR,"Could not find fix_modify pressure ID {}", id_press); if (pressure->pressflag == 0) - error->all(FLERR,"Fix_modify pressure ID does not compute pressure"); + error->all(FLERR,"Fix_modify pressure ID {} does not compute pressure", id_press); return 2; } diff --git a/src/EXTRA-FIX/fix_npt_cauchy.h b/src/EXTRA-FIX/fix_npt_cauchy.h index e7e6630208..43a944acb4 100644 --- a/src/EXTRA-FIX/fix_npt_cauchy.h +++ b/src/EXTRA-FIX/fix_npt_cauchy.h @@ -73,9 +73,8 @@ class FixNPTCauchy : public Fix { double drag, tdrag_factor; // drag factor on particle thermostat double pdrag_factor; // drag factor on barostat int kspace_flag; // 1 if KSpace invoked, 0 if not - int nrigid; // number of rigid fixes int dilate_group_bit; // mask for dilation group - int *rfix; // indices of rigid fixes + std::vector rfix; // indices of rigid fixes char *id_dilate; // group name to dilate class Irregular *irregular; // for migrating atoms after box flips diff --git a/src/QTB/fix_qbmsst.cpp b/src/QTB/fix_qbmsst.cpp index b5fb5ca77c..2450561363 100644 --- a/src/QTB/fix_qbmsst.cpp +++ b/src/QTB/fix_qbmsst.cpp @@ -228,7 +228,6 @@ FixQBMSST::FixQBMSST(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) pressure = nullptr; pe = nullptr; old_velocity = nullptr; - rfix = nullptr; gfactor = nullptr; random = nullptr; omega_H = nullptr; @@ -263,17 +262,16 @@ FixQBMSST::FixQBMSST(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) FixQBMSST::~FixQBMSST() { - delete [] rfix; - delete [] gfactor; + delete[] gfactor; delete random; // delete temperature and pressure if fix created them if (tflag) modify->delete_compute(id_temp); if (pflag) modify->delete_compute(id_press); if (peflag) modify->delete_compute(id_pe); - delete [] id_temp; - delete [] id_press; - delete [] id_pe; + delete[] id_temp; + delete[] id_press; + delete[] id_pe; memory->destroy(old_velocity); memory->destroy(fran); @@ -385,18 +383,10 @@ void FixQBMSST::init() else kspace_flag = 0; // detect if any fix rigid exist so rigid bodies move when box is dilated - // rfix[] = indices to each fix rigid - nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (utils::strmatch(modify->fix[i]->style,"^rigid") || - (strcmp(modify->fix[i]->style,"poems") == 0)) nrigid++; - if (nrigid > 0) { - rfix = new int[nrigid]; - nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (utils::strmatch(modify->fix[i]->style,"^rigid") || - (strcmp(modify->fix[i]->style,"poems") == 0)) rfix[nrigid++] = i; - } + + rfix.clear(); + for (auto &ifix : modify->get_fix_list()) + if (ifix->rigid_flag) rfix.push_back(ifix); } /* ---------------------------------------------------------------------- @@ -787,9 +777,7 @@ void FixQBMSST::remap(int flag) domain->x2lamda(n); - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(0); + for (auto &ifix : rfix) ifix->deform(0); // reset global and local box to new size/shape @@ -810,9 +798,7 @@ void FixQBMSST::remap(int flag) domain->lamda2x(n); - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(1); + for (auto &ifix : rfix) ifix->deform(1); for (i = 0; i < n; i++) { v[i][direction] = v[i][direction] * @@ -868,7 +854,7 @@ int FixQBMSST::modify_param(int narg, char **arg) modify->delete_compute(id_temp); tflag = 0; } - delete [] id_temp; + delete[] id_temp; id_temp = utils::strdup(arg[1]); int icompute = modify->find_compute(id_temp); @@ -888,7 +874,7 @@ int FixQBMSST::modify_param(int narg, char **arg) modify->delete_compute(id_press); pflag = 0; } - delete [] id_press; + delete[] id_press; id_press = utils::strdup(arg[1]); int icompute = modify->find_compute(id_press); diff --git a/src/QTB/fix_qbmsst.h b/src/QTB/fix_qbmsst.h index ecfa5abf8e..cccb4e6a17 100644 --- a/src/QTB/fix_qbmsst.h +++ b/src/QTB/fix_qbmsst.h @@ -78,8 +78,7 @@ class FixQBMSST : public Fix { double omega[3]; // Time derivative of the volume. double total_mass; // Mass of the computational cell int kspace_flag; // 1 if KSpace invoked, 0 if not - int nrigid; // number of rigid fixes - int *rfix; // indices of rigid fixes + std::vector rfix; // indices of rigid fixes double p_current[3]; // pressure double velocity_sum; // Sum of the velocities squared. double lagrangian_position; // Lagrangian location of computational cell diff --git a/src/SHOCK/fix_msst.cpp b/src/SHOCK/fix_msst.cpp index a4c9db3fd7..55842250ec 100644 --- a/src/SHOCK/fix_msst.cpp +++ b/src/SHOCK/fix_msst.cpp @@ -41,11 +41,10 @@ using namespace FixConst; /* ---------------------------------------------------------------------- */ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), old_velocity(nullptr), rfix(nullptr), - id_temp(nullptr), id_press(nullptr), id_pe(nullptr), temperature(nullptr), - pressure(nullptr), pe(nullptr) + Fix(lmp, narg, arg), old_velocity(nullptr), id_temp(nullptr), id_press(nullptr), id_pe(nullptr), + temperature(nullptr), pressure(nullptr), pe(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal fix msst command"); + if (narg < 4) error->all(FLERR, "Illegal fix msst command"); restart_global = 1; time_integrate = 1; @@ -80,95 +79,103 @@ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) : dftb = 0; beta = 0.0; - if (strcmp(arg[3],"x") == 0) { + if (strcmp(arg[3], "x") == 0) { direction = 0; box_change |= BOX_CHANGE_X; - } else if (strcmp(arg[3],"y") == 0) { + } else if (strcmp(arg[3], "y") == 0) { direction = 1; box_change |= BOX_CHANGE_Y; - } else if (strcmp(arg[3],"z") == 0) { + } else if (strcmp(arg[3], "z") == 0) { direction = 2; box_change |= BOX_CHANGE_Z; - } else error->all(FLERR,"Illegal fix msst command"); + } else + error->all(FLERR, "Illegal fix msst command"); - velocity = utils::numeric(FLERR,arg[4],false,lmp); - if (velocity < 0) error->all(FLERR,"Illegal fix msst command"); + velocity = utils::numeric(FLERR, arg[4], false, lmp); + if (velocity < 0) error->all(FLERR, "Illegal fix msst command"); // optional args int iarg = 5; while (iarg < narg) { - if (strcmp(arg[iarg],"q") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command"); - qmass = utils::numeric(FLERR,arg[iarg+1],false,lmp); + if (strcmp(arg[iarg], "q") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix msst command"); + qmass = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"mu") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command"); - mu = utils::numeric(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "mu") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix msst command"); + mu = utils::numeric(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"p0") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command"); - p0 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "p0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix msst command"); + p0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); p0_set = 1; iarg += 2; - } else if (strcmp(arg[iarg],"v0") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command"); - v0 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "v0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix msst command"); + v0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); v0_set = 1; iarg += 2; - } else if (strcmp(arg[iarg],"e0") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command"); - e0 = utils::numeric(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "e0") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix msst command"); + e0 = utils::numeric(FLERR, arg[iarg + 1], false, lmp); e0_set = 1; iarg += 2; - } else if (strcmp(arg[iarg],"tscale") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command"); - tscale = utils::numeric(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "tscale") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix msst command"); + tscale = utils::numeric(FLERR, arg[iarg + 1], false, lmp); if (tscale < 0.0 || tscale > 1.0) - error->all(FLERR,"Fix msst tscale must satisfy 0 <= tscale < 1"); + error->all(FLERR, "Fix msst tscale must satisfy 0 <= tscale < 1"); iarg += 2; - } else if (strcmp(arg[iarg],"dftb") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command"); - dftb = utils::logical(FLERR,arg[iarg+1],false,lmp); + } else if (strcmp(arg[iarg], "dftb") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix msst command"); + dftb = utils::logical(FLERR, arg[iarg + 1], false, lmp); iarg += 2; - } else if (strcmp(arg[iarg],"beta") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix msst command"); - beta = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (beta < 0.0 || beta > 1.0) - error->all(FLERR,"Illegal fix msst command"); + } else if (strcmp(arg[iarg], "beta") == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Illegal fix msst command"); + beta = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (beta < 0.0 || beta > 1.0) error->all(FLERR, "Illegal fix msst command"); iarg += 2; - } else error->all(FLERR,"Illegal fix msst command"); + } else + error->all(FLERR, "Illegal fix msst command"); } // output MSST info if (comm->me == 0) { std::string mesg = "MSST parameters:\n"; - if (direction == 0) mesg += " Shock in x direction\n"; - else if (direction == 1) mesg += " Shock in y direction\n"; - else if (direction == 2) mesg += " Shock in z direction\n"; + if (direction == 0) + mesg += " Shock in x direction\n"; + else if (direction == 1) + mesg += " Shock in y direction\n"; + else if (direction == 2) + mesg += " Shock in z direction\n"; mesg += fmt::format(" Cell mass-like parameter qmass " - "(units of mass^2/length^4) = {:.8g}\n", qmass); + "(units of mass^2/length^4) = {:.8g}\n", + qmass); mesg += fmt::format(" Shock velocity = {:.8g}\n", velocity); mesg += fmt::format(" Artificial viscosity (units of mass/length/time) = {:.8g}\n", mu); if (p0_set) mesg += fmt::format(" Initial pressure specified to be {:.8g}\n", p0); - else mesg += " Initial pressure calculated on first step\n"; + else + mesg += " Initial pressure calculated on first step\n"; if (v0_set) mesg += fmt::format(" Initial volume specified to be {:.8g}\n", v0); - else mesg += " Initial volume calculated on first step\n"; + else + mesg += " Initial volume calculated on first step\n"; if (e0_set) mesg += fmt::format(" Initial energy specified to be {:.8g}\n", e0); - else mesg += " Initial energy calculated on first step\n"; - utils::logmesg(lmp,mesg); + else + mesg += " Initial energy calculated on first step\n"; + utils::logmesg(lmp, mesg); } // check for periodicity in controlled dimensions - if (domain->nonperiodic) error->all(FLERR,"Fix msst requires a periodic box"); + if (domain->nonperiodic) error->all(FLERR, "Fix msst requires a periodic box"); // create a new temperature compute // id = fix-ID + "MSST_temp" @@ -200,8 +207,6 @@ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) : // initialize the time derivative of the volume omega[0] = omega[1] = omega[2] = 0.0; - nrigid = 0; - rfix = nullptr; maxold = -1; old_velocity = nullptr; @@ -211,17 +216,15 @@ FixMSST::FixMSST(LAMMPS *lmp, int narg, char **arg) : FixMSST::~FixMSST() { - delete [] rfix; - // delete temperature and pressure if fix created them if (tflag) modify->delete_compute(id_temp); if (pflag) modify->delete_compute(id_press); if (peflag) modify->delete_compute(id_pe); - delete [] id_temp; - delete [] id_press; - delete [] id_pe; + delete[] id_temp; + delete[] id_press; + delete[] id_pe; memory->destroy(old_velocity); } @@ -240,22 +243,20 @@ int FixMSST::setmask() void FixMSST::init() { - if (atom->mass == nullptr) - error->all(FLERR,"Cannot use fix msst without per-type mass defined"); + if (atom->mass == nullptr) error->all(FLERR, "Cannot use fix msst without per-type mass defined"); // set compute ptrs int itemp = modify->find_compute(id_temp); int ipress = modify->find_compute(id_press); int ipe = modify->find_compute(id_pe); - if (itemp < 0 || ipress < 0|| ipe < 0) - error->all(FLERR,"Could not find fix msst compute ID"); + if (itemp < 0 || ipress < 0 || ipe < 0) error->all(FLERR, "Could not find fix msst compute ID"); if (modify->compute[itemp]->tempflag == 0) - error->all(FLERR,"Fix msst compute ID does not compute temperature"); + error->all(FLERR, "Fix msst compute ID does not compute temperature"); if (modify->compute[ipress]->pressflag == 0) - error->all(FLERR,"Fix msst compute ID does not compute pressure"); + error->all(FLERR, "Fix msst compute ID does not compute pressure"); if (modify->compute[ipe]->peflag == 0) - error->all(FLERR,"Fix msst compute ID does not compute potential energy"); + error->all(FLERR, "Fix msst compute ID does not compute potential energy"); temperature = modify->compute[itemp]; pressure = modify->compute[ipress]; @@ -271,37 +272,27 @@ void FixMSST::init() double mass = 0.0; for (int i = 0; i < atom->nlocal; i++) mass += atom->mass[atom->type[i]]; - MPI_Allreduce(&mass,&total_mass,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&mass, &total_mass, 1, MPI_DOUBLE, MPI_SUM, world); - if (force->kspace) kspace_flag = 1; - else kspace_flag = 0; + if (force->kspace) + kspace_flag = 1; + else + kspace_flag = 0; // detect if any fix rigid exist so rigid bodies move when box is dilated - // rfix[] = indices to each fix rigid - delete [] rfix; - nrigid = 0; - rfix = nullptr; - - for (int i = 0; i < modify->nfix; i++) - if (utils::strmatch(modify->fix[i]->style,"^rigid") || - utils::strmatch(modify->fix[i]->style,"^poems$")) nrigid++; - if (nrigid) { - rfix = new int[nrigid]; - nrigid = 0; - for (int i = 0; i < modify->nfix; i++) - if (utils::strmatch(modify->fix[i]->style,"^rigid") || - utils::strmatch(modify->fix[i]->style,"^poems$")) rfix[nrigid++] = i; - } + rfix.clear(); + for (auto &ifix : modify->get_fix_list()) + if (ifix->rigid_flag) rfix.push_back(ifix); // find fix external being used to drive LAMMPS from DFTB+ if (dftb) { for (int i = 0; i < modify->nfix; i++) - if (utils::strmatch(modify->fix[i]->style,"^external$")) + if (utils::strmatch(modify->fix[i]->style, "^external$")) fix_external = dynamic_cast(modify->fix[i]); if (fix_external == nullptr) - error->all(FLERR,"Fix msst dftb cannot be used w/out fix external"); + error->all(FLERR, "Fix msst dftb cannot be used w/out fix external"); } } @@ -321,29 +312,26 @@ void FixMSST::setup(int /*vflag*/) if (v0_set == 0) { v0 = compute_vol(); v0_set = 1; - if (comm->me == 0) - utils::logmesg(lmp,"Fix MSST v0 = {:.8g}\n", v0); + if (comm->me == 0) utils::logmesg(lmp, "Fix MSST v0 = {:.8g}\n", v0); } if (p0_set == 0) { p0 = p_current[direction]; p0_set = 1; - if (comm->me == 0) - utils::logmesg(lmp,"Fix MSST p0 = {:.8g}\n", p0); + if (comm->me == 0) utils::logmesg(lmp, "Fix MSST p0 = {:.8g}\n", p0); } if (e0_set == 0) { e0 = compute_etotal(); e0_set = 1; - if (comm->me == 0) - utils::logmesg(lmp,"Fix MSST e0 = {:.8g}\n", e0); + if (comm->me == 0) utils::logmesg(lmp, "Fix MSST e0 = {:.8g}\n", e0); } temperature->compute_vector(); double *ke_tensor = temperature->vector; - double ke_temp = ke_tensor[0]+ke_tensor[1]+ke_tensor[2]; + double ke_temp = ke_tensor[0] + ke_tensor[1] + ke_tensor[2]; if (ke_temp > 0.0 && tscale > 0.0) { // transfer energy from atom velocities to cell volume motion @@ -351,30 +339,30 @@ void FixMSST::setup(int /*vflag*/) double **v = atom->v; int *mask = atom->mask; - double sqrt_initial_temperature_scaling = sqrt(1.0-tscale); + double sqrt_initial_temperature_scaling = sqrt(1.0 - tscale); - double fac1 = tscale*total_mass/qmass*ke_temp/force->mvv2e; + double fac1 = tscale * total_mass / qmass * ke_temp / force->mvv2e; - omega[direction]=-1*sqrt(fac1); - double fac2 = omega[direction]/v0; + omega[direction] = -1 * sqrt(fac1); + double fac2 = omega[direction] / v0; - if ( comm->me == 0 && tscale != 1.0) - utils::logmesg(lmp,"Fix MSST initial strain rate of {:.8g} " + if (comm->me == 0 && tscale != 1.0) + utils::logmesg(lmp, + "Fix MSST initial strain rate of {:.8g} " "established by reducing temperature by factor " - "of {:.8g}\n",fac2,tscale); + "of {:.8g}\n", + fac2, tscale); for (int i = 0; i < atom->nlocal; i++) { if (mask[i] & groupbit) { - for (int k = 0; k < 3; k++) { - v[i][k]*=sqrt_initial_temperature_scaling; - } + for (int k = 0; k < 3; k++) { v[i][k] *= sqrt_initial_temperature_scaling; } } } } // trigger virial computation on next timestep - pe->addstep(update->ntimestep+1); - pressure->addstep(update->ntimestep+1); + pe->addstep(update->ntimestep + 1); + pressure->addstep(update->ntimestep + 1); } /* ---------------------------------------------------------------------- @@ -383,8 +371,8 @@ void FixMSST::setup(int /*vflag*/) void FixMSST::initial_integrate(int /*vflag*/) { - int i,k; - double p_msst; // MSST driving pressure + int i, k; + double p_msst; // MSST driving pressure double vol; int nlocal = atom->nlocal; @@ -401,7 +389,7 @@ void FixMSST::initial_integrate(int /*vflag*/) if (nlocal > maxold) { memory->destroy(old_velocity); maxold = atom->nmax; - memory->create(old_velocity,maxold,3,"msst:old_velocity"); + memory->create(old_velocity, maxold, 3, "msst:old_velocity"); } // for DFTB, extract TS_dftb from fix external @@ -409,14 +397,14 @@ void FixMSST::initial_integrate(int /*vflag*/) if (dftb) { const double TS_dftb = fix_external->compute_vector(0); - const double TS = force->ftm2v*TS_dftb; + const double TS = force->ftm2v * TS_dftb; // update S_elec terms and compute TS_dot via finite differences S_elec_2 = S_elec_1; S_elec_1 = S_elec; const double Temp = temperature->compute_scalar(); - S_elec = TS/Temp; - TS_dot = Temp*(3.0*S_elec-4.0*S_elec_1+S_elec_2)/(2.0*update->dt); - TS_int += (update->dt*TS_dot); + S_elec = TS / Temp; + TS_dot = Temp * (3.0 * S_elec - 4.0 * S_elec_1 + S_elec_2) / (2.0 * update->dt); + TS_int += (update->dt * TS_dot); if (update->ntimestep == 1) T0S0 = TS; } @@ -434,11 +422,9 @@ void FixMSST::initial_integrate(int /*vflag*/) // propagate the time derivative of // the volume 1/2 step at fixed vol, r, rdot - p_msst = nktv2p * mvv2e * velocity * velocity * total_mass * - ( v0 - vol)/( v0 * v0); - double A = total_mass * ( p_current[sd] - p0 - p_msst ) / - (qmass * nktv2p * mvv2e); - double B = total_mass * mu / ( qmass * vol ); + p_msst = nktv2p * mvv2e * velocity * velocity * total_mass * (v0 - vol) / (v0 * v0); + double A = total_mass * (p_current[sd] - p0 - p_msst) / (qmass * nktv2p * mvv2e); + double B = total_mass * mu / (qmass * vol); // prevent blow-up of the volume @@ -447,11 +433,10 @@ void FixMSST::initial_integrate(int /*vflag*/) // use Taylor expansion to avoid singularity at B = 0 if (B * dthalf > 1.0e-06) { - omega[sd] = ( omega[sd] + A * ( exp(B * dthalf) - 1.0 ) / B ) - * exp(-B * dthalf); + omega[sd] = (omega[sd] + A * (exp(B * dthalf) - 1.0) / B) * exp(-B * dthalf); } else { omega[sd] = omega[sd] + (A - B * omega[sd]) * dthalf + - 0.5 * (B * B * omega[sd] - A * B ) * dthalf * dthalf; + 0.5 * (B * B * omega[sd] - A * B) * dthalf * dthalf; } // propagate velocity sum 1/2 step by @@ -464,20 +449,19 @@ void FixMSST::initial_integrate(int /*vflag*/) if (mask[i] & groupbit) { for (k = 0; k < 3; k++) { const double C = f[i][k] * force->ftm2v / mass[type[i]]; - const double TS_term = TS_dot/(mass[type[i]]*velocity_sum); - const double 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 ); + const double TS_term = TS_dot / (mass[type[i]] * velocity_sum); + const double 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); D += escale_term - TS_term; old_velocity[i][k] = v[i][k]; if (k == direction) D -= 2.0 * omega[sd] / vol; if (fabs(dthalf * D) > 1.0e-06) { const double expd = exp(D * dthalf); - v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D; + v[i][k] = expd * (C + D * v[i][k] - C / expd) / D; } else { - v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf + - 0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf; + v[i][k] = v[i][k] + (C + D * v[i][k]) * dthalf + + 0.5 * (D * D * v[i][k] + C * D) * dthalf * dthalf; } } } @@ -487,18 +471,15 @@ void FixMSST::initial_integrate(int /*vflag*/) if (mask[i] & groupbit) { for (k = 0; k < 3; k++) { const double C = f[i][k] * force->ftm2v / mass[type[i]]; - double D = mu * omega[sd] * omega[sd] / - (velocity_sum * mass[type[i]] * vol ); + double D = mu * omega[sd] * omega[sd] / (velocity_sum * mass[type[i]] * vol); old_velocity[i][k] = v[i][k]; - if (k == direction) { - D -= 2.0 * omega[sd] / vol; - } + if (k == direction) { D -= 2.0 * omega[sd] / vol; } if (fabs(dthalf * D) > 1.0e-06) { const double expd = exp(D * dthalf); - v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D; + v[i][k] = expd * (C + D * v[i][k] - C / expd) / D; } else { - v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf + - 0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf; + v[i][k] = v[i][k] + (C + D * v[i][k]) * dthalf + + 0.5 * (D * D * v[i][k] + C * D) * dthalf * dthalf; } } } @@ -524,19 +505,18 @@ void FixMSST::initial_integrate(int /*vflag*/) if (mask[i] & groupbit) { for (k = 0; k < 3; k++) { const double C = f[i][k] * force->ftm2v / mass[type[i]]; - const double TS_term = TS_dot/(mass[type[i]]*velocity_sum); - const double 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 ); + const double TS_term = TS_dot / (mass[type[i]] * velocity_sum); + const double 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); D += escale_term - TS_term; if (k == direction) D -= 2.0 * omega[sd] / vol; if (fabs(dthalf * D) > 1.0e-06) { const double expd = exp(D * dthalf); - v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D; + v[i][k] = expd * (C + D * v[i][k] - C / expd) / D; } else { - v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf + - 0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf; + v[i][k] = v[i][k] + (C + D * v[i][k]) * dthalf + + 0.5 * (D * D * v[i][k] + C * D) * dthalf * dthalf; } } } @@ -546,17 +526,14 @@ void FixMSST::initial_integrate(int /*vflag*/) if (mask[i] & groupbit) { for (k = 0; k < 3; k++) { const double C = f[i][k] * force->ftm2v / mass[type[i]]; - double D = mu * omega[sd] * omega[sd] / - (velocity_sum * mass[type[i]] * vol ); - if (k == direction) { - D -= 2.0 * omega[sd] / vol; - } + double D = mu * omega[sd] * omega[sd] / (velocity_sum * mass[type[i]] * vol); + if (k == direction) { D -= 2.0 * omega[sd] / vol; } if (fabs(dthalf * D) > 1.0e-06) { const double expd = exp(D * dthalf); - v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D; + v[i][k] = expd * (C + D * v[i][k] - C / expd) / D; } else { - v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf + - 0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf; + v[i][k] = v[i][k] + (C + D * v[i][k]) * dthalf + + 0.5 * (D * D * v[i][k] + C * D) * dthalf * dthalf; } } } @@ -569,7 +546,7 @@ void FixMSST::initial_integrate(int /*vflag*/) // rescale positions and change box size - dilation[sd] = vol1/vol; + dilation[sd] = vol1 / vol; remap(0); // propagate particle positions 1 time step @@ -588,7 +565,7 @@ void FixMSST::initial_integrate(int /*vflag*/) // rescale positions and change box size - dilation[sd] = vol2/vol1; + dilation[sd] = vol2 / vol1; remap(0); if (kspace_flag) force->kspace->setup(); @@ -601,7 +578,7 @@ void FixMSST::initial_integrate(int /*vflag*/) void FixMSST::final_integrate() { int i; - double p_msst; // MSST driving pressure + double p_msst; // MSST driving pressure // v update only for atoms in MSST group @@ -624,14 +601,14 @@ void FixMSST::final_integrate() if (dftb) { const double TS_dftb = fix_external->compute_vector(0); - const double TS = force->ftm2v*TS_dftb; + const double TS = force->ftm2v * TS_dftb; S_elec_2 = S_elec_1; S_elec_1 = S_elec; const double Temp = temperature->compute_scalar(); // update S_elec terms and compute TS_dot via finite differences - S_elec = TS/Temp; - TS_dot = Temp*(3.0*S_elec-4.0*S_elec_1+S_elec_2)/(2.0*update->dt); - TS_int += (update->dt*TS_dot); + S_elec = TS / Temp; + TS_dot = Temp * (3.0 * S_elec - 4.0 * S_elec_1 + S_elec_2) / (2.0 * update->dt); + TS_int += (update->dt * TS_dot); if (update->ntimestep == 1) T0S0 = TS; } @@ -642,19 +619,18 @@ void FixMSST::final_integrate() if (mask[i] & groupbit) { for (int k = 0; k < 3; k++) { const double C = f[i][k] * force->ftm2v / mass[type[i]]; - const double TS_term = TS_dot/(mass[type[i]]*velocity_sum); - const double 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 ); + const double TS_term = TS_dot / (mass[type[i]] * velocity_sum); + const double 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); D += escale_term - TS_term; if (k == direction) D -= 2.0 * omega[sd] / vol; if (fabs(dthalf * D) > 1.0e-06) { const double expd = exp(D * dthalf); - v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D; + v[i][k] = expd * (C + D * v[i][k] - C / expd) / D; } else { - v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf + - 0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf; + v[i][k] = v[i][k] + (C + D * v[i][k]) * dthalf + + 0.5 * (D * D * v[i][k] + C * D) * dthalf * dthalf; } } } @@ -664,17 +640,14 @@ void FixMSST::final_integrate() if (mask[i] & groupbit) { for (int k = 0; k < 3; k++) { const double C = f[i][k] * force->ftm2v / mass[type[i]]; - double D = mu * omega[sd] * omega[sd] / - (velocity_sum * mass[type[i]] * vol ); - if (k == direction) { - D -= 2.0 * omega[sd] / vol; - } + double D = mu * omega[sd] * omega[sd] / (velocity_sum * mass[type[i]] * vol); + if (k == direction) { D -= 2.0 * omega[sd] / vol; } if (fabs(dthalf * D) > 1.0e-06) { const double expd = exp(D * dthalf); - v[i][k] = expd * ( C + D * v[i][k] - C / expd ) / D; + v[i][k] = expd * (C + D * v[i][k] - C / expd) / D; } else { - v[i][k] = v[i][k] + ( C + D * v[i][k] ) * dthalf + - 0.5 * (D * D * v[i][k] + C * D ) * dthalf * dthalf; + v[i][k] = v[i][k] + (C + D * v[i][k]) * dthalf + + 0.5 * (D * D * v[i][k] + C * D) * dthalf * dthalf; } } } @@ -692,11 +665,9 @@ void FixMSST::final_integrate() // propagate the time derivative of the volume 1/2 step at fixed V, r, rdot - p_msst = nktv2p * mvv2e * velocity * velocity * total_mass * - ( v0 - vol )/( v0 * v0 ); - double A = total_mass * ( p_current[sd] - p0 - p_msst ) / - ( qmass * nktv2p * mvv2e ); - const double B = total_mass * mu / ( qmass * vol ); + p_msst = nktv2p * mvv2e * velocity * velocity * total_mass * (v0 - vol) / (v0 * v0); + double A = total_mass * (p_current[sd] - p0 - p_msst) / (qmass * nktv2p * mvv2e); + const double B = total_mass * mu / (qmass * vol); // prevent blow-up of the volume @@ -705,21 +676,20 @@ void FixMSST::final_integrate() // use taylor expansion to avoid singularity at B == 0. if (B * dthalf > 1.0e-06) { - omega[sd] = ( omega[sd] + A * - ( exp(B * dthalf) - 1.0 ) / B ) * exp(-B * dthalf); + omega[sd] = (omega[sd] + A * (exp(B * dthalf) - 1.0) / B) * exp(-B * dthalf); } else { omega[sd] = omega[sd] + (A - B * omega[sd]) * dthalf + - 0.5 * (B * B * omega[sd] - A * B ) * dthalf * dthalf; + 0.5 * (B * B * omega[sd] - A * B) * dthalf * dthalf; } // calculate Lagrangian position of computational cell - lagrangian_position -= velocity*vol/v0*update->dt; + lagrangian_position -= velocity * vol / v0 * update->dt; // trigger energy and virial computation on next timestep - pe->addstep(update->ntimestep+1); - pressure->addstep(update->ntimestep+1); + pe->addstep(update->ntimestep + 1); + pressure->addstep(update->ntimestep + 1); } /* ---------------------------------------------------------------------- */ @@ -741,20 +711,20 @@ void FixMSST::couple() void FixMSST::remap(int flag) { - int i,n; - double oldlo,oldhi,ctr; + int i, n; + double oldlo, oldhi, ctr; double **v = atom->v; - if (flag) n = atom->nlocal + atom->nghost; - else n = atom->nlocal; + if (flag) + n = atom->nlocal + atom->nghost; + else + n = atom->nlocal; // convert pertinent atoms and rigid bodies to lamda coords domain->x2lamda(n); - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(0); + for (auto &ifix : rfix) ifix->deform(0); // reset global and local box to new size/shape @@ -763,8 +733,8 @@ void FixMSST::remap(int flag) oldlo = domain->boxlo[i]; oldhi = domain->boxhi[i]; ctr = 0.5 * (oldlo + oldhi); - domain->boxlo[i] = (oldlo-ctr)*dilation[i] + ctr; - domain->boxhi[i] = (oldhi-ctr)*dilation[i] + ctr; + domain->boxlo[i] = (oldlo - ctr) * dilation[i] + ctr; + domain->boxhi[i] = (oldhi - ctr) * dilation[i] + ctr; } } @@ -775,14 +745,9 @@ void FixMSST::remap(int flag) domain->lamda2x(n); - if (nrigid) - for (i = 0; i < nrigid; i++) - modify->fix[rfix[i]]->deform(1); + for (auto &ifix : rfix) ifix->deform(1); - for (i = 0; i < n; i++) { - v[i][direction] = v[i][direction] * - dilation[direction]; - } + for (i = 0; i < n; i++) v[i][direction] = v[i][direction] * dilation[direction]; } /* ---------------------------------------------------------------------- @@ -800,8 +765,8 @@ void FixMSST::write_restart(FILE *fp) list[n++] = TS_int; if (comm->me == 0) { int size = n * sizeof(double); - fwrite(&size,sizeof(int),1,fp); - fwrite(&list,sizeof(double),n,fp); + fwrite(&size, sizeof(int), 1, fp); + fwrite(&list, sizeof(double), n, fp); } } @@ -818,7 +783,7 @@ void FixMSST::restart(char *buf) v0 = list[n++]; p0 = list[n++]; TS_int = list[n++]; - tscale = 0.0; // set tscale to zero for restart + tscale = 0.0; // set tscale to zero for restart p0_set = 1; v0_set = 1; e0_set = 1; @@ -828,43 +793,43 @@ void FixMSST::restart(char *buf) int FixMSST::modify_param(int narg, char **arg) { - if (strcmp(arg[0],"temp") == 0) { - if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (strcmp(arg[0], "temp") == 0) { + if (narg < 2) error->all(FLERR, "Illegal fix_modify command"); if (tflag) { modify->delete_compute(id_temp); tflag = 0; } - delete [] id_temp; + delete[] id_temp; id_temp = utils::strdup(arg[1]); int icompute = modify->find_compute(id_temp); - if (icompute < 0) - error->all(FLERR,"Could not find fix_modify temperature ID"); + if (icompute < 0) error->all(FLERR, "Could not find fix_modify temperature ID"); temperature = modify->compute[icompute]; if (temperature->tempflag == 0) - error->all(FLERR,"Fix_modify temperature ID does not " + error->all(FLERR, + "Fix_modify temperature ID does not " "compute temperature"); if (temperature->igroup != 0 && comm->me == 0) - error->warning(FLERR,"Temperature for MSST is not for group all"); + error->warning(FLERR, "Temperature for MSST is not for group all"); return 2; - } else if (strcmp(arg[0],"press") == 0) { - if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + } else if (strcmp(arg[0], "press") == 0) { + if (narg < 2) error->all(FLERR, "Illegal fix_modify command"); if (pflag) { modify->delete_compute(id_press); pflag = 0; } - delete [] id_press; + delete[] id_press; id_press = utils::strdup(arg[1]); int icompute = modify->find_compute(id_press); - if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID"); + if (icompute < 0) error->all(FLERR, "Could not find fix_modify pressure ID"); pressure = modify->compute[icompute]; if (pressure->pressflag == 0) - error->all(FLERR,"Fix_modify pressure ID does not compute pressure"); + error->all(FLERR, "Fix_modify pressure ID does not compute pressure"); return 2; } return 0; @@ -887,10 +852,9 @@ double FixMSST::compute_scalar() i = direction; energy = qmass * omega[i] * omega[i] / (2.0 * total_mass) * mvv2e; - energy -= 0.5 * total_mass * velocity * velocity * - (1.0 - volume/ v0) * - (1.0 - volume/ v0) * mvv2e; - energy -= p0 * ( v0 - volume ) / nktv2p; + energy -= + 0.5 * total_mass * velocity * velocity * (1.0 - volume / v0) * (1.0 - volume / v0) * mvv2e; + energy -= p0 * (v0 - volume) / nktv2p; // subtract off precomputed TS_int integral value // TS_int = 0 for non DFTB calculations @@ -938,8 +902,7 @@ double FixMSST::compute_hugoniot() v = compute_vol(); - dhugo = (0.5 * (p + p0 ) * ( v0 - v)) / - force->nktv2p + e0 - e; + dhugo = (0.5 * (p + p0) * (v0 - v)) / force->nktv2p + e0 - e; dhugo /= temperature->dof * force->boltz; return dhugo; @@ -964,8 +927,7 @@ double FixMSST::compute_rayleigh() v = compute_vol(); drayleigh = p - p0 - - total_mass * velocity * velocity * force->mvv2e * - (1.0 - v / v0 ) * force->nktv2p / v0; + total_mass * velocity * velocity * force->mvv2e * (1.0 - v / v0) * force->nktv2p / v0; return drayleigh; } @@ -978,7 +940,7 @@ double FixMSST::compute_rayleigh() double FixMSST::compute_lagrangian_speed() { double v = compute_vol(); - return velocity*(1.0-v/v0); + return velocity * (1.0 - v / v0); } /* ---------------------------------------------------------------------- @@ -988,7 +950,7 @@ double FixMSST::compute_lagrangian_speed() double FixMSST::compute_lagrangian_position() { - return lagrangian_position; + return lagrangian_position; } /* ---------------------------------------------------------------------- */ @@ -997,11 +959,11 @@ double FixMSST::compute_etotal() { if (!pe) return 0.0; - double epot,ekin,etot; + double epot, ekin, etot; epot = pe->compute_scalar(); ekin = temperature->compute_scalar(); ekin *= 0.5 * temperature->dof * force->boltz; - etot = epot+ekin; + etot = epot + ekin; return etot; } @@ -1028,12 +990,10 @@ double FixMSST::compute_vsum() double t = 0.0; for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { - t += (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]) ; - } + if (mask[i] & groupbit) { t += (v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]); } } - MPI_Allreduce(&t,&vsum,1,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&t, &vsum, 1, MPI_DOUBLE, MPI_SUM, world); return vsum; } @@ -1043,6 +1003,6 @@ double FixMSST::compute_vsum() double FixMSST::memory_usage() { - double bytes = 3*atom->nmax * sizeof(double); + double bytes = 3 * atom->nmax * sizeof(double); return bytes; } diff --git a/src/SHOCK/fix_msst.h b/src/SHOCK/fix_msst.h index 8cd3f79a89..c7d4983dc4 100644 --- a/src/SHOCK/fix_msst.h +++ b/src/SHOCK/fix_msst.h @@ -64,9 +64,8 @@ class FixMSST : public Fix { double **old_velocity; // saved velocities - int kspace_flag; // 1 if KSpace invoked, 0 if not - int nrigid; // number of rigid fixes - int *rfix; // indices of rigid fixes + int kspace_flag; // 1 if KSpace invoked, 0 if not + std::vector rfix; // indices of rigid fixes char *id_temp, *id_press; // strings with identifiers of char *id_pe; // created computes diff --git a/src/fix_press_berendsen.cpp b/src/fix_press_berendsen.cpp index 05e523abae..40dcdeeb10 100644 --- a/src/fix_press_berendsen.cpp +++ b/src/fix_press_berendsen.cpp @@ -223,17 +223,12 @@ FixPressBerendsen::FixPressBerendsen(LAMMPS *lmp, int narg, char **arg) : id_press = utils::strdup(std::string(id) + "_press"); modify->add_compute(fmt::format("{} all pressure {}", id_press, id_temp)); pflag = 1; - - nrigid = 0; - rfix = nullptr; } /* ---------------------------------------------------------------------- */ FixPressBerendsen::~FixPressBerendsen() { - delete[] rfix; - // delete temperature and pressure if fix created them if (tflag) modify->delete_compute(id_temp); @@ -291,20 +286,10 @@ void FixPressBerendsen::init() kspace_flag = 0; // detect if any rigid fixes exist so rigid bodies move when box is remapped - // rfix[] = indices to each fix rigid - delete[] rfix; - nrigid = 0; - rfix = nullptr; - - for (const auto &ifix : modify->get_fix_list()) - if (ifix->rigid_flag) nrigid++; - if (nrigid > 0) { - rfix = new Fix *[nrigid]; - nrigid = 0; - for (auto &ifix : modify->get_fix_list()) - if (ifix->rigid_flag) rfix[nrigid++] = ifix; - } + rfix.clear(); + for (auto &ifix : modify->get_fix_list()) + if (ifix->rigid_flag) rfix.push_back(ifix); } /* ---------------------------------------------------------------------- @@ -409,8 +394,7 @@ void FixPressBerendsen::remap() if (mask[i] & groupbit) domain->x2lamda(x[i], x[i]); } - if (nrigid) - for (i = 0; i < nrigid; i++) rfix[i]->deform(0); + for (auto &ifix : rfix) ifix->deform(0); // reset global and local box to new size/shape @@ -436,8 +420,7 @@ void FixPressBerendsen::remap() if (mask[i] & groupbit) domain->lamda2x(x[i], x[i]); } - if (nrigid) - for (i = 0; i < nrigid; i++) rfix[i]->deform(1); + for (auto &ifix : rfix) ifix->deform(1); } /* ---------------------------------------------------------------------- */ diff --git a/src/fix_press_berendsen.h b/src/fix_press_berendsen.h index 9e83533746..ddbe31e7ed 100644 --- a/src/fix_press_berendsen.h +++ b/src/fix_press_berendsen.h @@ -44,9 +44,8 @@ class FixPressBerendsen : public Fix { double p_period[3], p_target[3]; double p_current[3], dilation[3]; double factor[3]; - int kspace_flag; // 1 if KSpace invoked, 0 if not - int nrigid; // number of rigid fixes - class Fix **rfix; // indices of rigid fixes + int kspace_flag; // 1 if KSpace invoked, 0 if not + std::vector rfix; // indices of rigid fixes char *id_temp, *id_press; class Compute *temperature, *pressure; diff --git a/src/fix_press_langevin.cpp b/src/fix_press_langevin.cpp index 2f6e765cd5..752f826dfe 100644 --- a/src/fix_press_langevin.cpp +++ b/src/fix_press_langevin.cpp @@ -376,9 +376,6 @@ FixPressLangevin::FixPressLangevin(LAMMPS *lmp, int narg, char **arg) : (1.0 + p_alpha[i] * update->dt / 2.0 / p_mass[i]); gjfb[i] = 1. / (1.0 + p_alpha[i] * update->dt / 2.0 / p_mass[i]); } - - nrigid = 0; - rfix = nullptr; } /* ---------------------------------------------------------------------- */ @@ -386,7 +383,6 @@ FixPressLangevin::FixPressLangevin(LAMMPS *lmp, int narg, char **arg) : FixPressLangevin::~FixPressLangevin() { delete random; - delete[] rfix; delete irregular; // delete temperature and pressure if fix created them @@ -437,20 +433,10 @@ void FixPressLangevin::init() kspace_flag = 0; // detect if any rigid fixes exist so rigid bodies move when box is remapped - // rfix[] = indices to each fix rigid - delete[] rfix; - nrigid = 0; - rfix = nullptr; - - for (const auto &ifix : modify->get_fix_list()) - if (ifix->rigid_flag) nrigid++; - if (nrigid > 0) { - rfix = new Fix *[nrigid]; - nrigid = 0; - for (auto &ifix : modify->get_fix_list()) - if (ifix->rigid_flag) rfix[nrigid++] = ifix; - } + rfix.clear(); + for (auto &ifix : modify->get_fix_list()) + if (ifix->rigid_flag) rfix.push_back(ifix); // Nullifies piston derivatives and forces so that it is not integrated at // the start of a second run. @@ -680,8 +666,7 @@ void FixPressLangevin::remap() if (mask[i] & groupbit) domain->x2lamda(x[i], x[i]); } - if (nrigid) - for (i = 0; i < nrigid; i++) rfix[i]->deform(0); + for (auto &ifix : rfix) ifix->deform(0); // reset global and local box to new size/shape @@ -719,8 +704,7 @@ void FixPressLangevin::remap() if (mask[i] & groupbit) domain->lamda2x(x[i], x[i]); } - if (nrigid) - for (i = 0; i < nrigid; i++) rfix[i]->deform(1); + for (auto &ifix : rfix) ifix->deform(1); } /* ---------------------------------------------------------------------- diff --git a/src/fix_press_langevin.h b/src/fix_press_langevin.h index 868993b1f4..498f9443a7 100644 --- a/src/fix_press_langevin.h +++ b/src/fix_press_langevin.h @@ -51,9 +51,8 @@ class FixPressLangevin : public Fix { double p_deriv[6], dilation[6]; double f_piston[6], f_old_piston[6]; double gjfa[6], gjfb[6], fran[6]; - int kspace_flag; // 1 if KSpace invoked, 0 if not - int nrigid; // number of rigid fixes - class Fix **rfix; // list of rigid fixes + int kspace_flag; // 1 if KSpace invoked, 0 if not + std::vector rfix; // indices of rigid fixes char *id_temp, *id_press; class Compute *temperature, *pressure;