diff --git a/doc/src/fix_bd_sphere.rst b/doc/src/fix_bd_sphere.rst new file mode 100644 index 0000000000..2ab184e4fb --- /dev/null +++ b/doc/src/fix_bd_sphere.rst @@ -0,0 +1,168 @@ +.. index:: fix bd/sphere + +fix bd/sphere command +====================== + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group-ID bd/sphere gamma_t gamma_r diff_t diff_r seed keyword args + +* ID, group-ID are documented in :doc:`fix ` command +* bd/sphere = style name of this fix command +* gamma_t = translational friction coefficient +* gamma_r = rotational friction coefficient +* diff_t = translational diffusion coefficient +* diff_r = rotational diffusion coefficient +* zero or more keyword/value pairs may be appended +* keyword = *rng* or *rotate_planar* + + .. parsed-literal:: + + *rng* arg = *uniform* or *gaussian* or *none* + uniform = use uniform random number generator + gaussian = use gaussian random number generator + none = turn off noise + *rotate_planar* arg = *yes* or *no* + yes = only allow rotations in 2D plane (2D simulation only) + no = allow for full 3D rotations + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all bd/sphere 1.0 1.0 1.0 1.0 1294019 + fix 1 all bd/sphere 1.0 1.0 1.0 1.0 19581092 rng none + fix 1 all bd/sphere 1.0 1.0 1.0 1.0 19581092 rotate_planar yes + fix 1 all bd/sphere 1.0 1.0 1.0 1.0 19581092 rotate_planar no rng gaussian + + +Description +""""""""""" + +Perform Brownian Dynamics integration to update position, velocity, +angular velocity, and dipole moment for finite-size spherical particles +in the group each timestep. Brownian Dynamics uses Newton's laws of +motion in the limit that inertial forces are negligible compared to +viscous forces. The stochastic equations of motion are + +.. math:: + + dr = \frac{F}{\gamma_t}dt+\sqrt{2D_t}dW_t, \\ + d\Omega = \frac{T}{\gamma_r}dt + \sqrt{2D_r}dW_r, + +where :math:`d\Omega` is an infinitesimal rotation vector (see e.g. +Chapter 4 of :ref:`(GoldsteinCM) `), :math:`dW_t` and +:math:`dW_r` are Wiener processes (see e.g. :ref:`(GardinerC) `). +The dipole vectors :math:`e_i` are updated using the rotation matrix + +.. math:: + + e_i(t+dt) = e^{\theta_X} e_i(t),\\ + +where :math:`\omega=d\Omega/dt` is the angular velocity, +:math:`\Delta\theta=|\omega|dt` is the rotation angle about +the :math:`\omega` axis, and +:math:`(\theta_X)_{ij}=-\epsilon_{ijk}d\Omega_k` is the +infinitesimal rotation matrix (see e.g. :ref:`(Callegari1) `, +section 7.4). + +IMPORTANT NOTE: This integrator is designed for generic non-equilibrium +simulations with additive noise. There are two important cases which +(conceptually) reduce the number of free parameters in this fix. +(a) In equilibrium simulations +(where fluctuation dissipation theorems are obeyed), one can define +the thermal energy :math:`k_bT=D_t\gamma_t=D_r\gamma_r`. +(b) When a no-slip boundary condition is expected between the spheres and +the surrounding medium, the condition +:math:`\gamma_t=3\gamma_r/\sigma^2` must be explicitly +accounted for (e.g. by setting *gamma_t* to 3 and *gamma_r* to 1) where +:math:`sigma` is the particle diameter. +If both (a) and (b) are true, then one must ensure this explicitly via +the above relationships. + +If the *rng* keyword is used with the *uniform* value, then the noise +is generated from a uniform distribution (see +:ref:`(Dunweg) ` for why this works). This is the same method +of noise generation as used in :doc:`fix_langevin `. + +If the *rng* keyword is used with the *gaussian* value, then the noise +is generated from a gaussian distribution. Typically this added +complexity is unnecessary, and one should be fine using the *uniform* +value for reasons argued in :ref:`(Dunweg) `. + +If the *rng* keyword is used with the *none* value, then the noise +terms are set to zero. + +If the *rotate_planar* keyword is used with the *yes* value, then only +two dimensional rotational diffusion occurs (i.e. only the z-component +of the angular momentum is non-zero). This option is only available to +2D simulations. + +If the *rotate_planar* keyword is used with the *no* value, then three +dimensional rotational diffusion occurs regardless of the simulation +dimension. + +---------- + +.. include:: accel_styles.rst + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files `. +None of the :doc:`fix_modify ` options +are relevant to this fix. No global or per-atom quantities are stored +by this fix for access by various :doc:`output commands `. +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. This fix is not invoked during +:doc:`energy minimization `. + +Restrictions +"""""""""""" + +This fix requires that atoms store torque and angular velocity (omega) +as defined by the :doc:`atom_style sphere ` command. +They must also store a dipole moment as defined by the +:doc:`atom_style dipole ` command. + +All particles in the group must be finite-size spheres. They cannot +be point particles. + +Related commands +"""""""""""""""" + +:doc:`fix langevin `, :doc:`fix nve/sphere `, +:doc:`atom style ` + +Default +""""""" + +The default for *rng* is *uniform* and the default for *rotate_planar* +is *no*. + +---------- + +.. _GoldsteinCM: + +**(GoldsteinCM)** Goldstein, Poole, and Safko, Classical Mechanics, 3rd Ed. (2001). + +.. _GardinerC: + +**(GardinerC)** Gardiner, A Handbook for the Natural and Social Sciences 4th Ed. (2009). + +.. _Callegari1: + +**(Callegari1)** Callegari and Volpe, *Numerical Simulations of Active Brownian +Particles*, Flowing Matter, 211-238 (2019). + +.. _Dunweg1: + +**(Dunweg)** Dunweg and Paul, Int J of Modern Physics C, 2, 817-27 (1991). + + diff --git a/src/fix_bd_sphere.cpp b/src/fix_bd_sphere.cpp new file mode 100644 index 0000000000..a30e8150f6 --- /dev/null +++ b/src/fix_bd_sphere.cpp @@ -0,0 +1,264 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp. + + Contributing author: Sam Cameron (University of Bristol) +------------------------------------------------------------------------- */ + +#include +#include +#include +#include "fix_bd_sphere.h" +#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" + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixBdSphere::FixBdSphere(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + time_integrate = 1; + + if (narg != 8 && narg != 10 && narg != 12) + error->all(FLERR,"Illegal fix bd/sphere command."); + + if (!atom->sphere_flag) + error->all(FLERR,"Fix bd/sphere requires atom style sphere"); + if (!atom->mu_flag) + error->all(FLERR,"Fix bd/sphere requires atom attribute mu"); + + gamma_t = force->numeric(FLERR,arg[3]); + if (gamma_t <= 0.0) + error->all(FLERR,"Fix bd/sphere translational viscous drag " + "coefficient must be > 0."); + + gamma_r = force->numeric(FLERR,arg[4]); + if (gamma_t <= 0.0) + error->all(FLERR,"Fix bd/sphere rotational viscous drag " + "coefficient must be > 0."); + + diff_t = force->numeric(FLERR,arg[5]); + if (diff_t <= 0.0) + error->all(FLERR,"Fix bd/sphere translational diffusion " + "coefficient must be > 0."); + + diff_r = force->numeric(FLERR,arg[6]); + if (diff_r <= 0.0) + error->all(FLERR,"Fix bd/sphere rotational diffusion " + "coefficient must be > 0."); + + seed = force->inumeric(FLERR,arg[7]); + if (seed <= 0) error->all(FLERR,"Fix bd/sphere seed must be > 0."); + + noise_flag = 1; + gaussian_noise_flag = 0; + rotate_planar_flag = 0; + + int iarg == 8; + + while (iarg < narg) { + if (strcmp(arg[iarg],"rng") == 0) { + 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/bd/sphere command."); + } + } else if (strcmp(arg[iarg],"rotate_planar") == 0) { + + if (strcmp(arg[iarg + 1],"yes") == 0) { + rotate_planar_flag = 1; + if (domain->dimension != 2) { + error->all(FLERR,"Cannot constrain rotational degrees of freedom " + "to the xy plane if the simulation is in 3D " + "(in fix/bd/sphere)."); + } + } else if (strcmp(arg[iarg + 1],"no") == 0) { + rotate_planar_flag = 0; + } else { + error->all(FLERR,"Illegal fix/bd/sphere command."); + } + } else { + error->all(FLERR,"Illegal fix/bd/sphere command."); + } + iarg = iarg + 2; + } + + + // initialize Marsaglia RNG with processor-unique seed + random = new RanMars(lmp,seed + comm->me); + + +} + +/* ---------------------------------------------------------------------- */ + +int FixBdSphere::setmask() +{ + int mask = 0; + mask |= INITIAL_INTEGRATE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +FixBdSphere::~FixBdSphere() +{ + + delete random; + +} + + +/* ---------------------------------------------------------------------- */ + +void FixBdSphere::init() +{ + + g1 = force->ftm2v/gamma_t; + g3 = force->ftm2v/gamma_r; + if (noise_flag == 0) { + g2 = 0; + g4 = 0; + rng_func = &RanMars::zero_rng; + } else if (gaussian_noise_flag == 1) { + g2 = sqrt(2 * diff_t); + g4 = sqrt(2 * diff_r); + rng_func = &RanMars::gaussian; + } else { + g2 = sqrt( 24 * diff_t); + g4 = sqrt( 24 * diff_r ); + rng_func = &RanMars::uniform_middle; + } + + if (domain->dimension == 2 && rotate_planar_flag == 0) { + error->warning(FLERR,"Using a 2D simulation, but allowing for " + "full (3D) rotation (in fix/bd/sphere)."); + } + + +} + +/* ---------------------------------------------------------------------- */ + +void FixBdSphere::initial_integrate(int /* vflag */) +{ + double **x = atom->x; + double **v = atom->v; + double **mu = atom->mu; + double **f = atom->f; + double **omega = atom->omega; + int *mask = atom->mask; + int nlocal = atom->nlocal; + double dx,dy,dz; + double dtheta; + double mux,muy,muz,mu_tmp,wx,wy,wz; + double prefac_1, prefac_2; + + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // set timestep here since dt may have changed + dt = update->dt; + sqrtdt = sqrt(dt); + + int d3rot; // whether to compute angular momentum in xy plane + + if (rotate_planar_flag) { + d3rot = 0; + } else { + d3rot = 1; + } + + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + + dx = (dt * g1 * f[i][0] + + sqrtdt * g2 * (random->*rng_func)()); + x[i][0] += dx; + v[i][0] = dx/dt; + dy = (dt * g1 * f[i][1] + + sqrtdt * g2 * (random->*rng_func)()); + x[i][1] += dy; + v[i][1] = dy/dt; + + dz = (dt * g1 * f[i][2] + + sqrtdt * g2 * (random->*rng_func)()); + x[i][2] += dz; + v[i][2] = dz/dt; + + + omega[i][0] = d3rot*(g3* torque[i][0] + + g4 * (random->*rng_func)()/sqrtdt); + omega[i][1] = d3rot*(g3* torque[i][1] + + g4 * (random->*rng_func)()/sqrtdt); + omega[i][2] = (g3* torque[i][2] + + g4 * (random->*rng_func)()/sqrtdt); + + dtheta = sqrt((omega[i][0]*dt)**2+(omega[i][1]*dt)**2+(omega[i][2]*dt)**2); + + if (abs(dtheta) < 1e-14) { + prefac_1 = dt; + prefac_2 = 0.5*dt*dt; + } else { + prefac_1 = dt*sin(dtheta)/dtheta; + prefac_2 = dt*dt*(1-cos(dtheta))/(dtheta*dtheta); + } + + mux = mu[i][0]; + muy = mu[i][1]; + muz = mu[i][2]; + + wx = omega[i][0]; + wy = omega[i][1]; + wz = omega[i][2]; + + mu[i][0] = (mux + prefac_1 * ( -wz*muy + wy*muz ) + + prefac_2 * ( -1*( wz*wz + wy*wy ) * mux + + ( wz*muz + wy*muy ) * wx)); + + mu[i][1] = (muy + prefac_1 * ( wz*mux - wx*muz ) + + prefac_2 * ( -1*(wz*wz + wx*wx) * muy + + ( wz*muz + wx*mux ) * wy)); + + mu[i][2] = (muz + prefac_1 * ( -wy*mux + wx*muy ) + + prefac_2 * ( -1*( wx*wx + wy*wy ) * muz + + ( wy*muy + wx*mux ) * wz)); + + mu_tmp = sqrt(mu[i][0]*mu[i][0]+mu[i][1]*mu[i][1]+mu[i][2]*mu[i][2]); + + mu[i][0] = mu[i][0]/mu_tmp; + mu[i][1] = mu[i][1]/mu_tmp; + mu[i][2] = mu[i][2]/mu_tmp; + + } + } + + return; +} diff --git a/src/fix_bd_sphere.h b/src/fix_bd_sphere.h new file mode 100644 index 0000000000..2777427bc8 --- /dev/null +++ b/src/fix_bd_sphere.h @@ -0,0 +1,104 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + http://lammps.sandia.gov, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef FIX_CLASS + +FixStyle(bd/sphere,FixBdSphere) + +#else + +#ifndef LMP_FIX_BD_SPHERE_H +#define LMP_FIX_BD_SPHERE_H + +#include "fix_nve.h" + +namespace LAMMPS_NS { + +class FixBdSphere : public Fix { + public: + FixBdSphere(class LAMMPS *, int, char **); + virtual ~FixBdSphere(); + void init(); + void initial_integrate(int); + int setmask(); + + private: + typedef double (RanMars::*rng_member)(); + rng_member rng_func; // placeholder for RNG function + int seed; // RNG seed + + double dt, sqrtdt; // time step interval and its sqrt + + + double gamma_t,gamma_r; // translational and rotational damping params + double diff_t,diff_r; // translational and rotational diffusion coeffs + + double g1,g2, g3, g4; // prefactors in time stepping + int noise_flag; // 0/1 for noise off/on + int gaussian_noise_flag; // 0/1 for uniform/gaussian noise + int rotate_planar_flag; // 0/1 for 2D/3D rotational dof + +protected: + class RanMars *random; + + int activity_flag; +}; + +} +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal fix bd/sphere command. + +Wrong number/type of input arguments. + +E: Compute bd/sphere requires atom style sphere + +Self-explanatory. + +E: Compute bd/sphere requires atom attribute mu + +Self-explanatory. + +E: Fix bd/sphere translational viscous drag coefficient must be > 0. + +Self-explanatory. + +E: Fix bd/sphere rotational viscous drag coefficient must be > 0. + +Self-explanatory. + +E: Fix bd/sphere translational diffusion coefficient must be > 0. + +Self-explanatory. + +E: Fix bd/sphere rotational diffusion coefficient must be > 0. + +Self-explanatory. + +E: Fix bd/sphere seed must be > 0. + +E: Cannot constrain rotational degrees of freedom + to the xy plane if the simulation is in 3D (in fix/bd/sphere). + +Self-explanatory. + +W: Using a 2D simulation, but allowing for full (3D) rotation + (in fix/bd/sphere). + +Self-explanatory. + + +*/ diff --git a/src/random_mars.cpp b/src/random_mars.cpp index 621d7d3008..f31a92624a 100644 --- a/src/random_mars.cpp +++ b/src/random_mars.cpp @@ -94,6 +94,22 @@ double RanMars::uniform() return uni; } +/* ---------------------------------------------------------------------- + uniform RN shifted to be symmetric about zero (for fix bd/sphere). +------------------------------------------------------------------------- */ +double RanMars::uniform_middle() +{ + return uniform()-0.5; +} + +/* ---------------------------------------------------------------------- + Return 0 (for fix/bd/sphere). +------------------------------------------------------------------------- */ +double RanMars::zero_rng() +{ + return 0.0; +} + /* ---------------------------------------------------------------------- gaussian RN ------------------------------------------------------------------------- */ diff --git a/src/random_mars.h b/src/random_mars.h index 1bcd16b051..7028ef54d2 100644 --- a/src/random_mars.h +++ b/src/random_mars.h @@ -23,6 +23,8 @@ class RanMars : protected Pointers { RanMars(class LAMMPS *, int); ~RanMars(); double uniform(); + double uniform_middle(); + double zero_rng(); double gaussian(); double gaussian(double mu, double sigma); double rayleigh(double sigma);