Merge pull request #1847 from pdebuyl/fix-active

Add fix propel/self

doc/src/Commands_fix.rst

@@ -156,6 +156,7 @@ OPT.
    * :doc:`precession/spin <fix_precession_spin>`
    * :doc:`press/berendsen <fix_press_berendsen>`
    * :doc:`print <fix_print>`
+   * :doc:`propel/self <fix_propel_self>`
    * :doc:`property/atom (k) <fix_property_atom>`
    * :doc:`python/invoke <fix_python_invoke>`
    * :doc:`python/move <fix_python_move>`
@@ -243,4 +244,3 @@ OPT.
    *
    *
    *
-   *

doc/src/fix.rst

@@ -308,6 +308,7 @@ accelerated styles exist.
 * :doc:`precession/spin <fix_precession_spin>` -
 * :doc:`press/berendsen <fix_press_berendsen>` - pressure control by Berendsen barostat
 * :doc:`print <fix_print>` - print text and variables during a simulation
+* :doc:`propel/self <fix_propel_self>` - model self-propelled particles
 * :doc:`property/atom <fix_property_atom>` - add customized per-atom values
 * :doc:`python/invoke <fix_python_invoke>` - call a Python function during a simulation
 * :doc:`python/move <fix_python_move>` -  call a Python function during a simulation run

doc/src/fix_propel_self.rst

@@ -0,0 +1,113 @@
+.. index:: fix propel/self
+
+fix propel/self command
+=======================
+
+Syntax
+""""""
+
+fix ID group-ID propel/self mode magnitude keyword values ...
+
+* ID, group-ID are documented in :doc:`fix <fix>` command
+* propel/self = style name of this fix command
+* mode = velocity or quat
+* magnitude = magnitude of the active force
+* one or more keyword/value pairs may be appended to args
+* keyword = *types*
+
+  *types* values = one or more atom types
+
+
+
+Examples
+""""""""
+
+
+.. parsed-literal::
+
+   fix active_group all propel/self velocity 1.0
+   fix constant_velocity all viscous 1.0
+
+   fix active_group all propel/self quat 1.0
+
+   fix active all propel/self quat 1.0 types 1 2 4
+
+Description
+"""""""""""
+
+Adds a force of a constant magnitude to each atom in the group. The nature in
+which the force is added depends on the mode.
+
+For *mode* = *velocity*, the active force acts along the velocity vector of
+each atom. This can be interpreted as a velocity-dependent friction,
+such as proposed by :ref:`(Erdmann) <Erdmann>`.
+
+For *mode* = *quat* the force is applied along the axis obtained
+by rotating the x-axis along the atom's quaternion. In other words, the
+force is along the x-axis in the atom's body frame. This mode requires
+all atoms in the group to have a quaternion, so atom\_style should
+either be ellipsoid or body.  In combination with Langevin thermostat
+for translation and rotation in the overdamped regime, the quaternion
+mode corresponds to the active Brownian particle model introduced by
+:ref:`(Henkes) <Henkes>`, :ref:`(Bialke) <Bialke>` and :ref:`(Fily)
+<Fily>`.
+
+By default, this fix is applied to all atoms in the group. You can
+override this behavior by specifying the atom types the fix should work
+on through the *types* keyword.
+
+
+----------
+
+
+**Restart, fix\_modify, output, run start/stop, minimize info:**
+
+No information about this fix is written to :doc:`binary restart files <restart>`.
+
+This fix is not imposed  during minimization.
+
+Restrictions
+""""""""""""
+
+
+In quat mode, this fix makes use of per-atom quaternions to take
+into account the fact that the orientation can rotate and hence the
+direction of the active force can change. The quat mode
+of this fix only works with atom\_style ellipsoid.
+
+Related commands
+""""""""""""""""
+
+:doc:`fix setforce <fix_setforce>`, :doc:`fix addforce <fix_addforce>`
+
+.. _Erdmann:
+
+
+
+**(Erdmann)** U. Erdmann , W. Ebeling, L. Schimansky-Geier, and F. Schweitzer,
+Eur. Phys. J. B 15, 105-113, 2000.
+
+.. _Henkes:
+
+
+
+**(Henkes)** Henkes, S, Fily, Y., and Marchetti, M. C. Phys. Rev. E, 84, 040301(R), 2011.
+
+.. _Bialke:
+
+
+
+**(Bialke)** J. Bialke, T. Speck, and H Loewen, Phys. Rev. Lett. 108, 168301, 2012.
+
+.. _Fily:
+
+
+
+**(Fily)** Y. Fily and M.C. Marchetti, Phys. Rev. Lett. 108, 235702, 2012.
+
+**Default:** types
+
+
+.. _lws: http://lammps.sandia.gov
+.. _ld: Manual.html
+.. _lc: Commands_all.html

doc/utils/sphinx-config/false_positives.txt

@@ -220,6 +220,7 @@ Bext
 Bfrac
 bgq
 Bh
+Bialke
 Biersack
 bigbig
 bigint
@@ -691,6 +692,7 @@ eangle
 eatom
 Eb
 Eba
+Ebeling
 ebond
 ebook
 ebt
@@ -789,6 +791,7 @@ equilibration
 Equilibria
 equilization
 Ercolessi
+Erdmann
 eradius
 erate
 erc
@@ -879,6 +882,7 @@ filename
 Filename
 filenames
 Filenames
+Fily
 fileper
 filesystem
 Fincham
@@ -975,6 +979,7 @@ gcc
 gcmc
 gdot
 GeC
+Geier
 gencode
 georg
 Georg
@@ -1087,6 +1092,7 @@ Heenen
 Hendrik
 Henin
 Henkelman
+Henkes
 henrich
 Henrich
 Herrmann
@@ -1534,6 +1540,7 @@ lmptype
 ln
 localTemp
 localvectors
+Loewen
 logfile
 logfreq
 logicals
@@ -1598,6 +1605,7 @@ manpages
 manybody
 MANYBODY
 Maras
+Marchetti
 Marrink
 Marroquin
 Marsaglia
@@ -2093,6 +2101,7 @@ outfile
 outmost
 outputss
 Ouyang
+overdamped
 overlayed
 Ovito
 oxdna
@@ -2327,6 +2336,7 @@ quadratically
 quadrupolar
 Quant
 quartic
+quat
 quaternion
 quaternions
 quati
@@ -2529,6 +2539,7 @@ scalexz
 scaleyz
 Schaik
 Schilfgarde
+Schimansky
 Schiotz
 Schlitter
 Schmid

examples/USER/misc/propel_self/2d_velocity/in.2d_langevin

@@ -0,0 +1,54 @@
+dimension 2
+boundary p p p
+
+variable L equal 20
+region total block -$L $L -$L $L -0.5 0.5
+lattice hex 0.3
+create_box 2 total
+create_atoms 1 box
+
+# Set random fraction to passive:
+set type 1 type/fraction 2 0.5 1337
+
+# Purely repulsive particles:
+variable rc equal "2^(1.0/6.0)"
+pair_style lj/cut ${rc}
+pair_coeff * * 1.0 1.0
+pair_modify shift yes
+
+mass * 1.0
+
+fix step all nve
+fix temp all langevin 1.0 1.0 1.0 13
+fix twod all enforce2d
+
+neighbor 0.6 bin
+
+dump traj all custom 250 2d_active.dump.bin id type x y z
+
+thermo_style custom time step pe ke etotal temp
+thermo 1000
+run 5000
+
+group one type 1
+group two type 2
+
+compute ke1 one ke
+compute ke2 two ke
+
+thermo_style custom step pe ke etotal temp c_ke1 c_ke2
+
+fix active all propel/self velocity 1.0
+
+# With active force there is more motion so increase bin size:
+neighbor 1.0 bin
+run 10000
+
+# Only make type 1 active:
+fix active all propel/self velocity 1.0 types 1
+
+# With active force there is more motion so increase bin size:
+neighbor 1.0 bin
+run 10000
+
+

examples/USER/misc/propel_self/2d_velocity/in.2d_viscous

@@ -0,0 +1,37 @@
+dimension 2
+boundary p p p
+
+variable L equal 20
+region total block -$L $L -$L $L -0.5 0.5
+lattice hex 0.3
+create_box 2 total
+create_atoms 1 box
+
+# Set random fraction to passive:
+set type 1 type/fraction 2 0.5 1337
+
+# Purely repulsive particles:
+variable rc equal "2^(1.0/6.0)"
+pair_style lj/cut ${rc}
+pair_coeff * * 1.0 1.0
+pair_modify shift yes
+
+mass * 1.0
+
+fix step all nve
+fix twod all enforce2d
+
+neighbor 0.6 bin
+
+dump traj all custom 250 2d_active.dump.bin id type x y z
+
+thermo_style custom step pe ke etotal temp
+thermo 1000
+run 10000
+
+fix active all propel/self velocity 1.0
+fix fric all viscous 1.0
+
+# With active force there is more motion so increase bin size:
+neighbor 1.0 bin
+run 10000

examples/USER/misc/propel_self/3d_quaternion/in.3d_quaternion

@@ -0,0 +1,40 @@
+dimension 3
+boundary p p p
+
+atom_style ellipsoid
+variable L equal 20
+region total block -$L $L -$L $L -$L $L
+lattice sc 0.1
+create_box 2 total
+create_atoms 1 box
+
+# Set random fraction to passive:
+set type 1 type/fraction 2 0.5 1337
+
+# Purely repulsive particles:
+variable rc equal "2^(1.0/6.0)"
+pair_style lj/cut ${rc}
+pair_coeff * * 1.0 1.0
+pair_modify shift yes
+
+# mass * 1.0
+set type * shape 1.0 1.0 1.0
+set type * density 1.9098593171027443
+set type * quat 0 0 1 0
+
+fix step all nve/asphere
+fix temp all langevin 1.0 1.0 1.0 13 angmom 3.333333333
+
+neighbor 0.6 bin
+
+dump traj all custom 100 3d_active.dump.bin id type x y z fx fy fz
+
+thermo_style custom step pe ke etotal temp
+thermo 100
+run 500
+
+fix active all propel/self quat 1.0
+
+# With active force there is more motion so increase bin size:
+neighbor 1.0 bin
+run 500

src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi

@@ -6,7 +6,7 @@ SHELL = /bin/sh
 # compiler/linker settings
 # specify flags and libraries needed for your compiler
 
-KOKKOS_ABSOLUTE_PATH = $(shell cd $(KOKKOS_PATH); pwd)
+KOKKOS_ABSOLUTE_PATH = /home/stefan/projects/lammps-mine/lib/kokkos
 export MPICH_CXX = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper
 export OMPI_CXX = $(KOKKOS_ABSOLUTE_PATH)/bin/nvcc_wrapper
 CC =		mpicxx
@@ -33,7 +33,7 @@ KOKKOS_ARCH = Kepler35
 # LAMMPS ifdef settings
 # see possible settings in Section 2.2 (step 4) of manual
 
-LMP_INC =	-DLAMMPS_GZIP
+LMP_INC =	-DLAMMPS_GZIP -DLAMMPS_PNG -DLAMMPS_JPEG -DLAMMPS_FFMPEG
 
 # MPI library
 # see discussion in Section 2.2 (step 5) of manual
@@ -55,9 +55,9 @@ MPI_LIB =
 # PATH = path for FFT library
 # LIB = name of FFT library
 
-FFT_INC =    	
+FFT_INC =    	-I/usr/include/
-FFT_PATH = 
+FFT_PATH = 	-L/usr/lib/
-FFT_LIB =	
+FFT_LIB =	-lfftw3
 
 # JPEG and/or PNG library
 # see discussion in Section 2.2 (step 7) of manual
@@ -68,7 +68,7 @@ FFT_LIB =
 
 JPG_INC =       
 JPG_PATH = 	
-JPG_LIB =	
+JPG_LIB =	-ljpeg -lpng
 
 # ---------------------------------------------------------------------
 # build rules and dependencies

src/USER-MISC/README

@@ -60,6 +60,7 @@ fix ipi, Michele Ceriotti (EPFL Lausanne), michele.ceriotti at gmail.com, 24 Nov
 fix npt/cauchy, R. E. Miller (Carleton University), F. Pavia and S. Pattamatta, 12 Jan 2020
 fix nvk, Efrem Braun (UC Berkeley), efrem.braun at gmail.com, https://github.com/lammps/lammps/pull/310
 fix pimd, Yuxing Peng (U Chicago), yuxing at uchicago.edu, 24 Nov 2014
+fix propel/self, Stefan Paquay (Brandeis U), stefanpaquay at gmail.com, 20 Jan 2020
 fix rhok, Ulf Pedersen (Roskilde U), ulf at urp.dk, 25 Sep 2017
 fix smd, Axel Kohlmeyer, akohlmey at gmail.com, 19 May 2008
 fix ti/spring, Rodrigo Freitas (Unicamp/Brazil), rodrigohb at gmail.com, 7 Nov 2013

src/USER-MISC/fix_propel_self.cpp

@@ -0,0 +1,268 @@
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+/* -----------------------------------------------------------------------
+   Contributed by Stefan Paquay @ Brandeis University
+
+   Thanks to Liesbeth Janssen @ Eindhoven University for useful discussions!
+----------------------------------------------------------------------- */
+
+#include "fix_propel_self.h"
+
+#include <cstdio>
+#include <cstring>
+#include <string>
+
+#include "atom.h"
+#include "atom_vec_ellipsoid.h"
+#include "citeme.h"
+#include "comm.h"
+#include "error.h"
+#include "force.h"
+#include "group.h"
+#include "math.h"
+#include "math_const.h"
+#include "math_extra.h"
+#include "math_vector.h"
+#include "modify.h"
+#include "random_mars.h"
+#include "respa.h"
+#include "update.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+using namespace MathConst;
+
+#define PRINT_DEBUG_OUTPUT 0
+
+/* ---------------------------------------------------------------------- */
+
+FixPropelSelf::FixPropelSelf( LAMMPS *lmp, int narg, char **argv )
+  : Fix(lmp, narg, argv), magnitude(0.0),
+    mode(VELOCITY), n_types_filter(0), apply_to_type(NULL)
+{
+  if (narg < 5) error->all(FLERR, "Illegal fix propel/self command");
+
+  // The fix is to support the following cases:
+  // 1. Simple atoms, in which case the force points along the velocity
+  // 2. Aspherical particles with an orientation.
+  // The first argument (mode) is used to differentiate between these.
+
+  // args: fix ID all propel/self mode magnitude
+  // Optional args are
+
+  const char *mode_str = argv[3];
+
+  if (strcmp(mode_str, "velocity") == 0) {
+    mode = VELOCITY;
+
+  } else if (strcmp(mode_str, "quat") == 0) {
+
+    // This mode should only be supported if the atom style has
+    // a quaternion (and if all atoms in the group have it)
+
+    if (!atoms_have_quaternion()) {
+      error->all(FLERR, "All fix atoms need to be extended particles");
+    }
+    mode = QUATERNION;
+
+  } else {
+    char msg[2048];
+    sprintf(msg, "Illegal mode \"%s\" for fix propel/self", mode_str);
+    error->all(FLERR, msg);
+  }
+
+  magnitude = force->numeric( FLERR, argv[4] );
+
+  // Handle rest of args:
+
+  int iarg = 5;
+  while (iarg < narg) {
+
+    if (strcmp(argv[iarg],"types") == 0) {
+
+      apply_to_type = new int[atom->ntypes+1];
+      memset(apply_to_type, 0, atom->ntypes * sizeof(int));
+
+      // consume all following numerical arguments as types
+
+      iarg++;
+      int flag=0;
+      while (iarg < narg) {
+        if (isdigit(argv[iarg][0])) {
+          int thistype = force->inumeric(FLERR,argv[iarg]);
+          if ((thistype < 1) || (thistype > atom->ntypes))
+            error->all(FLERR,"Illegal atom type to types keyword");
+          apply_to_type[thistype] = 1;
+          flag = 1;
+          iarg++;
+        } else break;
+      }
+      if (!flag)
+        error->all(FLERR,"'types' keyword requires at least one type");
+      else
+        n_types_filter = 1;
+
+    } else {
+      error->all(FLERR,"Illegal fix propel/self command.");
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixPropelSelf::~FixPropelSelf()
+{
+  delete[] apply_to_type;
+}
+/* ---------------------------------------------------------------------- */
+
+int FixPropelSelf::setmask()
+{
+  int mask = 0;
+  mask |= POST_FORCE;
+
+  return mask;
+}
+
+/* ---------------------------------------------------------------------- */
+
+double FixPropelSelf::memory_usage()
+{
+  // magnitude + thermostat_orient + mode + n_types_filter + apply_to_type
+  double bytes = sizeof(double) + 3*sizeof(int) + sizeof(int*);
+  bytes += sizeof(int)*atom->ntypes*n_types_filter;
+
+  return bytes;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixPropelSelf::post_force(int vflag )
+{
+  switch(mode) {
+  case QUATERNION:
+    if (n_types_filter) post_force_quaternion<1>(vflag);
+    else                post_force_quaternion<0>(vflag);
+    break;
+  case VELOCITY:
+    if (n_types_filter) post_force_velocity<1>(vflag);
+    else                post_force_velocity<0>(vflag);
+    break;
+  default:
+    ;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+template <int filter_by_type>
+void FixPropelSelf::post_force_quaternion(int /* vflag */ )
+{
+  double **f = atom->f;
+  AtomVecEllipsoid *av = static_cast<AtomVecEllipsoid*>(atom->avec);
+
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int *type = atom->type;
+  int* ellipsoid = atom->ellipsoid;
+
+  AtomVecEllipsoid::Bonus *bonus = av->bonus;
+
+  // Add the active force to the atom force:
+
+  for( int i = 0; i < nlocal; ++i ){
+    if( mask[i] & groupbit ){
+      if (filter_by_type && !apply_to_type[type[i]]) {
+        continue;
+      }
+
+      double f_act[3] = { 1.0, 0.0, 0.0 };
+      double f_rot[3];
+
+      double *quat  = bonus[ellipsoid[i]].quat;
+
+      double Q[3][3];
+      MathExtra::quat_to_mat( quat, Q );
+      MathExtra::matvec( Q, f_act, f_rot );
+
+      f[i][0] += magnitude * f_rot[0];
+      f[i][1] += magnitude * f_rot[1];
+      f[i][2] += magnitude * f_rot[2];
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+template <int filter_by_type>
+void FixPropelSelf::post_force_velocity(int /*vflag*/)
+{
+  double **f = atom->f;
+  double **v = atom->v;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  int *type = atom->type;
+
+  // Add the active force to the atom force:
+
+  for(int i = 0; i < nlocal; ++i) {
+    if( mask[i] & groupbit ){
+      if (filter_by_type && !apply_to_type[type[i]]) {
+        continue;
+      }
+
+      const double *vi = v[i];
+      double f_act[3] = { vi[0], vi[1], vi[2] };
+      double nv2 = vi[0]*vi[0] + vi[1]*vi[1] + vi[2]*vi[2];
+      double fnorm = 0.0;
+      const double TOL = 1e-14;
+
+      if (nv2 > TOL) {
+
+        // Without this check you can run into numerical
+        // issues because fnorm will blow up.
+
+        fnorm = magnitude / sqrt(nv2);
+      }
+
+      f[i][0] += fnorm * f_act[0];
+      f[i][1] += fnorm * f_act[1];
+      f[i][2] += fnorm * f_act[2];
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int FixPropelSelf::atoms_have_quaternion()
+{
+  if (!atom->ellipsoid_flag) {
+    error->all(FLERR, "Mode 'quat' requires atom style ellipsoid");
+    return 0;
+  }
+
+  int *mask = atom->mask;
+  int flag=0,flagall=0;
+
+  // Make sure all atoms have ellipsoid data:
+
+  for (int i = 0; i < atom->nlocal; ++i)
+    if (mask[i] & groupbit)
+      if (atom->ellipsoid[i] < 0) ++flag;
+
+  MPI_Allreduce(&flag,&flagall,1,MPI_INT,MPI_SUM,world);
+  if (flagall > 0) return 0;
+
+  return 1;
+}

src/USER-MISC/fix_propel_self.h

@@ -0,0 +1,71 @@
+/* -*- 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(propel/self,FixPropelSelf)
+
+#else
+
+#ifndef LMP_FIX_PROPEL_SELF_H
+#define LMP_FIX_PROPEL_SELF_H
+
+#include "fix.h"
+
+namespace LAMMPS_NS {
+
+class FixPropelSelf : public Fix {
+ public:
+
+  FixPropelSelf(class LAMMPS *, int, char **);
+  virtual ~FixPropelSelf();
+  virtual int setmask();
+  virtual void post_force(int);
+
+  double memory_usage();
+
+ protected:
+  enum operation_modes {
+    VELOCITY = 0,
+    QUATERNION = 1
+  };
+
+private:
+  double magnitude;
+  int mode;
+
+  // If 0, apply fix to everything in group. If > 0, apply only to those
+  // types i for which i <= n_types_filter _and_ apply_to_type[i] == 1:
+  int n_types_filter;
+  int *apply_to_type; //< Specifies, per type, if the fix applies to it or not.
+
+
+  int atoms_have_quaternion();
+
+  template <int filter_by_type> void post_force_velocity(int);
+  template <int filter_by_type> void post_force_quaternion(int);
+};
+}
+
+#endif
+#endif
+
+/* ERROR/WARNING messages:
+
+E: Illegal ... command
+
+Self-explanatory.  Check the input script syntax and compare to the
+documentation for the command.  You can use -echo screen as a
+command-line option when running LAMMPS to see the offending line.
+
+*/