From 8a507cc7da50ff76232f74f81d7e3b2fb4f574c1 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 11 May 2021 19:32:57 -0400 Subject: [PATCH] reformat with clang-format and related changes --- src/USER-BROWNIAN/fix_brownian.cpp | 136 ++++----- src/USER-BROWNIAN/fix_brownian.h | 14 +- src/USER-BROWNIAN/fix_brownian_asphere.cpp | 325 ++++++++------------- src/USER-BROWNIAN/fix_brownian_asphere.h | 21 +- src/USER-BROWNIAN/fix_brownian_base.cpp | 208 ++++++------- src/USER-BROWNIAN/fix_brownian_base.h | 40 +-- src/USER-BROWNIAN/fix_brownian_sphere.cpp | 221 ++++++-------- src/USER-BROWNIAN/fix_brownian_sphere.h | 18 +- src/USER-BROWNIAN/fix_propel_self.cpp | 165 +++++------ src/USER-BROWNIAN/fix_propel_self.h | 14 +- 10 files changed, 468 insertions(+), 694 deletions(-) diff --git a/src/USER-BROWNIAN/fix_brownian.cpp b/src/USER-BROWNIAN/fix_brownian.cpp index 128bd92042..92a969661b 100644 --- a/src/USER-BROWNIAN/fix_brownian.cpp +++ b/src/USER-BROWNIAN/fix_brownian.cpp @@ -12,86 +12,67 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp. + Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp. Contributing author: Sam Cameron (University of Bristol) ------------------------------------------------------------------------- */ #include "fix_brownian.h" -#include -#include -#include "math_extra.h" #include "atom.h" -#include "force.h" -#include "update.h" #include "comm.h" #include "domain.h" -#include "random_mars.h" -#include "memory.h" #include "error.h" +#include "force.h" +#include "math_extra.h" +#include "memory.h" +#include "random_mars.h" +#include "update.h" +#include +#include using namespace LAMMPS_NS; using namespace FixConst; /* ---------------------------------------------------------------------- */ -FixBrownian::FixBrownian(LAMMPS *lmp, int narg, char **arg) : - FixBrownianBase(lmp, narg, arg) +FixBrownian::FixBrownian(LAMMPS *lmp, int narg, char **arg) : FixBrownianBase(lmp, narg, arg) { - if (dipole_flag || gamma_t_eigen_flag || gamma_r_eigen_flag || gamma_r_flag) { - error->all(FLERR,"Illegal fix brownian command."); + error->all(FLERR, "Illegal fix brownian command."); } - - if (!gamma_t_flag) { - error->all(FLERR,"Illegal fix brownian command."); - } - + if (!gamma_t_flag) { error->all(FLERR, "Illegal fix brownian command."); } } - -/* ---------------------------------------------------------------------- */ - -FixBrownian::~FixBrownian() -{ -} - - - /* ---------------------------------------------------------------------- */ void FixBrownian::init() { - FixBrownianBase::init(); - - g1 /= gamma_t; - g2 *= sqrt(gamma_t); - + g1 /= gamma_t; + g2 *= sqrt(gamma_t); } /* ---------------------------------------------------------------------- */ void FixBrownian::initial_integrate(int /*vflag */) { - if (domain->dimension == 2) { if (!noise_flag) { - initial_integrate_templated<0,0,1>(); + initial_integrate_templated<0, 0, 1>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0,1,1>(); + initial_integrate_templated<0, 1, 1>(); } else { - initial_integrate_templated<1,0,1>(); + initial_integrate_templated<1, 0, 1>(); } } else { if (!noise_flag) { - initial_integrate_templated<0,0,0>(); + initial_integrate_templated<0, 0, 0>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0,1,0>(); + initial_integrate_templated<0, 1, 0>(); } else { - initial_integrate_templated<1,0,0>(); + initial_integrate_templated<1, 0, 0>(); } } return; @@ -99,8 +80,7 @@ void FixBrownian::initial_integrate(int /*vflag */) /* ---------------------------------------------------------------------- */ -template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D > -void FixBrownian::initial_integrate_templated() +template void FixBrownian::initial_integrate_templated() { double **x = atom->x; double **v = atom->v; @@ -110,52 +90,46 @@ void FixBrownian::initial_integrate_templated() if (igroup == atom->firstgroup) nlocal = atom->nfirst; - double dx,dy,dz; - + double dx, dy, dz; + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - if (Tp_2D) { - dz = 0; - if (Tp_UNIFORM) { - dx = dt * (g1 * f[i][0] + g2 * (random->uniform()-0.5)); - dy = dt * (g1 * f[i][1] + g2 * (random->uniform()-0.5)); - - } else if (Tp_GAUSS) { - dx = dt * (g1 * f[i][0] + g2 * random->gaussian()); - dy = dt * (g1 * f[i][1] + g2 * random->gaussian()); - } else { - dx = dt * g1 * f[i][0]; - dy = dt * g1 * f[i][1]; - } + dz = 0; + if (Tp_UNIFORM) { + dx = dt * (g1 * f[i][0] + g2 * (rng->uniform() - 0.5)); + dy = dt * (g1 * f[i][1] + g2 * (rng->uniform() - 0.5)); + } else if (Tp_GAUSS) { + dx = dt * (g1 * f[i][0] + g2 * rng->gaussian()); + dy = dt * (g1 * f[i][1] + g2 * rng->gaussian()); + } else { + dx = dt * g1 * f[i][0]; + dy = dt * g1 * f[i][1]; + } } else { - if (Tp_UNIFORM) { - dx = dt * (g1 * f[i][0] + g2 * (random->uniform()-0.5)); - dy = dt * (g1 * f[i][1] + g2 * (random->uniform()-0.5)); - dz = dt * (g1 * f[i][2] + g2 * (random->uniform()-0.5)); - } else if (Tp_GAUSS) { - dx = dt * (g1 * f[i][0] + g2 * random->gaussian()); - dy = dt * (g1 * f[i][1] + g2 * random->gaussian()); - dz = dt * (g1 * f[i][2] + g2 * random->gaussian()); - } else { - dx = dt * g1 * f[i][0]; - dy = dt * g1 * f[i][1]; - dz = dt * g1 * f[i][2]; - } + if (Tp_UNIFORM) { + dx = dt * (g1 * f[i][0] + g2 * (rng->uniform() - 0.5)); + dy = dt * (g1 * f[i][1] + g2 * (rng->uniform() - 0.5)); + dz = dt * (g1 * f[i][2] + g2 * (rng->uniform() - 0.5)); + } else if (Tp_GAUSS) { + dx = dt * (g1 * f[i][0] + g2 * rng->gaussian()); + dy = dt * (g1 * f[i][1] + g2 * rng->gaussian()); + dz = dt * (g1 * f[i][2] + g2 * rng->gaussian()); + } else { + dx = dt * g1 * f[i][0]; + dy = dt * g1 * f[i][1]; + dz = dt * g1 * f[i][2]; + } } - - x[i][0] += dx; - v[i][0] = dx/dt; - - - x[i][1] += dy; - v[i][1] = dy/dt; - - - x[i][2] += dz; - v[i][2] = dz/dt; - - + + x[i][0] += dx; + v[i][0] = dx / dt; + + x[i][1] += dy; + v[i][1] = dy / dt; + + x[i][2] += dz; + v[i][2] = dz / dt; } } return; diff --git a/src/USER-BROWNIAN/fix_brownian.h b/src/USER-BROWNIAN/fix_brownian.h index 7260eb15e3..97705f4f48 100644 --- a/src/USER-BROWNIAN/fix_brownian.h +++ b/src/USER-BROWNIAN/fix_brownian.h @@ -12,9 +12,9 @@ ------------------------------------------------------------------------- */ #ifdef FIX_CLASS - -FixStyle(brownian,FixBrownian) - +// clang-format off +FixStyle(brownian,FixBrownian); +// clang-format on #else #ifndef LMP_FIX_BROWNIAN_H @@ -27,18 +27,16 @@ namespace LAMMPS_NS { class FixBrownian : public FixBrownianBase { public: FixBrownian(class LAMMPS *, int, char **); - virtual ~FixBrownian(); + virtual ~FixBrownian() {}; void init(); void initial_integrate(int); private: - template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D > + template void initial_integrate_templated(); - - }; -} +} // namespace LAMMPS_NS #endif #endif diff --git a/src/USER-BROWNIAN/fix_brownian_asphere.cpp b/src/USER-BROWNIAN/fix_brownian_asphere.cpp index 98c201230a..e7b6eb2728 100644 --- a/src/USER-BROWNIAN/fix_brownian_asphere.cpp +++ b/src/USER-BROWNIAN/fix_brownian_asphere.cpp @@ -12,155 +12,127 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp. + Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp. Contributing author: Sam Cameron (University of Bristol) ------------------------------------------------------------------------- */ #include "fix_brownian_asphere.h" -#include -#include -#include "math_extra.h" #include "atom.h" #include "atom_vec_ellipsoid.h" -#include "force.h" -#include "update.h" #include "comm.h" #include "domain.h" -#include "random_mars.h" -#include "memory.h" #include "error.h" +#include "force.h" +#include "math_extra.h" +#include "memory.h" +#include "random_mars.h" +#include "update.h" +#include +#include using namespace LAMMPS_NS; using namespace FixConst; - /* ---------------------------------------------------------------------- */ FixBrownianAsphere::FixBrownianAsphere(LAMMPS *lmp, int narg, char **arg) : - FixBrownianBase(lmp, narg, arg) + FixBrownianBase(lmp, narg, arg) { - - if (!gamma_t_eigen_flag || !gamma_r_eigen_flag) { - error->all(FLERR,"Illegal fix brownian command."); + error->all(FLERR, "Illegal fix brownian command."); } - if (gamma_t_flag || gamma_r_flag) { - error->all(FLERR,"Illegal fix brownian command."); - } - - + if (gamma_t_flag || gamma_r_flag) error->all(FLERR, "Illegal fix brownian command."); if (dipole_flag && !atom->mu_flag) - error->all(FLERR,"Fix brownian/asphere dipole requires atom attribute mu"); - - + error->all(FLERR, "Fix brownian/asphere dipole requires atom attribute mu"); } -/* ---------------------------------------------------------------------- */ - -FixBrownianAsphere::~FixBrownianAsphere() -{ -} - - - /* ---------------------------------------------------------------------- */ void FixBrownianAsphere::init() { - - avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); - if (!avec) - error->all(FLERR,"Compute brownian/asphere requires " - "atom style ellipsoid"); - + if (!avec) error->all(FLERR, "Compute brownian/asphere requires atom style ellipsoid"); + // check that all particles are finite-size ellipsoids // no point particles allowed, spherical is OK - + int *ellipsoid = atom->ellipsoid; int *mask = atom->mask; int nlocal = atom->nlocal; - + for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) - if (ellipsoid[i] < 0) - error->one(FLERR,"Fix brownian/asphere requires extended particles"); - - + if (ellipsoid[i] < 0) error->one(FLERR, "Fix brownian/asphere requires extended particles"); + if (dipole_flag) { - + double f_rot[3]; double *quat; int *ellipsoid = atom->ellipsoid; AtomVecEllipsoid::Bonus *bonus = avec->bonus; - + double Q[3][3]; - double **mu = atom->mu; - + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - quat = bonus[ellipsoid[i]].quat; - MathExtra::quat_to_mat( quat, Q ); - MathExtra::matvec( Q, dipole_body, f_rot ); - - mu[i][0] = f_rot[0]; - mu[i][1] = f_rot[1]; - mu[i][2] = f_rot[2]; - + quat = bonus[ellipsoid[i]].quat; + MathExtra::quat_to_mat(quat, Q); + MathExtra::matvec(Q, dipole_body, f_rot); + + mu[i][0] = f_rot[0]; + mu[i][1] = f_rot[1]; + mu[i][2] = f_rot[2]; } - } + } } FixBrownianBase::init(); - - } /* ---------------------------------------------------------------------- */ void FixBrownianAsphere::initial_integrate(int /*vflag */) { - if (domain->dimension == 2) { if (dipole_flag) { if (!noise_flag) { - initial_integrate_templated<0,0,1,1>(); + initial_integrate_templated<0, 0, 1, 1>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0,1,1,1>(); + initial_integrate_templated<0, 1, 1, 1>(); } else { - initial_integrate_templated<1,0,1,1>(); + initial_integrate_templated<1, 0, 1, 1>(); } } else { if (!noise_flag) { - initial_integrate_templated<0,0,0,1>(); + initial_integrate_templated<0, 0, 0, 1>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0,1,0,1>(); + initial_integrate_templated<0, 1, 0, 1>(); } else { - initial_integrate_templated<1,0,0,1>(); + initial_integrate_templated<1, 0, 0, 1>(); } } } else { if (dipole_flag) { if (!noise_flag) { - initial_integrate_templated<0,0,1,0>(); + initial_integrate_templated<0, 0, 1, 0>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0,1,1,0>(); + initial_integrate_templated<0, 1, 1, 0>(); } else { - initial_integrate_templated<1,0,1,0>(); + initial_integrate_templated<1, 0, 1, 0>(); } } else { if (!noise_flag) { - initial_integrate_templated<0,0,0,0>(); + initial_integrate_templated<0, 0, 0, 0>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0,1,0,0>(); + initial_integrate_templated<0, 1, 0, 0>(); } else { - initial_integrate_templated<1,0,0,0>(); + initial_integrate_templated<1, 0, 0, 0>(); } } } @@ -169,7 +141,7 @@ void FixBrownianAsphere::initial_integrate(int /*vflag */) /* ---------------------------------------------------------------------- */ -template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D > +template void FixBrownianAsphere::initial_integrate_templated() { double **x = atom->x; @@ -179,180 +151,121 @@ void FixBrownianAsphere::initial_integrate_templated() int nlocal = atom->nlocal; AtomVecEllipsoid::Bonus *bonus = avec->bonus; - - + double **mu = atom->mu; double **torque = atom->torque; - double qw[4]; - double *quat; int *ellipsoid = atom->ellipsoid; - + if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - // project dipole along x axis of quat double f_rot[3]; - double rotationmatrix_transpose[3][3]; - double tmp[3]; double dv[3]; - + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - + // update orientation first - + quat = bonus[ellipsoid[i]].quat; - - MathExtra::quat_to_mat_trans( quat, rotationmatrix_transpose ); - + MathExtra::quat_to_mat_trans(quat, rotationmatrix_transpose); + // tmp holds angular velocity in body frame - MathExtra::matvec(rotationmatrix_transpose,torque[i],tmp); + MathExtra::matvec(rotationmatrix_transpose, torque[i], tmp); if (Tp_2D) { - tmp[0] = tmp[1] = 0.0; - if (Tp_UNIFORM) { - tmp[2] = g1*tmp[2]*gamma_r_inv[2] + gamma_r_invsqrt[2]*(random->uniform()-0.5)*g2; - - } else if (Tp_GAUSS) { - tmp[2] = g1*tmp[2]*gamma_r_inv[2] + gamma_r_invsqrt[2]*random->gaussian()*g2; - - } else { - tmp[2] = g1*tmp[2]*gamma_r_inv[2]; - - } + tmp[0] = tmp[1] = 0.0; + if (Tp_UNIFORM) { + tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g2; + } else if (Tp_GAUSS) { + tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g2; + } else { + tmp[2] = g1 * tmp[2] * gamma_r_inv[2]; + } } else { - if (Tp_UNIFORM) { - - tmp[0] = g1*tmp[0]*gamma_r_inv[0] + gamma_r_invsqrt[0]*(random->uniform()-0.5)*g2; - tmp[1] = g1*tmp[1]*gamma_r_inv[1] + gamma_r_invsqrt[1]*(random->uniform()-0.5)*g2; - tmp[2] = g1*tmp[2]*gamma_r_inv[2] + gamma_r_invsqrt[2]*(random->uniform()-0.5)*g2; - - } else if (Tp_GAUSS) { - - tmp[0] = g1*tmp[0]*gamma_r_inv[0] + gamma_r_invsqrt[0]*random->gaussian()*g2; - tmp[1] = g1*tmp[1]*gamma_r_inv[1] + gamma_r_invsqrt[1]*random->gaussian()*g2; - tmp[2] = g1*tmp[2]*gamma_r_inv[2] + gamma_r_invsqrt[2]*random->gaussian()*g2; - - } else { - tmp[0] = g1*tmp[0]*gamma_r_inv[0]; - tmp[1] = g1*tmp[1]*gamma_r_inv[1]; - tmp[2] = g1*tmp[2]*gamma_r_inv[2]; - - } + if (Tp_UNIFORM) { + tmp[0] = g1 * tmp[0] * gamma_r_inv[0] + gamma_r_invsqrt[0] * (rng->uniform() - 0.5) * g2; + tmp[1] = g1 * tmp[1] * gamma_r_inv[1] + gamma_r_invsqrt[1] * (rng->uniform() - 0.5) * g2; + tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g2; + } else if (Tp_GAUSS) { + tmp[0] = g1 * tmp[0] * gamma_r_inv[0] + gamma_r_invsqrt[0] * rng->gaussian() * g2; + tmp[1] = g1 * tmp[1] * gamma_r_inv[1] + gamma_r_invsqrt[1] * rng->gaussian() * g2; + tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g2; + } else { + tmp[0] = g1 * tmp[0] * gamma_r_inv[0]; + tmp[1] = g1 * tmp[1] * gamma_r_inv[1]; + tmp[2] = g1 * tmp[2] * gamma_r_inv[2]; + } } - + // convert body frame angular velocity to quaternion - MathExtra::quatvec(quat,tmp,qw); - quat[0] = quat[0] + 0.5*dt*qw[0]; - quat[1] = quat[1] + 0.5*dt*qw[1]; - quat[2] = quat[2] + 0.5*dt*qw[2]; - quat[3] = quat[3] + 0.5*dt*qw[3]; - + MathExtra::quatvec(quat, tmp, qw); + quat[0] = quat[0] + 0.5 * dt * qw[0]; + quat[1] = quat[1] + 0.5 * dt * qw[1]; + quat[2] = quat[2] + 0.5 * dt * qw[2]; + quat[3] = quat[3] + 0.5 * dt * qw[3]; + // normalisation introduces the stochastic drift term // to recover the Boltzmann distribution for the case of conservative torques MathExtra::qnormalize(quat); - + // next, update centre of mass positions and velocities - + // tmp now holds force in body frame - MathExtra::matvec(rotationmatrix_transpose,f[i],tmp); + MathExtra::matvec(rotationmatrix_transpose, f[i], tmp); // and then converts to gamma_t^{-1} * F (velocity) in body frame if (Tp_2D) { - tmp[2] = 0.0; - if (Tp_UNIFORM) { - tmp[0] = g1*tmp[0]*gamma_t_inv[0] - + gamma_t_invsqrt[0]*(random->uniform()-0.5)*g2; - tmp[1] = g1*tmp[1]*gamma_t_inv[1] - + gamma_t_invsqrt[1]*(random->uniform()-0.5)*g2; - - } else if (Tp_GAUSS) { - tmp[0] = g1*tmp[0]*gamma_t_inv[0] - + gamma_t_invsqrt[0]*random->gaussian()*g2; - tmp[1] = g1*tmp[1]*gamma_t_inv[1] - + gamma_t_invsqrt[1]*random->gaussian()*g2; - - } else { - tmp[0] = g1*tmp[0]*gamma_t_inv[0]; - tmp[1] = g1*tmp[1]*gamma_t_inv[1]; - - } + tmp[2] = 0.0; + if (Tp_UNIFORM) { + tmp[0] = g1 * tmp[0] * gamma_t_inv[0] + gamma_t_invsqrt[0] * (rng->uniform() - 0.5) * g2; + tmp[1] = g1 * tmp[1] * gamma_t_inv[1] + gamma_t_invsqrt[1] * (rng->uniform() - 0.5) * g2; + } else if (Tp_GAUSS) { + tmp[0] = g1 * tmp[0] * gamma_t_inv[0] + gamma_t_invsqrt[0] * rng->gaussian() * g2; + tmp[1] = g1 * tmp[1] * gamma_t_inv[1] + gamma_t_invsqrt[1] * rng->gaussian() * g2; + } else { + tmp[0] = g1 * tmp[0] * gamma_t_inv[0]; + tmp[1] = g1 * tmp[1] * gamma_t_inv[1]; + } } else { - if (Tp_UNIFORM) { - - tmp[0] = g1*tmp[0]*gamma_t_inv[0] - + gamma_t_invsqrt[0]*(random->uniform()-0.5)*g2; - tmp[1] = g1*tmp[1]*gamma_t_inv[1] - + gamma_t_invsqrt[1]*(random->uniform()-0.5)*g2; - tmp[2] = g1*tmp[2]*gamma_t_inv[2] - + gamma_t_invsqrt[2]*(random->uniform()-0.5)*g2; - - } else if (Tp_GAUSS) { - - tmp[0] = g1*tmp[0]*gamma_t_inv[0] - + gamma_t_invsqrt[0]*random->gaussian()*g2; - tmp[1] = g1*tmp[1]*gamma_t_inv[1] - + gamma_t_invsqrt[1]*random->gaussian()*g2; - tmp[2] = g1*tmp[2]*gamma_t_inv[2] - + gamma_t_invsqrt[2]*random->gaussian()*g2; - - - } else { - - tmp[0] = g1*tmp[0]*gamma_t_inv[0]; - tmp[1] = g1*tmp[1]*gamma_t_inv[1]; - tmp[2] = g1*tmp[2]*gamma_t_inv[2]; - } - } - - - - // finally, convert this back to lab-frame velocity and store in dv - MathExtra::transpose_matvec(rotationmatrix_transpose, tmp, dv ); - - - /* - if (Tp_UNIFORM) { - dv[0] = dt * (g1 * f[i][0] + g2 * (random->uniform()-0.5)); - dv[1] = dt * (g1 * f[i][1] + g2 * (random->uniform()-0.5)); - dv[2] = dt * (g1 * f[i][2] + g2 * (random->uniform()-0.5)); - - } else if (Tp_GAUSS) { - dv[0] = dt * (g1 * f[i][0] + g2 * random->gaussian()); - dv[1] = dt * (g1 * f[i][1] + g2 * random->gaussian()); - dv[2] = dt * (g1 * f[i][2] + g2 * random->gaussian()); - - } else { - dv[0] = dt * g1 * f[i][0]; - dv[1] = dt * g1 * f[i][1]; - dv[2] = dt * g1 * f[i][2]; - + if (Tp_UNIFORM) { + tmp[0] = g1 * tmp[0] * gamma_t_inv[0] + gamma_t_invsqrt[0] * (rng->uniform() - 0.5) * g2; + tmp[1] = g1 * tmp[1] * gamma_t_inv[1] + gamma_t_invsqrt[1] * (rng->uniform() - 0.5) * g2; + tmp[2] = g1 * tmp[2] * gamma_t_inv[2] + gamma_t_invsqrt[2] * (rng->uniform() - 0.5) * g2; + } else if (Tp_GAUSS) { + tmp[0] = g1 * tmp[0] * gamma_t_inv[0] + gamma_t_invsqrt[0] * rng->gaussian() * g2; + tmp[1] = g1 * tmp[1] * gamma_t_inv[1] + gamma_t_invsqrt[1] * rng->gaussian() * g2; + tmp[2] = g1 * tmp[2] * gamma_t_inv[2] + gamma_t_invsqrt[2] * rng->gaussian() * g2; + } else { + tmp[0] = g1 * tmp[0] * gamma_t_inv[0]; + tmp[1] = g1 * tmp[1] * gamma_t_inv[1]; + tmp[2] = g1 * tmp[2] * gamma_t_inv[2]; + } } - */ + + // finally, convert this back to lab-frame velocity and store in dv + MathExtra::transpose_matvec(rotationmatrix_transpose, tmp, dv); + v[i][0] = dv[0]; v[i][1] = dv[1]; v[i][2] = dv[2]; - - x[i][0] += dv[0]*dt; - x[i][1] += dv[1]*dt; - x[i][2] += dv[2]*dt; - + + x[i][0] += dv[0] * dt; + x[i][1] += dv[1] * dt; + x[i][2] += dv[2] * dt; + if (Tp_DIPOLE) { - - MathExtra::quat_to_mat_trans( quat, rotationmatrix_transpose ); - MathExtra::transpose_matvec(rotationmatrix_transpose, dipole_body, f_rot ); - - mu[i][0] = f_rot[0]; - mu[i][1] = f_rot[1]; - mu[i][2] = f_rot[2]; + MathExtra::quat_to_mat_trans(quat, rotationmatrix_transpose); + MathExtra::transpose_matvec(rotationmatrix_transpose, dipole_body, f_rot); + mu[i][0] = f_rot[0]; + mu[i][1] = f_rot[1]; + mu[i][2] = f_rot[2]; } } - } return; diff --git a/src/USER-BROWNIAN/fix_brownian_asphere.h b/src/USER-BROWNIAN/fix_brownian_asphere.h index c229ff1c8b..23a88ddaba 100644 --- a/src/USER-BROWNIAN/fix_brownian_asphere.h +++ b/src/USER-BROWNIAN/fix_brownian_asphere.h @@ -12,9 +12,9 @@ ------------------------------------------------------------------------- */ #ifdef FIX_CLASS - -FixStyle(brownian/asphere,FixBrownianAsphere) - +// clang-format off +FixStyle(brownian/asphere,FixBrownianAsphere); +// clang-format on #else #ifndef LMP_FIX_BROWNIAN_ASPHERE_H @@ -27,25 +27,18 @@ namespace LAMMPS_NS { class FixBrownianAsphere : public FixBrownianBase { public: FixBrownianAsphere(class LAMMPS *, int, char **); - virtual ~FixBrownianAsphere(); + virtual ~FixBrownianAsphere() {}; void initial_integrate(int); - void init(); protected: class AtomVecEllipsoid *avec; - + private: - template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D > + template void initial_integrate_templated(); - - - - - }; - -} +} // namespace LAMMPS_NS #endif #endif diff --git a/src/USER-BROWNIAN/fix_brownian_base.cpp b/src/USER-BROWNIAN/fix_brownian_base.cpp index 1a571c7a74..93990681bd 100644 --- a/src/USER-BROWNIAN/fix_brownian_base.cpp +++ b/src/USER-BROWNIAN/fix_brownian_base.cpp @@ -12,35 +12,33 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp. + Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp. Contributing author: Sam Cameron (University of Bristol) ------------------------------------------------------------------------- */ #include "fix_brownian.h" -#include -#include -#include "math_extra.h" #include "atom.h" -#include "force.h" -#include "update.h" #include "comm.h" #include "domain.h" -#include "random_mars.h" -#include "memory.h" #include "error.h" +#include "force.h" +#include "math_extra.h" +#include "memory.h" +#include "random_mars.h" +#include "update.h" +#include +#include using namespace LAMMPS_NS; using namespace FixConst; /* ---------------------------------------------------------------------- */ -FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) +FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { - time_integrate = 1; noise_flag = 1; @@ -48,142 +46,123 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) : gamma_t_flag = gamma_r_flag = 0; gamma_t_eigen_flag = gamma_r_eigen_flag = 0; dipole_flag = 0; - - if (narg < 5) - error->all(FLERR,"Illegal fix brownian command."); - temp = utils::numeric(FLERR,arg[3],false,lmp); - if (temp <= 0) error->all(FLERR,"Fix brownian temp must be > 0."); - - seed = utils::inumeric(FLERR,arg[4],false,lmp); - if (seed <= 0) error->all(FLERR,"Fix brownian seed must be > 0."); + if (narg < 5) error->all(FLERR, "Illegal fix brownian command."); + temp = utils::numeric(FLERR, arg[3], false, lmp); + if (temp <= 0) error->all(FLERR, "Fix brownian temp must be > 0."); + + seed = utils::inumeric(FLERR, arg[4], false, lmp); + if (seed <= 0) error->all(FLERR, "Fix brownian seed must be > 0."); int iarg = 5; while (iarg < narg) { - if (strcmp(arg[iarg],"rng") == 0) { - if (narg == iarg + 1) { - error->all(FLERR,"Illegal fix brownian command."); - } - if (strcmp(arg[iarg + 1],"uniform") == 0) { - noise_flag = 1; - } else if (strcmp(arg[iarg + 1],"gaussian") == 0) { - noise_flag = 1; - gaussian_noise_flag = 1; - } else if (strcmp(arg[iarg + 1],"none") == 0) { - noise_flag = 0; + if (strcmp(arg[iarg], "rng") == 0) { + if (narg == iarg + 1) error->all(FLERR, "Illegal fix brownian command."); + if (strcmp(arg[iarg + 1], "uniform") == 0) { + noise_flag = 1; + } else if (strcmp(arg[iarg + 1], "gaussian") == 0) { + noise_flag = 1; + gaussian_noise_flag = 1; + } else if (strcmp(arg[iarg + 1], "none") == 0) { + noise_flag = 0; } else { - error->all(FLERR,"Illegal fix brownian command."); + error->all(FLERR, "Illegal fix brownian command."); } iarg = iarg + 2; - } else if (strcmp(arg[iarg],"dipole") == 0) { - if (narg == iarg + 3) { - error->all(FLERR,"Illegal fix brownian command."); - } + } else if (strcmp(arg[iarg], "dipole") == 0) { + if (narg == iarg + 3) error->all(FLERR, "Illegal fix brownian command."); dipole_flag = 1; dipole_body = new double[3]; - - dipole_body[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp); - dipole_body[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp); - dipole_body[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp); + + dipole_body[0] = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + dipole_body[1] = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + dipole_body[2] = utils::numeric(FLERR, arg[iarg + 3], false, lmp); iarg = iarg + 4; - - } else if (strcmp(arg[iarg],"gamma_t_eigen") == 0) { - if (narg == iarg + 3) { - error->all(FLERR,"Illegal fix brownian command."); - } + + } else if (strcmp(arg[iarg], "gamma_t_eigen") == 0) { + if (narg == iarg + 3) error->all(FLERR, "Illegal fix brownian command."); gamma_t_eigen_flag = 1; gamma_t_inv = new double[3]; gamma_t_invsqrt = new double[3]; - gamma_t_inv[0] = 1./utils::numeric(FLERR,arg[iarg+1],false,lmp); - gamma_t_inv[1] = 1./utils::numeric(FLERR,arg[iarg+2],false,lmp); + gamma_t_inv[0] = 1. / utils::numeric(FLERR, arg[iarg + 1], false, lmp); + gamma_t_inv[1] = 1. / utils::numeric(FLERR, arg[iarg + 2], false, lmp); if (domain->dimension == 2) { - if (strcmp(arg[iarg+3],"inf") != 0) { - error->all(FLERR,"Fix brownian gamma_t_eigen third value must be inf for 2D system."); - } - gamma_t_inv[2] = 0; - } else { - gamma_t_inv[2] = 1./utils::numeric(FLERR,arg[iarg+3],false,lmp); + if (strcmp(arg[iarg + 3], "inf") != 0) { + error->all(FLERR, "Fix brownian gamma_t_eigen third value must be inf for 2D system."); + } + gamma_t_inv[2] = 0; + } else { + gamma_t_inv[2] = 1.0 / utils::numeric(FLERR, arg[iarg + 3], false, lmp); } - if (gamma_t_inv[0] < 0 || gamma_t_inv[1] < 0 || gamma_t_inv[2] < 0) - error->all(FLERR,"Fix brownian gamma_t_eigen values must be > 0."); - + if (gamma_t_inv[0] < 0 || gamma_t_inv[1] < 0 || gamma_t_inv[2] < 0) + error->all(FLERR, "Fix brownian gamma_t_eigen values must be > 0."); + gamma_t_invsqrt[0] = sqrt(gamma_t_inv[0]); gamma_t_invsqrt[1] = sqrt(gamma_t_inv[1]); - gamma_t_invsqrt[2] = sqrt(gamma_t_inv[2]); + gamma_t_invsqrt[2] = sqrt(gamma_t_inv[2]); iarg = iarg + 4; - } else if (strcmp(arg[iarg],"gamma_r_eigen") == 0) { - if (narg == iarg + 3) { - error->all(FLERR,"Illegal fix brownian command."); - } + } else if (strcmp(arg[iarg], "gamma_r_eigen") == 0) { + if (narg == iarg + 3) error->all(FLERR, "Illegal fix brownian command."); - gamma_r_eigen_flag = 1; + gamma_r_eigen_flag = 1; gamma_r_inv = new double[3]; gamma_r_invsqrt = new double[3]; - if (domain->dimension == 2) { - if (strcmp(arg[iarg+1],"inf") != 0) { - error->all(FLERR,"Fix brownian gamma_r_eigen first value must be inf for 2D system."); - } - gamma_r_inv[0] = 0; + if (strcmp(arg[iarg + 1], "inf") != 0) { + error->all(FLERR, "Fix brownian gamma_r_eigen first value must be inf for 2D system."); + } + gamma_r_inv[0] = 0; - if (strcmp(arg[iarg+2],"inf") != 0) { - error->all(FLERR,"Fix brownian gamma_r_eigen second value must be inf for 2D system."); - } - gamma_r_inv[1] = 0; + if (strcmp(arg[iarg + 2], "inf") != 0) { + error->all(FLERR, "Fix brownian gamma_r_eigen second value must be inf for 2D system."); + } + gamma_r_inv[1] = 0; } else { - gamma_r_inv[0] = 1./utils::numeric(FLERR,arg[iarg+1],false,lmp); - gamma_r_inv[1] = 1./utils::numeric(FLERR,arg[iarg+2],false,lmp); - + gamma_r_inv[0] = 1. / utils::numeric(FLERR, arg[iarg + 1], false, lmp); + gamma_r_inv[1] = 1. / utils::numeric(FLERR, arg[iarg + 2], false, lmp); } - - gamma_r_inv[2] = 1./utils::numeric(FLERR,arg[iarg+3],false,lmp); - if (gamma_r_inv[0] < 0 || gamma_r_inv[1] < 0 || gamma_r_inv[2] < 0) - error->all(FLERR,"Fix brownian gamma_r_eigen values must be > 0."); - + gamma_r_inv[2] = 1. / utils::numeric(FLERR, arg[iarg + 3], false, lmp); + + if (gamma_r_inv[0] < 0 || gamma_r_inv[1] < 0 || gamma_r_inv[2] < 0) + error->all(FLERR, "Fix brownian gamma_r_eigen values must be > 0."); + gamma_r_invsqrt[0] = sqrt(gamma_r_inv[0]); gamma_r_invsqrt[1] = sqrt(gamma_r_inv[1]); gamma_r_invsqrt[2] = sqrt(gamma_r_inv[2]); iarg = iarg + 4; - - } else if (strcmp(arg[iarg],"gamma_t") == 0) { - if (narg == iarg + 1) { - error->all(FLERR,"Illegal fix brownian command."); - } - + + } else if (strcmp(arg[iarg], "gamma_t") == 0) { + if (narg == iarg + 1) { error->all(FLERR, "Illegal fix brownian command."); } + gamma_t_flag = 1; - gamma_t = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (gamma_t <= 0) - error->all(FLERR,"Fix brownian gamma_t must be > 0."); + gamma_t = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (gamma_t <= 0) error->all(FLERR, "Fix brownian gamma_t must be > 0."); iarg = iarg + 2; - } else if (strcmp(arg[iarg],"gamma_r") == 0) { - if (narg == iarg + 1) { - error->all(FLERR,"Illegal fix brownian command."); - } + } else if (strcmp(arg[iarg], "gamma_r") == 0) { + if (narg == iarg + 1) { error->all(FLERR, "Illegal fix brownian command."); } gamma_r_flag = 1; - gamma_r = utils::numeric(FLERR,arg[iarg+1],false,lmp); - if (gamma_r <= 0) - error->all(FLERR,"Fix brownian gamma_r must be > 0."); + gamma_r = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + if (gamma_r <= 0) error->all(FLERR, "Fix brownian gamma_r must be > 0."); iarg = iarg + 2; } else { - error->all(FLERR,"Illegal fix brownian command."); + error->all(FLERR, "Illegal fix brownian command."); } } - - // initialize Marsaglia RNG with processor-unique seed - random = new RanMars(lmp,seed + comm->me); + // initialize Marsaglia RNG with processor-unique seed + rng = new RanMars(lmp, seed + comm->me); } /* ---------------------------------------------------------------------- */ @@ -201,50 +180,39 @@ FixBrownianBase::~FixBrownianBase() { if (gamma_t_eigen_flag) { - delete [] gamma_t_inv; - delete [] gamma_t_invsqrt; + delete[] gamma_t_inv; + delete[] gamma_t_invsqrt; } if (gamma_r_eigen_flag) { - delete [] gamma_r_inv; - delete [] gamma_r_invsqrt; + delete[] gamma_r_inv; + delete[] gamma_r_invsqrt; } - if (dipole_flag) { - delete [] dipole_body; - } - - - delete random; + if (dipole_flag) { delete[] dipole_body; } + delete rng; } - - /* ---------------------------------------------------------------------- */ void FixBrownianBase::init() { dt = update->dt; sqrtdt = sqrt(dt); - - g1 = force->ftm2v; + + g1 = force->ftm2v; if (noise_flag == 0) { g2 = 0; } else if (gaussian_noise_flag == 1) { - g2 = sqrt(2 * force->boltz*temp/dt/force->mvv2e); + g2 = sqrt(2 * force->boltz * temp / dt / force->mvv2e); } else { - g2 = sqrt( 24 * force->boltz*temp/dt/force->mvv2e); + g2 = sqrt(24 * force->boltz * temp / dt / force->mvv2e); } - - } - - void FixBrownianBase::reset_dt() { double sqrtdt_old = sqrtdt; dt = update->dt; sqrtdt = sqrt(dt); - g2 *= sqrtdt_old/sqrtdt; - + g2 *= sqrtdt_old / sqrtdt; } diff --git a/src/USER-BROWNIAN/fix_brownian_base.h b/src/USER-BROWNIAN/fix_brownian_base.h index 356742f314..5d3b7fe69b 100644 --- a/src/USER-BROWNIAN/fix_brownian_base.h +++ b/src/USER-BROWNIAN/fix_brownian_base.h @@ -11,14 +11,12 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ - #ifndef LMP_FIX_BROWNIAN_BASE_H #define LMP_FIX_BROWNIAN_BASE_H #include "fix.h" namespace LAMMPS_NS { - class FixBrownianBase : public Fix { public: @@ -29,42 +27,32 @@ class FixBrownianBase : public Fix { void reset_dt(); protected: - int seed; // RNG seed - double dt, sqrtdt; // time step interval and its sqrt + int seed; // RNG seed + double dt, sqrtdt; // time step interval and its sqrt + int gamma_t_flag; // 0/1 if isotropic translational damping is unset/set + int gamma_r_flag; // 0/1 if isotropic rotational damping is unset/set + int gamma_t_eigen_flag; // 0/1 if anisotropic translational damping is unset/set + int gamma_r_eigen_flag; // 0/1 if anisotropic rotational damping is unset/set - int gamma_t_flag; // 0/1 if isotropic translational damping is unset/set - int gamma_r_flag; // 0/1 if isotropic rotational damping is unset/set - double gamma_t,gamma_r; // translational and rotational (isotropic) damping params - - - int gamma_t_eigen_flag; // 0/1 if anisotropic translational damping is unset/set - int gamma_r_eigen_flag; // 0/1 if anisotropic rotational damping is unset/set - - double *gamma_t_inv; // anisotropic damping parameter eigenvalues + double gamma_t, gamma_r; // translational and rotational (isotropic) damping params + double *gamma_t_inv; // anisotropic damping parameter eigenvalues double *gamma_r_inv; double *gamma_t_invsqrt; double *gamma_r_invsqrt; - int dipole_flag; // set if dipole is used for asphere double *dipole_body; // direction dipole is slaved to in body frame - int noise_flag; // 0/1 for noise off/on - int gaussian_noise_flag; // 0/1 for uniform/gaussian noise - - double temp; // temperature - - - - double g1,g2; // prefactors in time stepping - - - class RanMars *random; + int noise_flag; // 0/1 for noise off/on + int gaussian_noise_flag; // 0/1 for uniform/gaussian noise + double temp; // temperature + double g1, g2; // prefactors in time stepping + class RanMars *rng; }; -} +} // namespace LAMMPS_NS #endif /* ERROR/WARNING messages: diff --git a/src/USER-BROWNIAN/fix_brownian_sphere.cpp b/src/USER-BROWNIAN/fix_brownian_sphere.cpp index a01eaac719..00893d77e7 100644 --- a/src/USER-BROWNIAN/fix_brownian_sphere.cpp +++ b/src/USER-BROWNIAN/fix_brownian_sphere.cpp @@ -12,108 +12,81 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp. + Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp. Contributing author: Sam Cameron (University of Bristol) ------------------------------------------------------------------------- */ #include "fix_brownian_sphere.h" -#include -#include -#include "math_extra.h" #include "atom.h" -#include "force.h" -#include "update.h" #include "comm.h" #include "domain.h" -#include "random_mars.h" -#include "memory.h" #include "error.h" +#include "force.h" +#include "math_extra.h" +#include "memory.h" +#include "random_mars.h" +#include "update.h" +#include +#include using namespace LAMMPS_NS; using namespace FixConst; - - /* ---------------------------------------------------------------------- */ FixBrownianSphere::FixBrownianSphere(LAMMPS *lmp, int narg, char **arg) : - FixBrownianBase(lmp, narg, arg) + FixBrownianBase(lmp, narg, arg) { - if (gamma_t_eigen_flag || gamma_r_eigen_flag) { - error->all(FLERR,"Illegal fix brownian command."); + error->all(FLERR, "Illegal fix brownian command."); } - if (!gamma_t_flag || !gamma_r_flag) { - error->all(FLERR,"Illegal fix brownian command."); - } - - - if (!atom->mu_flag) - error->all(FLERR,"Fix brownian/sphere requires atom attribute mu"); - - + if (!gamma_t_flag || !gamma_r_flag) { error->all(FLERR, "Illegal fix brownian command."); } + if (!atom->mu_flag) error->all(FLERR, "Fix brownian/sphere requires atom attribute mu"); } - -/* ---------------------------------------------------------------------- */ - -FixBrownianSphere::~FixBrownianSphere() -{ -} - - - /* ---------------------------------------------------------------------- */ void FixBrownianSphere::init() { - FixBrownianBase::init(); - - g3 = g1/gamma_r; - g4 = g2/sqrt(gamma_r); - + + g3 = g1 / gamma_r; + g4 = g2 / sqrt(gamma_r); g1 /= gamma_t; g2 /= sqrt(gamma_t); - - } /* ---------------------------------------------------------------------- */ void FixBrownianSphere::initial_integrate(int /*vflag */) { - if (domain->dimension == 2) { - if (!noise_flag) { - initial_integrate_templated<0,0,1>(); + initial_integrate_templated<0, 0, 1>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0,1,1>(); + initial_integrate_templated<0, 1, 1>(); } else { - initial_integrate_templated<1,0,1>(); + initial_integrate_templated<1, 0, 1>(); } } else { if (!noise_flag) { - initial_integrate_templated<0,0,0>(); + initial_integrate_templated<0, 0, 0>(); } else if (gaussian_noise_flag) { - initial_integrate_templated<0,1,0>(); + initial_integrate_templated<0, 1, 0>(); } else { - initial_integrate_templated<1,0,0>(); + initial_integrate_templated<1, 0, 0>(); } } - return; } - /* ---------------------------------------------------------------------- */ -template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D > +template void FixBrownianSphere::initial_integrate_templated() { double **x = atom->x; @@ -121,112 +94,92 @@ void FixBrownianSphere::initial_integrate_templated() double **f = atom->f; int *mask = atom->mask; int nlocal = atom->nlocal; - double wx,wy,wz; + double wx, wy, wz; double **torque = atom->torque; double **mu = atom->mu; - double mux,muy,muz,mulen; - + double mux, muy, muz, mulen; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; - - - double dx,dy,dz; - + + double dx, dy, dz; + for (int i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { - if (Tp_2D) { - dz = 0; - wx = wy = 0; - if (Tp_UNIFORM) { - dx = dt * (g1 * f[i][0] + g2 * (random->uniform()-0.5)); - dy = dt * (g1 * f[i][1] + g2 * (random->uniform()-0.5)); - wz = (random->uniform()-0.5)*g4; - - } else if (Tp_GAUSS) { - dx = dt * (g1 * f[i][0] + g2 * random->gaussian()); - dy = dt * (g1 * f[i][1] + g2 * random->gaussian()); - wz = random->gaussian()*g4; - - - } else { - dx = dt * g1 * f[i][0]; - dy = dt * g1 * f[i][1]; - wz = 0; - - } + dz = 0; + wx = wy = 0; + if (Tp_UNIFORM) { + dx = dt * (g1 * f[i][0] + g2 * (rng->uniform() - 0.5)); + dy = dt * (g1 * f[i][1] + g2 * (rng->uniform() - 0.5)); + wz = (rng->uniform() - 0.5) * g4; + } else if (Tp_GAUSS) { + dx = dt * (g1 * f[i][0] + g2 * rng->gaussian()); + dy = dt * (g1 * f[i][1] + g2 * rng->gaussian()); + wz = rng->gaussian() * g4; + } else { + dx = dt * g1 * f[i][0]; + dy = dt * g1 * f[i][1]; + wz = 0; + } } else { - if (Tp_UNIFORM) { - dx = dt * (g1 * f[i][0] + g2 * (random->uniform()-0.5)); - dy = dt * (g1 * f[i][1] + g2 * (random->uniform()-0.5)); - dz = dt * (g1 * f[i][2] + g2 * (random->uniform()-0.5)); - wx = (random->uniform()-0.5)*g4; - wy = (random->uniform()-0.5)*g4; - wz = (random->uniform()-0.5)*g4; - - } else if (Tp_GAUSS) { - dx = dt * (g1 * f[i][0] + g2 * random->gaussian()); - dy = dt * (g1 * f[i][1] + g2 * random->gaussian()); - dz = dt * (g1 * f[i][2] + g2 * random->gaussian()); - wx = random->gaussian()*g4; - wy = random->gaussian()*g4; - wz = random->gaussian()*g4; - - - } else { - dx = dt * g1 * f[i][0]; - dy = dt * g1 * f[i][1]; - dz = dt * g1 * f[i][2]; - wx = wy = wz = 0; - - } + if (Tp_UNIFORM) { + dx = dt * (g1 * f[i][0] + g2 * (rng->uniform() - 0.5)); + dy = dt * (g1 * f[i][1] + g2 * (rng->uniform() - 0.5)); + dz = dt * (g1 * f[i][2] + g2 * (rng->uniform() - 0.5)); + wx = (rng->uniform() - 0.5) * g4; + wy = (rng->uniform() - 0.5) * g4; + wz = (rng->uniform() - 0.5) * g4; + } else if (Tp_GAUSS) { + dx = dt * (g1 * f[i][0] + g2 * rng->gaussian()); + dy = dt * (g1 * f[i][1] + g2 * rng->gaussian()); + dz = dt * (g1 * f[i][2] + g2 * rng->gaussian()); + wx = rng->gaussian() * g4; + wy = rng->gaussian() * g4; + wz = rng->gaussian() * g4; + } else { + dx = dt * g1 * f[i][0]; + dy = dt * g1 * f[i][1]; + dz = dt * g1 * f[i][2]; + wx = wy = wz = 0; + } } - - x[i][0] += dx; - v[i][0] = dx/dt; - - - x[i][1] += dy; - v[i][1] = dy/dt; - - - x[i][2] += dz; - v[i][2] = dz/dt; + x[i][0] += dx; + v[i][0] = dx / dt; + + x[i][1] += dy; + v[i][1] = dy / dt; + + x[i][2] += dz; + v[i][2] = dz / dt; + + wx += g3 * torque[i][0]; + wy += g3 * torque[i][1]; + wz += g3 * torque[i][2]; - - wx += g3*torque[i][0]; - wy += g3*torque[i][1]; - wz += g3*torque[i][2]; - - // store length of dipole as we need to convert it to a unit vector and // then back again - - mulen = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]); + + mulen = sqrt(mu[i][0] * mu[i][0] + mu[i][1] * mu[i][1] + mu[i][2] * mu[i][2]); // unit vector at time t - mux = mu[i][0]/mulen; - muy = mu[i][1]/mulen; - muz = mu[i][2]/mulen; - - + mux = mu[i][0] / mulen; + muy = mu[i][1] / mulen; + muz = mu[i][2] / mulen; // un-normalised unit vector at time t + dt - mu[i][0] = mux + (wy*muz - wz*muy)*dt; - mu[i][1] = muy + (wz*mux - wx*muz)*dt; - mu[i][2] = muz + (wx*muy - wy*mux)*dt; - + mu[i][0] = mux + (wy * muz - wz * muy) * dt; + mu[i][1] = muy + (wz * mux - wx * muz) * dt; + mu[i][2] = muz + (wx * muy - wy * mux) * dt; + // normalisation introduces the stochastic drift term due to changing from // Stratonovich to Ito interpretation MathExtra::norm3(mu[i]); // multiply by original magnitude to obtain dipole of same length - mu[i][0] = mu[i][0]*mulen; - mu[i][1] = mu[i][1]*mulen; - mu[i][2] = mu[i][2]*mulen; - + mu[i][0] = mu[i][0] * mulen; + mu[i][1] = mu[i][1] * mulen; + mu[i][2] = mu[i][2] * mulen; } } - - return; } diff --git a/src/USER-BROWNIAN/fix_brownian_sphere.h b/src/USER-BROWNIAN/fix_brownian_sphere.h index e518d28db3..19e86204ad 100644 --- a/src/USER-BROWNIAN/fix_brownian_sphere.h +++ b/src/USER-BROWNIAN/fix_brownian_sphere.h @@ -12,9 +12,9 @@ ------------------------------------------------------------------------- */ #ifdef FIX_CLASS - -FixStyle(brownian/sphere,FixBrownianSphere) - +// clang-format off +FixStyle(brownian/sphere,FixBrownianSphere); +// clang-format on #else #ifndef LMP_FIX_BROWNIAN_SPHERE_H @@ -27,18 +27,16 @@ namespace LAMMPS_NS { class FixBrownianSphere : public FixBrownianBase { public: FixBrownianSphere(class LAMMPS *, int, char **); - virtual ~FixBrownianSphere(); + virtual ~FixBrownianSphere() {}; void init(); - void initial_integrate(int /*vflag */); - + void initial_integrate(int); private: - template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D > + template void initial_integrate_templated(); - double g3,g4; + double g3, g4; }; - -} +} // namespace LAMMPS_NS #endif #endif diff --git a/src/USER-BROWNIAN/fix_propel_self.cpp b/src/USER-BROWNIAN/fix_propel_self.cpp index 067fec3bf8..9adb49d039 100644 --- a/src/USER-BROWNIAN/fix_propel_self.cpp +++ b/src/USER-BROWNIAN/fix_propel_self.cpp @@ -36,60 +36,54 @@ using namespace LAMMPS_NS; using namespace FixConst; -enum{DIPOLE,VELOCITY,QUAT}; +enum { DIPOLE, VELOCITY, QUAT }; #define TOL 1e-14 /* ---------------------------------------------------------------------- */ -FixPropelSelf::FixPropelSelf(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg) +FixPropelSelf::FixPropelSelf(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { virial_global_flag = virial_peratom_flag = 1; - if (narg != 5 && narg != 9) - error->all(FLERR,"Illegal fix propel/self command"); + if (narg != 5 && narg != 9) error->all(FLERR, "Illegal fix propel/self command"); - - if (strcmp(arg[3],"velocity") == 0) { + if (strcmp(arg[3], "velocity") == 0) { mode = VELOCITY; thermo_virial = 0; - } else if (strcmp(arg[3],"dipole") == 0) { + } else if (strcmp(arg[3], "dipole") == 0) { mode = DIPOLE; thermo_virial = 1; - } else if (strcmp(arg[3],"quat") == 0) { + } else if (strcmp(arg[3], "quat") == 0) { mode = QUAT; thermo_virial = 1; } else { - error->all(FLERR,"Illegal fix propel/self command"); + error->all(FLERR, "Illegal fix propel/self command"); } - magnitude = utils::numeric(FLERR,arg[4],false,lmp); + magnitude = utils::numeric(FLERR, arg[4], false, lmp); // check for keyword if (narg == 9) { - if (mode != QUAT) { - error->all(FLERR,"Illegal fix propel/self command"); - } - if (strcmp(arg[5],"qvector") == 0) { - sx = utils::numeric(FLERR,arg[6],false,lmp); - sy = utils::numeric(FLERR,arg[7],false,lmp); - sz = utils::numeric(FLERR,arg[8],false,lmp); - double snorm = sqrt(sx*sx + sy*sy + sz*sz); - sx = sx/snorm; - sy = sy/snorm; - sz = sz/snorm; + if (mode != QUAT) { error->all(FLERR, "Illegal fix propel/self command"); } + if (strcmp(arg[5], "qvector") == 0) { + sx = utils::numeric(FLERR, arg[6], false, lmp); + sy = utils::numeric(FLERR, arg[7], false, lmp); + sz = utils::numeric(FLERR, arg[8], false, lmp); + double snorm = sqrt(sx * sx + sy * sy + sz * sz); + sx = sx / snorm; + sy = sy / snorm; + sz = sz / snorm; } else { - error->all(FLERR,"Illegal fix propel/self command"); + error->all(FLERR, "Illegal fix propel/self command"); } } else { sx = 1.0; sy = 0.0; sz = 0.0; } - } /* ---------------------------------------------------------------------- */ @@ -101,27 +95,16 @@ int FixPropelSelf::setmask() return mask; } -/* ---------------------------------------------------------------------- */ - -FixPropelSelf::~FixPropelSelf() -{ - -} - - /* ---------------------------------------------------------------------- */ void FixPropelSelf::init() { if (mode == DIPOLE && !atom->mu_flag) - error->all(FLERR,"Fix propel/self requires atom attribute mu " - "with option dipole."); + error->all(FLERR, "Fix propel/self requires atom attribute mu with option dipole"); if (mode == QUAT) { avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); - if (!avec) - error->all(FLERR,"Fix propel/self requires " - "atom style ellipsoid with option quat."); + if (!avec) error->all(FLERR, "Fix propel/self requires atom style ellipsoid with option quat"); // check that all particles are finite-size ellipsoids // no point particles allowed, spherical is OK @@ -133,17 +116,19 @@ void FixPropelSelf::init() for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) if (ellipsoid[i] < 0) - error->one(FLERR,"Fix propel/self requires extended particles " - "with option quat."); + error->one(FLERR, "Fix propel/self requires extended particles with option quat"); } - } +/* ---------------------------------------------------------------------- */ + void FixPropelSelf::setup(int vflag) { post_force(vflag); } +/* ---------------------------------------------------------------------- */ + void FixPropelSelf::post_force(int vflag) { if (mode == DIPOLE) @@ -152,9 +137,9 @@ void FixPropelSelf::post_force(int vflag) post_force_velocity(vflag); else if (mode == QUAT) post_force_quaternion(vflag); - } +/* ---------------------------------------------------------------------- */ void FixPropelSelf::post_force_dipole(int vflag) { @@ -163,13 +148,14 @@ void FixPropelSelf::post_force_dipole(int vflag) int nlocal = atom->nlocal; double **x = atom->x; double **mu = atom->mu; - double fx,fy,fz; - + double fx, fy, fz; // energy and virial setup double vi[6]; - if (vflag) v_setup(vflag); - else evflag = 0; + if (vflag) + v_setup(vflag); + else + evflag = 0; // if domain has PBC, need to unwrap for virial double unwrap[3]; @@ -179,26 +165,27 @@ void FixPropelSelf::post_force_dipole(int vflag) for (int i = 0; i < nlocal; i++) if (mask[i] & groupbit) { - fx = magnitude*mu[i][0]; - fy = magnitude*mu[i][1]; - fz = magnitude*mu[i][2]; + fx = magnitude * mu[i][0]; + fy = magnitude * mu[i][1]; + fz = magnitude * mu[i][2]; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; if (evflag) { - domain->unmap(x[i],image[i],unwrap); - vi[0] = fx*unwrap[0]; - vi[1] = fy*unwrap[1]; - vi[2] = fz*unwrap[2]; - vi[3] = fx*unwrap[1]; - vi[4] = fx*unwrap[2]; - vi[5] = fy*unwrap[2]; + domain->unmap(x[i], image[i], unwrap); + vi[0] = fx * unwrap[0]; + vi[1] = fy * unwrap[1]; + vi[2] = fz * unwrap[2]; + vi[3] = fx * unwrap[1]; + vi[4] = fx * unwrap[2]; + vi[5] = fy * unwrap[2]; v_tally(i, vi); } } } +/* ---------------------------------------------------------------------- */ void FixPropelSelf::post_force_velocity(int vflag) { @@ -207,23 +194,24 @@ void FixPropelSelf::post_force_velocity(int vflag) double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; - double nv2,fnorm,fx,fy,fz; - + double nv2, fnorm, fx, fy, fz; // energy and virial setup double vi[6]; - if (vflag) v_setup(vflag); - else evflag = 0; + if (vflag) + v_setup(vflag); + else + evflag = 0; // if domain has PBC, need to unwrap for virial double unwrap[3]; imageint *image = atom->image; // Add the active force to the atom force: - for(int i = 0; i < nlocal; ++i) { - if(mask[i] & groupbit) { + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { - nv2 = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]; + nv2 = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]; fnorm = 0.0; if (nv2 > TOL) { @@ -242,69 +230,72 @@ void FixPropelSelf::post_force_velocity(int vflag) f[i][2] += fz; if (evflag) { - domain->unmap(x[i],image[i],unwrap); - vi[0] = fx*unwrap[0]; - vi[1] = fy*unwrap[1]; - vi[2] = fz*unwrap[2]; - vi[3] = fx*unwrap[1]; - vi[4] = fx*unwrap[2]; - vi[5] = fy*unwrap[2]; + domain->unmap(x[i], image[i], unwrap); + vi[0] = fx * unwrap[0]; + vi[1] = fy * unwrap[1]; + vi[2] = fz * unwrap[2]; + vi[3] = fx * unwrap[1]; + vi[4] = fx * unwrap[2]; + vi[5] = fy * unwrap[2]; v_tally(i, vi); } } } } +/* ---------------------------------------------------------------------- */ + void FixPropelSelf::post_force_quaternion(int vflag) { double **f = atom->f; double **x = atom->x; int *mask = atom->mask; int nlocal = atom->nlocal; - int* ellipsoid = atom->ellipsoid; + int *ellipsoid = atom->ellipsoid; // ellipsoidal properties AtomVecEllipsoid::Bonus *bonus = avec->bonus; - double f_act[3] = { sx, sy, sz }; + double f_act[3] = {sx, sy, sz}; double f_rot[3]; double *quat; double Q[3][3]; - double fx,fy,fz; - + double fx, fy, fz; // energy and virial setup double vi[6]; - if (vflag) v_setup(vflag); - else evflag = 0; + if (vflag) + v_setup(vflag); + else + evflag = 0; // if domain has PBC, need to unwrap for virial double unwrap[3]; imageint *image = atom->image; // Add the active force to the atom force: - for(int i = 0; i < nlocal; ++i) { - if(mask[i] & groupbit) { + for (int i = 0; i < nlocal; ++i) { + if (mask[i] & groupbit) { quat = bonus[ellipsoid[i]].quat; MathExtra::quat_to_mat(quat, Q); MathExtra::matvec(Q, f_act, f_rot); - fx = magnitude*f_rot[0]; - fy = magnitude*f_rot[1]; - fz = magnitude*f_rot[2]; + fx = magnitude * f_rot[0]; + fy = magnitude * f_rot[1]; + fz = magnitude * f_rot[2]; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; if (evflag) { - domain->unmap(x[i],image[i],unwrap); - vi[0] = fx*unwrap[0]; - vi[1] = fy*unwrap[1]; - vi[2] = fz*unwrap[2]; - vi[3] = fx*unwrap[1]; - vi[4] = fx*unwrap[2]; - vi[5] = fy*unwrap[2]; + domain->unmap(x[i], image[i], unwrap); + vi[0] = fx * unwrap[0]; + vi[1] = fy * unwrap[1]; + vi[2] = fz * unwrap[2]; + vi[3] = fx * unwrap[1]; + vi[4] = fx * unwrap[2]; + vi[5] = fy * unwrap[2]; v_tally(i, vi); } } diff --git a/src/USER-BROWNIAN/fix_propel_self.h b/src/USER-BROWNIAN/fix_propel_self.h index 154cef0e2d..9874635818 100644 --- a/src/USER-BROWNIAN/fix_propel_self.h +++ b/src/USER-BROWNIAN/fix_propel_self.h @@ -12,9 +12,9 @@ ------------------------------------------------------------------------- */ #ifdef FIX_CLASS - -FixStyle(propel/self,FixPropelSelf) - +// clang-format off +FixStyle(propel/self,FixPropelSelf); +// clang-format on #else #ifndef LMP_FIX_PROPEL_SELF_H @@ -26,7 +26,7 @@ namespace LAMMPS_NS { class FixPropelSelf : public Fix { public: FixPropelSelf(class LAMMPS *, int, char **); - virtual ~FixPropelSelf(); + virtual ~FixPropelSelf() {}; void init(); void post_force(int); void setup(int); @@ -34,7 +34,7 @@ class FixPropelSelf : public Fix { private: double magnitude; - double sx,sy,sz; + double sx, sy, sz; int mode; void post_force_dipole(int); @@ -42,10 +42,8 @@ class FixPropelSelf : public Fix { void post_force_quaternion(int); class AtomVecEllipsoid *avec; - }; - -} +} // namespace LAMMPS_NS #endif #endif