added brownian dynamics integrator fix bd/sphere.

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 <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) <GoldsteinCM>`), :math:`dW_t` and
+:math:`dW_r` are Wiener processes (see e.g. :ref:`(GardinerC) <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) <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) <Dunweg1>` for why this works). This is the same method
+of noise generation as used in :doc:`fix_langevin <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) <Dunweg1>`.
+
+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 <restart>`.
+None of the :doc:`fix_modify <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 <Howto_output>`.
+No parameter of this fix can be used with the *start/stop* keywords of
+the :doc:`run <run>` command.  This fix is not invoked during
+:doc:`energy minimization <minimize>`.
+
+Restrictions
+""""""""""""
+
+This fix requires that atoms store torque and angular velocity (omega)
+as defined by the :doc:`atom_style sphere <atom_style>` command.
+They must also store a dipole moment as defined by the
+:doc:`atom_style dipole <atom_style>` command.
+
+All particles in the group must be finite-size spheres.  They cannot
+be point particles.
+
+Related commands
+""""""""""""""""
+
+:doc:`fix langevin <fix_langevin>`, :doc:`fix nve/sphere <fix_nve_sphere>`,
+:doc:`atom style <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).
+
+

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 <math.h>
+#include <stdio.h>
+#include <string.h>
+#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;
+}

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.
+
+
+*/

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
 ------------------------------------------------------------------------- */

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