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

Merge pull request #3244 from akohlmey/shake-with-minimize

Browse Source 
Enable use of fix shake or fix rattle during minimization
This commit is contained in:
Axel Kohlmeyer
2022-08-17 18:02:40 -04:00
committed by GitHub
parent 20b87d3cd4 daa6b78c43
commit f199e45762
9 changed files with 448 additions and 256 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
src/KOKKOS/fix_shake_kokkos.cpp
View File

@ -12,12 +12,8 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <mpi.h>
#include <cmath>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include "fix_shake_kokkos.h"
#include "fix_rattle.h"
#include "atom_kokkos.h"
#include "atom_vec.h"
@ -38,6 +34,9 @@
#include "kokkos.h"
#include "atom_masks.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
@ -181,6 +180,16 @@ void FixShakeKokkos<DeviceType>::init()
}
/* ----------------------------------------------------------------------
run setup for minimization.
------------------------------------------------------------------------- */
template<class DeviceType>
void FixShakeKokkos<DeviceType>::min_setup(int /*vflag*/)
{
error->all(FLERR, "Cannot yet use fix {} during minimization with Kokkos", style);
}
/* ----------------------------------------------------------------------
build list of SHAKE clusters to constrain
if one or more atoms in cluster are on this proc,
@ -193,6 +202,7 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
// local copies of atom quantities
// used by SHAKE until next re-neighboring
ebond = 0.0;
x = atom->x;
v = atom->v;
f = atom->f;
@ -292,7 +302,7 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
if (h_error_flag() == 1) {
error->one(FLERR,"Shake atoms missing on proc "
"{} at step {}",me,update->ntimestep);
"{} at step {}",comm->me,update->ntimestep);
}
}
@ -303,6 +313,7 @@ void FixShakeKokkos<DeviceType>::pre_neighbor()
template<class DeviceType>
void FixShakeKokkos<DeviceType>::post_force(int vflag)
{
ebond = 0.0;
copymode = 1;
d_x = atomKK->k_x.view<DeviceType>();
@ -341,7 +352,7 @@ void FixShakeKokkos<DeviceType>::post_force(int vflag)
// communicate results if necessary
unconstrained_update();
if (nprocs > 1) comm->forward_comm(this);
if (comm->nprocs > 1) comm->forward_comm(this);
k_xshake.sync<DeviceType>();
// virial setup
@ -1702,7 +1713,7 @@ void FixShakeKokkos<DeviceType>::correct_coordinates(int vflag) {
double **xtmp = xshake;
xshake = x;
if (nprocs > 1) {
if (comm->nprocs > 1) {
forward_comm_device = 0;
comm->forward_comm(this);
forward_comm_device = 1;

1
src/KOKKOS/fix_shake_kokkos.h
View File

@ -51,6 +51,7 @@ class FixShakeKokkos : public FixShake, public KokkosBase {
FixShakeKokkos(class LAMMPS *, int, char **);
~FixShakeKokkos() override;
void init() override;
void min_setup(int) override;
void pre_neighbor() override;
void post_force(int) override;

104
src/RIGID/fix_rattle.cpp
View File

@ -81,7 +81,7 @@ FixRattle::~FixRattle()
{
memory->destroy(vp);
if (RATTLE_DEBUG) {
#if RATTLE_DEBUG
// communicate maximum distance error
@ -91,13 +91,11 @@ FixRattle::~FixRattle()
MPI_Reduce(&derr_max, &global_derr_max, 1 , MPI_DOUBLE, MPI_MAX, 0, world);
MPI_Reduce(&verr_max, &global_verr_max, 1 , MPI_DOUBLE, MPI_MAX, 0, world);
MPI_Comm_rank (world, &npid); // Find out process rank
if (npid == 0 && screen) {
if (comm->me == 0 && screen) {
fprintf(screen, "RATTLE: Maximum overall relative position error ( (r_ij-d_ij)/d_ij ): %.10g\n", global_derr_max);
fprintf(screen, "RATTLE: Maximum overall absolute velocity error (r_ij * v_ij): %.10g\n", global_verr_max);
}
}
#endif
}
/* ---------------------------------------------------------------------- */
@ -110,7 +108,10 @@ int FixRattle::setmask()
mask |= POST_FORCE_RESPA;
mask |= FINAL_INTEGRATE;
mask |= FINAL_INTEGRATE_RESPA;
if (RATTLE_DEBUG) mask |= END_OF_STEP;
mask |= MIN_POST_FORCE;
#if RATTLE_DEBUG
mask |= END_OF_STEP;
#endif
return mask;
}
@ -156,7 +157,7 @@ void FixRattle::post_force(int vflag)
// communicate the unconstrained velocities
if (nprocs > 1) {
if (comm->nprocs > 1) {
comm_mode = VP;
comm->forward_comm(this);
}
@ -188,7 +189,7 @@ void FixRattle::post_force_respa(int vflag, int ilevel, int /*iloop*/)
// communicate the unconstrained velocities
if (nprocs > 1) {
if (comm->nprocs > 1) {
comm_mode = VP;
comm->forward_comm(this);
}
@ -718,7 +719,7 @@ void FixRattle::unpack_forward_comm(int n, int first, double *buf)
void FixRattle::shake_end_of_step(int vflag) {
if (nprocs > 1) {
if (comm->nprocs > 1) {
comm_mode = V;
comm->forward_comm(this);
}
@ -738,7 +739,6 @@ void FixRattle::correct_coordinates(int vflag) {
FixShake::correct_coordinates(vflag);
}
/* ----------------------------------------------------------------------
Remove the velocity component along any bond.
------------------------------------------------------------------------- */
@ -759,7 +759,7 @@ void FixRattle::correct_velocities() {
// communicate the unconstrained velocities
if (nprocs > 1) {
if (comm->nprocs > 1) {
comm_mode = VP;
comm->forward_comm(this);
}
@ -769,14 +769,13 @@ void FixRattle::correct_velocities() {
int m;
for (int i = 0; i < nlist; i++) {
m = list[i];
if (shake_flag[m] == 2) vrattle2(m);
else if (shake_flag[m] == 3) vrattle3(m);
else if (shake_flag[m] == 4) vrattle4(m);
else vrattle3angle(m);
if (shake_flag[m] == 2) vrattle2(m);
else if (shake_flag[m] == 3) vrattle3(m);
else if (shake_flag[m] == 4) vrattle4(m);
else vrattle3angle(m);
}
}
/* ----------------------------------------------------------------------
DEBUGGING methods
The functions below allow you to check whether the
@ -788,16 +787,16 @@ void FixRattle::correct_velocities() {
void FixRattle::end_of_step()
{
if (nprocs > 1) {
comm_mode = V;
comm->forward_comm(this);
if (comm->nprocs > 1) {
comm_mode = V;
comm->forward_comm(this);
}
if (!check_constraints(v, RATTLE_TEST_POS, RATTLE_TEST_VEL) && RATTLE_RAISE_ERROR) {
#if RATTLE_RAISE_ERROR
if (!check_constraints(v, RATTLE_TEST_POS, RATTLE_TEST_VEL))
error->one(FLERR, "Rattle failed ");
}
#endif
}
/* ---------------------------------------------------------------------- */
bool FixRattle::check_constraints(double **v, bool checkr, bool checkv)
@ -807,12 +806,11 @@ bool FixRattle::check_constraints(double **v, bool checkr, bool checkv)
int i=0;
while (i < nlist && ret) {
m = list[i];
if (shake_flag[m] == 2) ret = check2(v,m,checkr,checkv);
else if (shake_flag[m] == 3) ret = check3(v,m,checkr,checkv);
else if (shake_flag[m] == 4) ret = check4(v,m,checkr,checkv);
else ret = check3angle(v,m,checkr,checkv);
if (shake_flag[m] == 2) ret = check2(v,m,checkr,checkv);
else if (shake_flag[m] == 3) ret = check3(v,m,checkr,checkv);
else if (shake_flag[m] == 4) ret = check4(v,m,checkr,checkv);
else ret = check3angle(v,m,checkr,checkv);
i++;
if (!RATTLE_RAISE_ERROR) ret = true;
}
return ret;
}
@ -834,14 +832,10 @@ bool FixRattle::check2(double **v, int m, bool checkr, bool checkv)
MathExtra::sub3(v[i1],v[i0],v01);
stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol));
if (!stat)
error->one(FLERR,"Coordinate constraints are not satisfied "
"up to desired tolerance ");
if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol));
if (!stat)
error->one(FLERR,"Velocity constraints are not satisfied "
"up to desired tolerance ");
if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance ");
return stat;
}
@ -871,15 +865,11 @@ bool FixRattle::check3(double **v, int m, bool checkr, bool checkv)
stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol ||
fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol));
if (!stat)
error->one(FLERR,"Coordinate constraints are not satisfied "
"up to desired tolerance ");
if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol ||
fabs(MathExtra::dot3(r02,v02)) > tol));
if (!stat)
error->one(FLERR,"Velocity constraints are not satisfied "
"up to desired tolerance ");
if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance ");
return stat;
}
@ -914,16 +904,12 @@ bool FixRattle::check4(double **v, int m, bool checkr, bool checkv)
stat = !(checkr && (fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1) > tol ||
fabs(sqrt(MathExtra::dot3(r02,r02))-bond2) > tol ||
fabs(sqrt(MathExtra::dot3(r03,r03))-bond3) > tol));
if (!stat)
error->one(FLERR,"Coordinate constraints are not satisfied "
"up to desired tolerance ");
if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
stat = !(checkv && (fabs(MathExtra::dot3(r01,v01)) > tol ||
fabs(MathExtra::dot3(r02,v02)) > tol ||
fabs(MathExtra::dot3(r03,v03)) > tol));
if (!stat)
error->one(FLERR,"Velocity constraints are not satisfied "
"up to desired tolerance ");
if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance ");
return stat;
}
@ -954,25 +940,19 @@ bool FixRattle::check3angle(double **v, int m, bool checkr, bool checkv)
MathExtra::sub3(v[i2],v[i0],v02);
MathExtra::sub3(v[i2],v[i1],v12);
double db1 = fabs(sqrt(MathExtra::dot3(r01,r01)) - bond1);
double db2 = fabs(sqrt(MathExtra::dot3(r02,r02))-bond2);
double db12 = fabs(sqrt(MathExtra::dot3(r12,r12))-bond12);
stat = !(checkr && (db1 > tol ||
db2 > tol ||
db12 > tol));
stat = !(checkr && (db1 > tol || db2 > tol || db12 > tol));
if (derr_max < db1/bond1) derr_max = db1/bond1;
if (derr_max < db2/bond2) derr_max = db2/bond2;
if (derr_max < db12/bond12) derr_max = db12/bond12;
if (!stat && RATTLE_RAISE_ERROR)
error->one(FLERR,"Coordinate constraints are not satisfied "
"up to desired tolerance ");
#if RATTLE_RAISE_ERROR
if (!stat) error->one(FLERR,"Coordinate constraints are not satisfied up to desired tolerance ");
#endif
double dv1 = fabs(MathExtra::dot3(r01,v01));
double dv2 = fabs(MathExtra::dot3(r02,v02));
@ -982,16 +962,10 @@ bool FixRattle::check3angle(double **v, int m, bool checkr, bool checkv)
if (verr_max < dv2) verr_max = dv2;
if (verr_max < dv12) verr_max = dv12;
stat = !(checkv && (dv1 > tol || dv2 > tol || dv12> tol));
stat = !(checkv && (dv1 > tol ||
dv2 > tol ||
dv12> tol));
if (!stat && RATTLE_RAISE_ERROR)
error->one(FLERR,"Velocity constraints are not satisfied "
"up to desired tolerance!");
#if RATTLE_RAISE_ERROR
if (!stat) error->one(FLERR,"Velocity constraints are not satisfied up to desired tolerance!");
#endif
return stat;
}

347
src/RIGID/fix_shake.cpp
View File

@ -41,8 +41,8 @@ using namespace MathConst;
#define RVOUS 1 // 0 for irregular, 1 for all2all
#define BIG 1.0e20
#define MASSDELTA 0.1
static constexpr double BIG = 1.0e20;
static constexpr double MASSDELTA = 0.1;
/* ---------------------------------------------------------------------- */
@ -52,26 +52,31 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
loop_respa(nullptr), step_respa(nullptr), x(nullptr), v(nullptr), f(nullptr), ftmp(nullptr),
vtmp(nullptr), mass(nullptr), rmass(nullptr), type(nullptr), shake_flag(nullptr),
shake_atom(nullptr), shake_type(nullptr), xshake(nullptr), nshake(nullptr), list(nullptr),
b_count(nullptr), b_count_all(nullptr), b_atom(nullptr), b_atom_all(nullptr), b_ave(nullptr), b_max(nullptr),
b_min(nullptr), b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr),
b_count(nullptr), b_count_all(nullptr), b_ave(nullptr), b_max(nullptr), b_min(nullptr),
b_ave_all(nullptr), b_max_all(nullptr), b_min_all(nullptr), a_count(nullptr),
a_count_all(nullptr), a_ave(nullptr), a_max(nullptr), a_min(nullptr), a_ave_all(nullptr),
a_max_all(nullptr), a_min_all(nullptr), atommols(nullptr), onemols(nullptr)
{
MPI_Comm_rank(world,&me);
MPI_Comm_size(world,&nprocs);
energy_global_flag = energy_peratom_flag = 1;
virial_global_flag = virial_peratom_flag = 1;
thermo_virial = 1;
thermo_energy = thermo_virial = 1;
create_attribute = 1;
dof_flag = 1;
scalar_flag = 1;
stores_ids = 1;
centroidstressflag = CENTROID_AVAIL;
next_output = -1;
// to avoid uninitialized access
vflag_post_force = 0;
eflag_pre_reverse = 0;
ebond = 0.0;
// error check
molecular = atom->molecular;
if (molecular == Atom::ATOMIC)
error->all(FLERR,"Cannot use fix shake with non-molecular system");
error->all(FLERR,"Cannot use fix {} with non-molecular system", style);
// perform initial allocation of atom-based arrays
// register with Atom class
@ -92,8 +97,9 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
comm_forward = 3;
// parse SHAKE args
auto mystyle = fmt::format("fix {}",style);
if (narg < 8) error->all(FLERR,"Illegal fix shake command");
if (narg < 8) utils::missing_cmd_args(FLERR,mystyle, error);
tolerance = utils::numeric(FLERR,arg[3],false,lmp);
max_iter = utils::inumeric(FLERR,arg[4],false,lmp);
@ -133,49 +139,55 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
else if (mode == 'b') {
int i = utils::inumeric(FLERR,arg[next],false,lmp);
if (i < 1 || i > atom->nbondtypes)
error->all(FLERR,"Invalid bond type index for fix shake");
error->all(FLERR,"Invalid bond type index for {}", mystyle);
bond_flag[i] = 1;
} else if (mode == 'a') {
int i = utils::inumeric(FLERR,arg[next],false,lmp);
if (i < 1 || i > atom->nangletypes)
error->all(FLERR,"Invalid angle type index for fix shake");
error->all(FLERR,"Invalid angle type index for {}", mystyle);
angle_flag[i] = 1;
} else if (mode == 't') {
int i = utils::inumeric(FLERR,arg[next],false,lmp);
if (i < 1 || i > atom->ntypes)
error->all(FLERR,"Invalid atom type index for fix shake");
error->all(FLERR,"Invalid atom type index for {}", mystyle);
type_flag[i] = 1;
} else if (mode == 'm') {
double massone = utils::numeric(FLERR,arg[next],false,lmp);
if (massone == 0.0) error->all(FLERR,"Invalid atom mass for fix shake");
if (massone == 0.0) error->all(FLERR,"Invalid atom mass for {}", mystyle);
if (nmass == atom->ntypes)
error->all(FLERR,"Too many masses for fix shake");
error->all(FLERR,"Too many masses for {}", mystyle);
mass_list[nmass++] = massone;
} else error->all(FLERR,"Illegal fix shake command");
} else error->all(FLERR,"Unknown {} command option: {}", mystyle, arg[next]);
next++;
}
// parse optional args
onemols = nullptr;
kbond = 1.0e6*force->boltz;
int iarg = next;
while (iarg < narg) {
if (strcmp(arg[next],"mol") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix shake command");
if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR,mystyle+" mol",error);
int imol = atom->find_molecule(arg[iarg+1]);
if (imol == -1)
error->all(FLERR,"Molecule template ID for fix shake does not exist");
if (atom->molecules[imol]->nset > 1 && comm->me == 0)
error->warning(FLERR,"Molecule template for fix shake has multiple molecules");
error->all(FLERR,"Molecule template ID {} for {} does not exist", mystyle, arg[iarg+1]);
if ((atom->molecules[imol]->nset > 1) && (comm->me == 0))
error->warning(FLERR,"Molecule template for {} has multiple molecules", mystyle);
onemols = &atom->molecules[imol];
nmol = onemols[0]->nset;
iarg += 2;
} else error->all(FLERR,"Illegal fix shake command");
} else if (strcmp(arg[iarg],"kbond") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR,mystyle+" kbond",error);
kbond = utils::numeric(FLERR, arg[iarg+1], false, lmp);
if (kbond < 0) error->all(FLERR,"Illegal {} kbond value {}. Must be >= 0.0", mystyle, kbond);
iarg += 2;
} else error->all(FLERR,"Unknown {} command option: {}", mystyle, arg[iarg]);
}
// error check for Molecule template
@ -183,7 +195,7 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
if (onemols) {
for (int i = 0; i < nmol; i++)
if (onemols[i]->shakeflag == 0)
error->all(FLERR,"Fix shake molecule template must have shake info");
error->all(FLERR,"Fix {} molecule template must have shake info", style);
}
// allocate bond and angle distance arrays, indexed from 1 to n
@ -197,8 +209,6 @@ FixShake::FixShake(LAMMPS *lmp, int narg, char **arg) :
int nb = atom->nbondtypes + 1;
b_count = new int[nb];
b_count_all = new int[nb];
b_atom = new int[nb];
b_atom_all = new int[nb];
b_ave = new double[nb];
b_ave_all = new double[nb];
b_max = new double[nb];
@ -291,8 +301,6 @@ FixShake::~FixShake()
if (output_every) {
delete[] b_count;
delete[] b_count_all;
delete[] b_atom;
delete[] b_atom_all;
delete[] b_ave;
delete[] b_ave_all;
delete[] b_max;
@ -321,6 +329,9 @@ int FixShake::setmask()
mask |= PRE_NEIGHBOR;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= MIN_PRE_REVERSE;
mask |= MIN_POST_FORCE;
mask |= POST_RUN;
return mask;
}
@ -335,28 +346,24 @@ void FixShake::init()
double rsq,angle;
// error if more than one shake fix
auto pattern = fmt::format("^{}",style);
int count = 0;
for (i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"shake") == 0) count++;
if (count > 1) error->all(FLERR,"More than one fix shake");
if (modify->get_fix_by_style(pattern).size() > 1)
error->all(FLERR,"More than one fix {} instance",style);
// cannot use with minimization since SHAKE turns off bonds
// that should contribute to potential energy
if (update->whichflag == 2)
error->all(FLERR,"Fix shake cannot be used with minimization");
if ((comm->me == 0) && (update->whichflag == 2))
error->warning(FLERR,"Using fix {} with minimization.\n Substituting constraints with "
"harmonic restraint forces using kbond={:.4g}", style, kbond);
// error if npt,nph fix comes before shake fix
for (i = 0; i < modify->nfix; i++) {
if (strcmp(modify->fix[i]->style,"npt") == 0) break;
if (strcmp(modify->fix[i]->style,"nph") == 0) break;
}
if (i < modify->nfix) {
for (int j = i; j < modify->nfix; j++)
if (strcmp(modify->fix[j]->style,"shake") == 0)
error->all(FLERR,"Shake fix must come before NPT/NPH fix");
// error if a fix changing the box comes before shake fix
bool boxflag = false;
for (auto ifix : modify->get_fix_list()) {
if (boxflag && utils::strmatch(ifix->style,pattern))
error->all(FLERR,"Fix {} must come before any box changing fix", style);
if (ifix->box_change) boxflag = true;
}
// if rRESPA, find associated fix that must exist
@ -379,7 +386,7 @@ void FixShake::init()
// set equilibrium bond distances
if (force->bond == nullptr)
error->all(FLERR,"Bond potential must be defined for SHAKE");
error->all(FLERR,"Bond style must be defined for fix {}",style);
for (i = 1; i <= atom->nbondtypes; i++)
bond_distance[i] = force->bond->equilibrium_distance(i);
@ -390,7 +397,7 @@ void FixShake::init()
for (i = 1; i <= atom->nangletypes; i++) {
if (angle_flag[i] == 0) continue;
if (force->angle == nullptr)
error->all(FLERR,"Angle potential must be defined for SHAKE");
error->all(FLERR,"Angle style must be defined for fix {}",style);
// scan all atoms for a SHAKE angle cluster
// extract bond types for the 2 bonds in the cluster
@ -417,7 +424,7 @@ void FixShake::init()
// error check for any bond types that are not the same
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_MAX,world);
if (flag_all) error->all(FLERR,"Shake angles have different bond types");
if (flag_all) error->all(FLERR,"Fix {} angles have different bond types", style);
// insure all procs have bond types
@ -494,6 +501,23 @@ void FixShake::setup(int vflag)
shake_end_of_step(vflag);
}
/* ----------------------------------------------------------------------
during minimization fix SHAKE adds strong bond forces
------------------------------------------------------------------------- */
void FixShake::min_setup(int vflag)
{
pre_neighbor();
min_post_force(vflag);
}
/* --------------------------------------------------------------------- */
void FixShake::setup_pre_reverse(int eflag, int vflag)
{
min_pre_reverse(eflag,vflag);
}
/* ----------------------------------------------------------------------
build list of SHAKE clusters to constrain
if one or more atoms in cluster are on this proc,
@ -502,6 +526,7 @@ void FixShake::setup(int vflag)
void FixShake::pre_neighbor()
{
ebond = 0.0;
int atom1,atom2,atom3,atom4;
// local copies of atom quantities
@ -533,19 +558,17 @@ void FixShake::pre_neighbor()
atom1 = atom->map(shake_atom[i][0]);
atom2 = atom->map(shake_atom[i][1]);
if (atom1 == -1 || atom2 == -1)
error->one(FLERR,"Shake atoms {} {} missing on proc "
"{} at step {}",shake_atom[i][0],
shake_atom[i][1],me,update->ntimestep);
error->one(FLERR,"Shake atoms {} {} missing on proc {} at step {}",shake_atom[i][0],
shake_atom[i][1],comm->me,update->ntimestep);
if (i <= atom1 && i <= atom2) list[nlist++] = i;
} else if (shake_flag[i] % 2 == 1) {
atom1 = atom->map(shake_atom[i][0]);
atom2 = atom->map(shake_atom[i][1]);
atom3 = atom->map(shake_atom[i][2]);
if (atom1 == -1 || atom2 == -1 || atom3 == -1)
error->one(FLERR,"Shake atoms {} {} {} missing on proc "
"{} at step {}",shake_atom[i][0],
error->one(FLERR,"Shake atoms {} {} {} missing on proc {} at step {}",shake_atom[i][0],
shake_atom[i][1],shake_atom[i][2],
me,update->ntimestep);
comm->me,update->ntimestep);
if (i <= atom1 && i <= atom2 && i <= atom3) list[nlist++] = i;
} else {
atom1 = atom->map(shake_atom[i][0]);
@ -553,10 +576,9 @@ void FixShake::pre_neighbor()
atom3 = atom->map(shake_atom[i][2]);
atom4 = atom->map(shake_atom[i][3]);
if (atom1 == -1 || atom2 == -1 || atom3 == -1 || atom4 == -1)
error->one(FLERR,"Shake atoms {} {} {} {} missing on "
"proc {} at step {}",shake_atom[i][0],
error->one(FLERR,"Shake atoms {} {} {} {} missing on proc {} at step {}",shake_atom[i][0],
shake_atom[i][1],shake_atom[i][2],
shake_atom[i][3],me,update->ntimestep);
shake_atom[i][3],comm->me,update->ntimestep);
if (i <= atom1 && i <= atom2 && i <= atom3 && i <= atom4)
list[nlist++] = i;
}
@ -575,11 +597,13 @@ void FixShake::post_force(int vflag)
// communicate results if necessary
unconstrained_update();
if (nprocs > 1) comm->forward_comm(this);
if (comm->nprocs > 1) comm->forward_comm(this);
// virial setup
v_init(vflag);
int eflag = eflag_pre_reverse;
ev_init(eflag, vflag);
ebond = 0.0;
// loop over clusters to add constraint forces
@ -619,7 +643,7 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
// communicate results if necessary
unconstrained_update_respa(ilevel);
if (nprocs > 1) comm->forward_comm(this);
if (comm->nprocs > 1) comm->forward_comm(this);
// virial setup only needed on last iteration of innermost level
// and if pressure is requested
@ -644,6 +668,59 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
vflag_post_force = vflag;
}
/* ----------------------------------------------------------------------
store eflag so it can be used in min_post_force
------------------------------------------------------------------------- */
void FixShake::min_pre_reverse(int eflag, int /*vflag*/)
{
eflag_pre_reverse = eflag;
}
/* ----------------------------------------------------------------------
substitute shake constraints with very strong bonds
------------------------------------------------------------------------- */
void FixShake::min_post_force(int vflag)
{
if (output_every) {
bigint ntimestep = update->ntimestep;
if (next_output == ntimestep) stats();
next_output = ntimestep + output_every;
if (ntimestep % output_every != 0)
next_output = (ntimestep/output_every)*output_every + output_every;
} else next_output = -1;
int eflag = eflag_pre_reverse;
ev_init(eflag, vflag);
x = atom->x;
f = atom->f;
nlocal = atom->nlocal;
ebond = 0.0;
// loop over shake clusters to add restraint forces
for (int i = 0; i < nlist; i++) {
int m = list[i];
if (shake_flag[m] == 2) {
bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
} else if (shake_flag[m] == 3) {
bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
} else if (shake_flag[m] == 4) {
bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
bond_force(shake_atom[m][0], shake_atom[m][3], bond_distance[shake_type[m][2]]);
} else {
bond_force(shake_atom[m][0], shake_atom[m][1], bond_distance[shake_type[m][0]]);
bond_force(shake_atom[m][0], shake_atom[m][2], bond_distance[shake_type[m][1]]);
bond_force(shake_atom[m][1], shake_atom[m][2], angle_distance[shake_type[m][2]]);
}
}
}
/* ----------------------------------------------------------------------
count # of degrees-of-freedom removed by SHAKE for atoms in igroup
------------------------------------------------------------------------- */
@ -690,10 +767,7 @@ void FixShake::find_clusters()
tagint tagprev;
double massone;
if ((me == 0) && screen) {
if (!rattle) fputs("Finding SHAKE clusters ...\n",screen);
else fputs("Finding RATTLE clusters ...\n",screen);
}
if (comm->me == 0) utils::logmesg(lmp, "Finding {} clusters ...\n",utils::uppercase(style));
atommols = atom->avec->onemols;
tagint *tag = atom->tag;
@ -805,7 +879,7 @@ void FixShake::find_clusters()
}
MPI_Allreduce(&flag,&flag_all,1,MPI_INT,MPI_SUM,world);
if (flag_all) error->all(FLERR,"Did not find fix shake partner info");
if (flag_all) error->all(FLERR,"Did not find fix {} partner info", style);
// -----------------------------------------------------
// identify SHAKEable bonds
@ -1015,7 +1089,7 @@ void FixShake::find_clusters()
tmp = count4;
MPI_Allreduce(&tmp,&count4,1,MPI_INT,MPI_SUM,world);
if (me == 0) {
if (comm->me == 0) {
utils::logmesg(lmp,"{:>8} = # of size 2 clusters\n"
"{:>8} = # of size 3 clusters\n"
"{:>8} = # of size 4 clusters\n"
@ -1043,8 +1117,8 @@ void FixShake::atom_owners()
// one datum for each owned atom: datum = owning proc, atomID
for (int i = 0; i < nlocal; i++) {
proclist[i] = tag[i] % nprocs;
idbuf[i].me = me;
proclist[i] = tag[i] % comm->nprocs;
idbuf[i].me = comm->me;
idbuf[i].atomID = tag[i];
}
@ -1089,7 +1163,7 @@ void FixShake::partner_info(int *npartner, tagint **partner_tag,
// set values in 4 partner arrays for all partner atoms I own
// also setup input buf to rendezvous comm
// input datums = pair of bonded atoms where I do not own partner
// owning proc for each datum = partner_tag % nprocs
// owning proc for each datum = partner_tag % comm->nprocs
// datum: atomID = partner_tag (off-proc), partnerID = tag (on-proc)
// 4 values for my owned atom
@ -1127,7 +1201,7 @@ void FixShake::partner_info(int *npartner, tagint **partner_tag,
}
} else {
proclist[nsend] = partner_tag[i][j] % nprocs;
proclist[nsend] = partner_tag[i][j] % comm->nprocs;
inbuf[nsend].atomID = partner_tag[i][j];
inbuf[nsend].partnerID = tag[i];
inbuf[nsend].mask = mask[i];
@ -1217,7 +1291,7 @@ void FixShake::nshake_info(int *npartner, tagint **partner_tag,
// set partner_nshake for all partner atoms I own
// also setup input buf to rendezvous comm
// input datums = pair of bonded atoms where I do not own partner
// owning proc for each datum = partner_tag % nprocs
// owning proc for each datum = partner_tag % comm->nprocs
// datum: atomID = partner_tag (off-proc), partnerID = tag (on-proc)
// nshake value for my owned atom
@ -1231,7 +1305,7 @@ void FixShake::nshake_info(int *npartner, tagint **partner_tag,
if (m >= 0 && m < nlocal) {
partner_nshake[i][j] = nshake[m];
} else {
proclist[nsend] = partner_tag[i][j] % nprocs;
proclist[nsend] = partner_tag[i][j] % comm->nprocs;
inbuf[nsend].atomID = partner_tag[i][j];
inbuf[nsend].partnerID = tag[i];
inbuf[nsend].nshake = nshake[i];
@ -1295,7 +1369,7 @@ void FixShake::shake_info(int *npartner, tagint **partner_tag,
// set 3 shake arrays for all partner atoms I own
// also setup input buf to rendezvous comm
// input datums = partner atom where I do not own partner
// owning proc for each datum = partner_tag % nprocs
// owning proc for each datum = partner_tag % comm->nprocs
// datum: atomID = partner_tag (off-proc)
// values in 3 shake arrays
@ -1317,7 +1391,7 @@ void FixShake::shake_info(int *npartner, tagint **partner_tag,
shake_type[m][2] = shake_type[i][2];
} else {
proclist[nsend] = partner_tag[i][j] % nprocs;
proclist[nsend] = partner_tag[i][j] % comm->nprocs;
inbuf[nsend].atomID = partner_tag[i][j];
inbuf[nsend].shake_flag = shake_flag[i];
inbuf[nsend].shake_atom[0] = shake_atom[i][0];
@ -2451,13 +2525,66 @@ void FixShake::shake3angle(int m)
}
}
/* ----------------------------------------------------------------------
apply bond force for minimization
------------------------------------------------------------------------- */
void FixShake::bond_force(tagint id1, tagint id2, double length)
{
int i1 = atom->map(id1);
int i2 = atom->map(id2);
if ((i1 < 0) || (i2 < 0)) return;
// distance vec between atoms, with PBC
double delx = x[i1][0] - x[i2][0];
double dely = x[i1][1] - x[i2][1];
double delz = x[i1][2] - x[i2][2];
domain->minimum_image(delx, dely, delz);
// compute and apply force
const double r = sqrt(delx * delx + dely * dely + delz * delz);
const double dr = r - length;
const double rk = kbond * dr;
const double fbond = (r > 0.0) ? -2.0 * rk / r : 0.0;
const double eb = rk*dr;
int list[2];
int nlist = 0;
if (i1 < nlocal) {
f[i1][0] += delx * fbond;
f[i1][1] += dely * fbond;
f[i1][2] += delz * fbond;
list[nlist++] = i1;
ebond += 0.5*eb;
}
if (i2 < nlocal) {
f[i2][0] -= delx * fbond;
f[i2][1] -= dely * fbond;
f[i2][2] -= delz * fbond;
list[nlist++] = i2;
ebond += 0.5*eb;
}
if (evflag) {
double v[6];
v[0] = 0.5 * delx * delx * fbond;
v[1] = 0.5 * dely * dely * fbond;
v[2] = 0.5 * delz * delz * fbond;
v[3] = 0.5 * delx * dely * fbond;
v[4] = 0.5 * delx * delz * fbond;
v[5] = 0.5 * dely * delz * fbond;
ev_tally(nlist, list, 2.0, eb, v);
}
}
/* ----------------------------------------------------------------------
print-out bond & angle statistics
------------------------------------------------------------------------- */
void FixShake::stats()
{
int i,j,m,n,iatom,jatom,katom;
double delx,dely,delz;
double r,r1,r2,r3,angle;
@ -2466,13 +2593,12 @@ void FixShake::stats()
int nb = atom->nbondtypes + 1;
int na = atom->nangletypes + 1;
for (i = 0; i < nb; i++) {
for (int i = 0; i < nb; i++) {
b_count[i] = 0;
b_ave[i] = b_max[i] = 0.0;
b_min[i] = BIG;
b_atom[i] = -1;
}
for (i = 0; i < na; i++) {
for (int i = 0; i < na; i++) {
a_count[i] = 0;
a_ave[i] = a_max[i] = 0.0;
a_min[i] = BIG;
@ -2484,25 +2610,26 @@ void FixShake::stats()
double **x = atom->x;
int nlocal = atom->nlocal;
for (i = 0; i < nlocal; i++) {
if (shake_flag[i] == 0) continue;
for (int ii = 0; ii < nlist; ++ii) {
int i = list[ii];
int n = shake_flag[i];
if (n == 0) continue;
// bond stats
n = shake_flag[i];
if (n == 1) n = 3;
iatom = atom->map(shake_atom[i][0]);
for (j = 1; j < n; j++) {
jatom = atom->map(shake_atom[i][j]);
int iatom = atom->map(shake_atom[i][0]);
for (int j = 1; j < n; j++) {
int jatom = atom->map(shake_atom[i][j]);
if (jatom >= nlocal) continue;
delx = x[iatom][0] - x[jatom][0];
dely = x[iatom][1] - x[jatom][1];
delz = x[iatom][2] - x[jatom][2];
domain->minimum_image(delx,dely,delz);
r = sqrt(delx*delx + dely*dely + delz*delz);
m = shake_type[i][j-1];
r = sqrt(delx*delx + dely*dely + delz*delz);
int m = shake_type[i][j-1];
b_count[m]++;
b_atom[m] = n;
b_ave[m] += r;
b_max[m] = MAX(b_max[m],r);
b_min[m] = MIN(b_min[m],r);
@ -2511,9 +2638,13 @@ void FixShake::stats()
// angle stats
if (shake_flag[i] == 1) {
iatom = atom->map(shake_atom[i][0]);
jatom = atom->map(shake_atom[i][1]);
katom = atom->map(shake_atom[i][2]);
int iatom = atom->map(shake_atom[i][0]);
int jatom = atom->map(shake_atom[i][1]);
int katom = atom->map(shake_atom[i][2]);
int n = 0;
if (iatom < nlocal) ++n;
if (jatom < nlocal) ++n;
if (katom < nlocal) ++n;
delx = x[iatom][0] - x[jatom][0];
dely = x[iatom][1] - x[jatom][1];
@ -2535,9 +2666,9 @@ void FixShake::stats()
angle = acos((r1*r1 + r2*r2 - r3*r3) / (2.0*r1*r2));
angle *= 180.0/MY_PI;
m = shake_type[i][2];
a_count[m]++;
a_ave[m] += angle;
int m = shake_type[i][2];
a_count[m] += n;
a_ave[m] += n*angle;
a_max[m] = MAX(a_max[m],angle);
a_min[m] = MIN(a_min[m],angle);
}
@ -2546,7 +2677,6 @@ void FixShake::stats()
// sum across all procs
MPI_Allreduce(b_count,b_count_all,nb,MPI_INT,MPI_SUM,world);
MPI_Allreduce(b_atom,b_atom_all,nb,MPI_INT,MPI_MAX,world);
MPI_Allreduce(b_ave,b_ave_all,nb,MPI_DOUBLE,MPI_SUM,world);
MPI_Allreduce(b_max,b_max_all,nb,MPI_DOUBLE,MPI_MAX,world);
MPI_Allreduce(b_min,b_min_all,nb,MPI_DOUBLE,MPI_MIN,world);
@ -2558,16 +2688,17 @@ void FixShake::stats()
// print stats only for non-zero counts
if (me == 0) {
if (comm->me == 0) {
const int width = log10((MAX(MAX(1,nb),na)))+2;
auto mesg = fmt::format("SHAKE stats (type/ave/delta/count) on step {}\n", update->ntimestep);
for (i = 1; i < nb; i++) {
auto mesg = fmt::format("{} stats (type/ave/delta/count) on step {}\n",
utils::uppercase(style), update->ntimestep);
for (int i = 1; i < nb; i++) {
const auto bcnt = b_count_all[i];
if (bcnt)
mesg += fmt::format("Bond: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width,
b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt/b_atom_all[i]);
b_ave_all[i]/bcnt,b_max_all[i]-b_min_all[i],bcnt);
}
for (i = 1; i < na; i++) {
for (int i = 1; i < na; i++) {
const auto acnt = a_count_all[i];
if (acnt)
mesg += fmt::format("Angle: {:>{}d} {:<9.6} {:<11.6} {:>8d}\n",i,width,
@ -3005,6 +3136,26 @@ void *FixShake::extract(const char *str, int &dim)
return nullptr;
}
/* ----------------------------------------------------------------------
energy due to restraint forces
------------------------------------------------------------------------- */
double FixShake::compute_scalar()
{
double all;
MPI_Allreduce(&ebond, &all, 1, MPI_DOUBLE, MPI_SUM, world);
return all;
}
/* ----------------------------------------------------------------------
print shake stats at the end of a minimization
------------------------------------------------------------------------- */
void FixShake::post_run()
{
if ((update->whichflag == 2) && (output_every > 0)) stats();
}
/* ----------------------------------------------------------------------
add coordinate constraining forces
this method is called at the end of a timestep
@ -3109,7 +3260,7 @@ void FixShake::correct_coordinates(int vflag) {
double **xtmp = xshake;
xshake = x;
if (nprocs > 1) {
if (comm->nprocs > 1) {
comm->forward_comm(this);
}
xshake = xtmp;

24
src/RIGID/fix_shake.h
View File

@ -34,9 +34,14 @@ class FixShake : public Fix {
int setmask() override;
void init() override;
void setup(int) override;
void setup_pre_reverse(int, int) override;
void min_setup(int) override;
void pre_neighbor() override;
void post_force(int) override;
void post_force_respa(int, int, int) override;
void min_pre_reverse(int, int) override;
void min_post_force(int) override;
void post_run() override;
double memory_usage() override;
void grow_arrays(int) override;
@ -57,16 +62,17 @@ class FixShake : public Fix {
int dof(int) override;
void reset_dt() override;
void *extract(const char *, int &) override;
double compute_scalar() override;
protected:
int vflag_post_force; // store the vflag of last post_force call
int eflag_pre_reverse; // store the eflag of last pre_reverse call
int respa; // 0 = vel. Verlet, 1 = respa
int me, nprocs;
int rattle; // 0 = SHAKE, 1 = RATTLE
double tolerance; // SHAKE tolerance
int max_iter; // max # of SHAKE iterations
int output_every; // SHAKE stat output every so often
bigint next_output; // timestep for next output
int rattle; // 0 = SHAKE, 1 = RATTLE
double tolerance; // SHAKE tolerance
int max_iter; // max # of SHAKE iterations
int output_every; // SHAKE stat output every so often
bigint next_output; // timestep for next output
// settings from input command
int *bond_flag, *angle_flag; // bond/angle types to constrain
@ -76,6 +82,8 @@ class FixShake : public Fix {
int molecular; // copy of atom->molecular
double *bond_distance, *angle_distance; // constraint distances
double kbond; // force constant for restraint
double ebond; // energy of bond restraints
class FixRespa *fix_respa; // rRESPA fix needed by SHAKE
int nlevels_respa; // copies of needed rRESPA variables
@ -108,8 +116,7 @@ class FixShake : public Fix {
int nlist, maxlist; // size and max-size of list
// stat quantities
int *b_count, *b_count_all, *b_atom,
*b_atom_all; // counts for each bond type, atoms in bond cluster
int *b_count, *b_count_all; // counts for each bond type, atoms in bond cluster
double *b_ave, *b_max, *b_min; // ave/max/min dist for each bond type
double *b_ave_all, *b_max_all, *b_min_all; // MPI summing arrays
int *a_count, *a_count_all; // ditto for angle types
@ -133,6 +140,7 @@ class FixShake : public Fix {
void shake3(int);
void shake4(int);
void shake3angle(int);
void bond_force(tagint, tagint, double);
void stats();
int bondtype_findset(int, tagint, tagint, int);
int angletype_findset(int, tagint, tagint, int);

66
src/compute_pe_atom.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -13,30 +12,31 @@
------------------------------------------------------------------------- */
#include "compute_pe_atom.h"
#include <cstring>
#include "atom.h"
#include "update.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "atom.h"
#include "bond.h"
#include "comm.h"
#include "dihedral.h"
#include "error.h"
#include "force.h"
#include "improper.h"
#include "kspace.h"
#include "modify.h"
#include "memory.h"
#include "error.h"
#include "modify.h"
#include "pair.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
energy(nullptr)
Compute(lmp, narg, arg), energy(nullptr)
{
if (narg < 3) error->all(FLERR,"Illegal compute pe/atom command");
if (narg < 3) error->all(FLERR, "Illegal compute pe/atom command");
peratom_flag = 1;
size_peratom_cols = 0;
@ -56,14 +56,22 @@ ComputePEAtom::ComputePEAtom(LAMMPS *lmp, int narg, char **arg) :
fixflag = 0;
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg],"pair") == 0) pairflag = 1;
else if (strcmp(arg[iarg],"bond") == 0) bondflag = 1;
else if (strcmp(arg[iarg],"angle") == 0) angleflag = 1;
else if (strcmp(arg[iarg],"dihedral") == 0) dihedralflag = 1;
else if (strcmp(arg[iarg],"improper") == 0) improperflag = 1;
else if (strcmp(arg[iarg],"kspace") == 0) kspaceflag = 1;
else if (strcmp(arg[iarg],"fix") == 0) fixflag = 1;
else error->all(FLERR,"Illegal compute pe/atom command");
if (strcmp(arg[iarg], "pair") == 0)
pairflag = 1;
else if (strcmp(arg[iarg], "bond") == 0)
bondflag = 1;
else if (strcmp(arg[iarg], "angle") == 0)
angleflag = 1;
else if (strcmp(arg[iarg], "dihedral") == 0)
dihedralflag = 1;
else if (strcmp(arg[iarg], "improper") == 0)
improperflag = 1;
else if (strcmp(arg[iarg], "kspace") == 0)
kspaceflag = 1;
else if (strcmp(arg[iarg], "fix") == 0)
fixflag = 1;
else
error->all(FLERR, "Illegal compute pe/atom command");
iarg++;
}
}
@ -86,7 +94,7 @@ void ComputePEAtom::compute_peratom()
invoked_peratom = update->ntimestep;
if (update->eflag_atom != invoked_peratom)
error->all(FLERR,"Per-atom energy was not tallied on needed timestep");
error->all(FLERR, "Per-atom energy was not tallied on needed timestep");
// grow local energy array if necessary
// needs to be atom->nmax in length
@ -94,7 +102,7 @@ void ComputePEAtom::compute_peratom()
if (atom->nmax > nmax) {
memory->destroy(energy);
nmax = atom->nmax;
memory->create(energy,nmax,"pe/atom:energy");
memory->create(energy, nmax, "pe/atom:energy");
vector_atom = energy;
}
@ -153,13 +161,11 @@ void ComputePEAtom::compute_peratom()
// add in per-atom contributions from relevant fixes
// always only for owned atoms, not ghost
if (fixflag && modify->n_energy_atom)
modify->energy_atom(nlocal,energy);
if (fixflag && modify->n_energy_atom) modify->energy_atom(nlocal, energy);
// communicate ghost energy between neighbor procs
if (force->newton || (force->kspace && force->kspace->tip4pflag))
comm->reverse_comm(this);
if (force->newton || (force->kspace && force->kspace->tip4pflag)) comm->reverse_comm(this);
// zero energy of atoms not in group
// only do this after comm since ghost contributions must be included
@ -174,7 +180,7 @@ void ComputePEAtom::compute_peratom()
int ComputePEAtom::pack_reverse_comm(int n, int first, double *buf)
{
int i,m,last;
int i, m, last;
m = 0;
last = first + n;
@ -186,7 +192,7 @@ int ComputePEAtom::pack_reverse_comm(int n, int first, double *buf)
void ComputePEAtom::unpack_reverse_comm(int n, int *list, double *buf)
{
int i,j,m;
int i, j, m;
m = 0;
for (i = 0; i < n; i++) {
@ -201,6 +207,6 @@ void ComputePEAtom::unpack_reverse_comm(int n, int *list, double *buf)
double ComputePEAtom::memory_usage()
{
double bytes = (double)nmax * sizeof(double);
double bytes = (double) nmax * sizeof(double);
return bytes;
}

21
src/fix_restrain.cpp
View File

@ -19,17 +19,18 @@
#include "fix_restrain.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "force.h"
#include "update.h"
#include "domain.h"
#include "comm.h"
#include "respa.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "respa.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -271,14 +272,12 @@ void FixRestrain::restrain_bond(int m)
if (newton_bond) {
if (i2 == -1 || i2 >= nlocal) return;
if (i1 == -1)
error->one(FLERR,"Restrain atoms {} {} missing on "
"proc {} at step {}", ids[m][0],ids[m][1],
error->one(FLERR,"Restrain atoms {} {} missing on proc {} at step {}", ids[m][0],ids[m][1],
comm->me,update->ntimestep);
} else {
if ((i1 == -1 || i1 >= nlocal) && (i2 == -1 || i2 >= nlocal)) return;
if (i1 == -1 || i2 == -1)
error->one(FLERR,"Restrain atoms {} {} missing on "
"proc {} at step {}", ids[m][0],ids[m][1],
error->one(FLERR,"Restrain atoms {} {} missing on proc {} at step {}", ids[m][0],ids[m][1],
comm->me,update->ntimestep);
}