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

reformat with clang-format and related changes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-05-11 19:32:57 -04:00
parent 0d4cca5af9
commit 8a507cc7da
10 changed files with 468 additions and 694 deletions
Download Patch File Download Diff File Expand all files Collapse all files

88
src/USER-BROWNIAN/fix_brownian.cpp
View File

@ -19,79 +19,60 @@
#include "fix_brownian.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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);
}
/* ---------------------------------------------------------------------- */
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 <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D> void FixBrownian::initial_integrate_templated()
{
double **x = atom->x;
double **v = atom->v;
@ -110,33 +90,31 @@ 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));
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 * random->gaussian());
dy = dt * (g1 * f[i][1] + g2 * random->gaussian());
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));
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 * random->gaussian());
dy = dt * (g1 * f[i][1] + g2 * random->gaussian());
dz = dt * (g1 * f[i][2] + g2 * random->gaussian());
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];
@ -145,17 +123,13 @@ void FixBrownian::initial_integrate_templated()
}
x[i][0] += dx;
v[i][0] = dx/dt;
v[i][0] = dx / dt;
x[i][1] += dy;
v[i][1] = dy/dt;
v[i][1] = dy / dt;
x[i][2] += dz;
v[i][2] = dz/dt;
v[i][2] = dz / dt;
}
}
return;

14
src/USER-BROWNIAN/fix_brownian.h
View File

@ -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 <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D>
void initial_integrate_templated();
};
}
} // namespace LAMMPS_NS
#endif
#endif

223
src/USER-BROWNIAN/fix_brownian_asphere.cpp
View File

@ -19,65 +19,44 @@
#include "fix_brownian_asphere.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixBrownianAsphere::FixBrownianAsphere(LAMMPS *lmp, int narg, char **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
@ -88,9 +67,7 @@ void FixBrownianAsphere::init()
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) {
@ -100,67 +77,62 @@ void FixBrownianAsphere::init()
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 );
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 <int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D>
void FixBrownianAsphere::initial_integrate_templated()
{
double **x = atom->x;
@ -180,24 +152,17 @@ void FixBrownianAsphere::initial_integrate_templated()
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];
@ -207,51 +172,42 @@ void FixBrownianAsphere::initial_integrate_templated()
// 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;
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]*random->gaussian()*g2;
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];
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;
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]*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;
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];
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
@ -260,99 +216,56 @@ void FixBrownianAsphere::initial_integrate_templated()
// 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;
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]*random->gaussian()*g2;
tmp[1] = g1*tmp[1]*gamma_t_inv[1]
+ gamma_t_invsqrt[1]*random->gaussian()*g2;
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];
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;
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]*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;
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];
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 );
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];
}
*/
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 );
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;

19
src/USER-BROWNIAN/fix_brownian_asphere.h
View File

@ -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 <int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D>
void initial_integrate_templated();
};
}
} // namespace LAMMPS_NS
#endif
#endif

154
src/USER-BROWNIAN/fix_brownian_base.cpp
View File

@ -19,28 +19,26 @@
#include "fix_brownian.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
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;
@ -49,141 +47,122 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) :
gamma_t_eigen_flag = gamma_r_eigen_flag = 0;
dipole_flag = 0;
if (narg < 5)
error->all(FLERR,"Illegal fix brownian command.");
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.");
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) {
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) {
} else if (strcmp(arg[iarg + 1], "gaussian") == 0) {
noise_flag = 1;
gaussian_noise_flag = 1;
} else if (strcmp(arg[iarg + 1],"none") == 0) {
} 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.");
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);
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.");
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]);
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_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.");
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.");
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);
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.");
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);
rng = new RanMars(lmp, seed + comm->me);
}
/* ---------------------------------------------------------------------- */
@ -201,24 +180,18 @@ 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()
@ -230,21 +203,16 @@ void FixBrownianBase::init()
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;
}

20
src/USER-BROWNIAN/fix_brownian_base.h
View File

@ -11,7 +11,6 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef LMP_FIX_BROWNIAN_BASE_H
#define LMP_FIX_BROWNIAN_BASE_H
@ -19,7 +18,6 @@
namespace LAMMPS_NS {
class FixBrownianBase : public Fix {
public:
FixBrownianBase(class LAMMPS *, int, char **);
@ -31,21 +29,17 @@ class FixBrownianBase : public Fix {
protected:
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
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, 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
@ -53,18 +47,12 @@ class FixBrownianBase : public Fix {
int gaussian_noise_flag; // 0/1 for uniform/gaussian noise
double temp; // temperature
double g1, g2; // prefactors in time stepping
double g1,g2; // prefactors in time stepping
class RanMars *random;
class RanMars *rng;
};
}
} // namespace LAMMPS_NS
#endif
/* ERROR/WARNING messages:

159
src/USER-BROWNIAN/fix_brownian_sphere.cpp
View File

@ -19,101 +19,74 @@
#include "fix_brownian_sphere.h"
#include <cmath>
#include <cstring>
#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 <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
/* ---------------------------------------------------------------------- */
FixBrownianSphere::FixBrownianSphere(LAMMPS *lmp, int narg, char **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 <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D>
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;
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 * random->gaussian());
dy = dt * (g1 * f[i][1] + g2 * random->gaussian());
wz = random->gaussian()*g4;
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;
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 * 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;
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;
v[i][0] = dx / dt;
x[i][1] += dy;
v[i][1] = dy/dt;
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];
v[i][2] = dz / dt;
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;
}

18
src/USER-BROWNIAN/fix_brownian_sphere.h
View File

@ -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 <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D>
void initial_integrate_templated();
double g3,g4;
double g3, g4;
};
}
} // namespace LAMMPS_NS
#endif
#endif

165
src/USER-BROWNIAN/fix_propel_self.cpp
View File

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

14
src/USER-BROWNIAN/fix_propel_self.h
View File

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