ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

modernize

Browse Source
This commit is contained in:
Axel Kohlmeyer
2023-11-11 12:07:05 -05:00
parent e6524b59fa
commit be6fcaa77f
14 changed files with 328 additions and 507 deletions
Download Patch File Download Diff File Expand all files Collapse all files

94
src/BOCS/fix_bocs.cpp
View File

@ -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<ComputePressureBocs *>(pressure))->send_cg_info(p_basis_type, N_p_match,
p_match_coeffs, N_mol, vavg);
if (p_match_flag) {
auto bocspress = dynamic_cast<ComputePressureBocs *>(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<ComputePressureBocs *>(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;
}

3
src/BOCS/fix_bocs.h
View File

@ -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<Fix *> rfix; // list of rigid fixes
char *id_dilate; // group name to dilate
class Irregular *irregular; // for migrating atoms after box flips

129
src/DRUDE/fix_tgnh_drude.cpp
View File

@ -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<FixDrude *>(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<FixDrude *>(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;
}

3
src/DRUDE/fix_tgnh_drude.h
View File

@ -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<Fix *> rfix; // indices of rigid fixes
class Irregular *irregular; // for migrating atoms after box flips
int nlevels_respa;

76
src/EXTRA-FIX/fix_npt_cauchy.cpp
View File

@ -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;
}

3
src/EXTRA-FIX/fix_npt_cauchy.h
View File

@ -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<Fix *> rfix; // indices of rigid fixes
char *id_dilate; // group name to dilate
class Irregular *irregular; // for migrating atoms after box flips

38
src/QTB/fix_qbmsst.cpp
View File

@ -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);

3
src/QTB/fix_qbmsst.h
View File

@ -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<Fix *> 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

418
src/SHOCK/fix_msst.cpp
View File

@ -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<FixExternal *>(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;
}

5
src/SHOCK/fix_msst.h
View File

@ -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<Fix *> rfix; // indices of rigid fixes
char *id_temp, *id_press; // strings with identifiers of
char *id_pe; // created computes

27
src/fix_press_berendsen.cpp
View File

@ -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);
}
/* ---------------------------------------------------------------------- */

5
src/fix_press_berendsen.h
View File

@ -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<Fix *> rfix; // indices of rigid fixes
char *id_temp, *id_press;
class Compute *temperature, *pressure;

26
src/fix_press_langevin.cpp
View File

@ -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);
}
/* ----------------------------------------------------------------------

5
src/fix_press_langevin.h
View File

@ -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<Fix *> rfix; // indices of rigid fixes
char *id_temp, *id_press;
class Compute *temperature, *pressure;