diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 335fa38d64..c1c97b6e59 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -156,6 +156,7 @@ OPT. * :doc:`precession/spin ` * :doc:`press/berendsen ` * :doc:`print ` + * :doc:`propel/self ` * :doc:`property/atom (k) ` * :doc:`python/invoke ` * :doc:`python/move ` @@ -243,4 +244,3 @@ OPT. * * * - * diff --git a/doc/src/fix.rst b/doc/src/fix.rst index ae0c3d625a..801479de6e 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -308,6 +308,7 @@ accelerated styles exist. * :doc:`precession/spin ` - * :doc:`press/berendsen ` - pressure control by Berendsen barostat * :doc:`print ` - print text and variables during a simulation +* :doc:`propel/self ` - model self-propelled particles * :doc:`property/atom ` - add customized per-atom values * :doc:`python/invoke ` - call a Python function during a simulation * :doc:`python/move ` - call a Python function during a simulation run diff --git a/doc/src/fix_propel_self.rst b/doc/src/fix_propel_self.rst new file mode 100644 index 0000000000..005657c759 --- /dev/null +++ b/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 ` 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) `. + +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) `, :ref:`(Bialke) ` and :ref:`(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 `. + +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 `, :doc:`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 diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index c9f1c0cd8e..43a4207fc1 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/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 diff --git a/examples/USER/misc/propel_self/2d_velocity/in.2d_langevin b/examples/USER/misc/propel_self/2d_velocity/in.2d_langevin new file mode 100644 index 0000000000..6a23e0d005 --- /dev/null +++ b/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 + + diff --git a/examples/USER/misc/propel_self/2d_velocity/in.2d_viscous b/examples/USER/misc/propel_self/2d_velocity/in.2d_viscous new file mode 100644 index 0000000000..1283a5d574 --- /dev/null +++ b/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 diff --git a/examples/USER/misc/propel_self/3d_quaternion/in.3d_quaternion b/examples/USER/misc/propel_self/3d_quaternion/in.3d_quaternion new file mode 100644 index 0000000000..d6b6aff7dd --- /dev/null +++ b/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 diff --git a/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi b/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi index d6c5c566aa..57fc37c10e 100644 --- a/src/MAKE/OPTIONS/Makefile.kokkos_cuda_mpi +++ b/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_PATH = -FFT_LIB = +FFT_INC = -I/usr/include/ +FFT_PATH = -L/usr/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 diff --git a/src/USER-MISC/README b/src/USER-MISC/README index ab646aab5a..9fdc18fbe5 100644 --- a/src/USER-MISC/README +++ b/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 diff --git a/src/USER-MISC/fix_propel_self.cpp b/src/USER-MISC/fix_propel_self.cpp new file mode 100644 index 0000000000..921ce9b9ef --- /dev/null +++ b/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 +#include +#include + +#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 +void FixPropelSelf::post_force_quaternion(int /* vflag */ ) +{ + double **f = atom->f; + AtomVecEllipsoid *av = static_cast(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 +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; +} diff --git a/src/USER-MISC/fix_propel_self.h b/src/USER-MISC/fix_propel_self.h new file mode 100644 index 0000000000..47d1e1255b --- /dev/null +++ b/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 void post_force_velocity(int); + template 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. + +*/