diff --git a/doc/src/fix_active.txt b/doc/src/fix_active.txt new file mode 100644 index 0000000000..3bc6c15326 --- /dev/null +++ b/doc/src/fix_active.txt @@ -0,0 +1,50 @@ +"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c + +:link(lws,http://lammps.sandia.gov) +:link(ld,Manual.html) +:link(lc,Section_commands.html#comm) + +:line + +fix active command :pre + +[Syntax:] + +fix ID group-ID magnitude + +ID, group-ID are documented in "fix"_fix.html command :ulb,l +addforce = style name of this fix command :l +magnitude = magnitude of the active force :l +:ule + +[Examples:] + +fix active_group all active 1.0 +fix constant_velocity all viscous 1.0 + +[Description:] + +Adds a force of a constant magnitude to each atom in the group. The direction +of the force is along the axis obtained by rotating the z-axis along the atom's +quaternion. In other words, the force is along the z-axis in the atom's body +frame. + +:line + +[Restart, fix_modify, output, run start/stop, minimize info:] + +No information about this fix is written to "binary restart +files"_restart.html. + +This fix is not imposed during minimization. + +[Restrictions:] + +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. Therefore, this fix only works with +ellipsoidal particles. + +[Related commands:] + +"fix setforce"_fix_setforce.html, "fix addforce"_fix_addforce.html diff --git a/examples/USER/misc/active/active_langevin.in b/examples/USER/misc/active/active_langevin.in new file mode 100644 index 0000000000..3854729fc5 --- /dev/null +++ b/examples/USER/misc/active/active_langevin.in @@ -0,0 +1,61 @@ +# Uses a combnation of fix active and fix langevin to achieve +# self-propelled particles. + +dimension 3 +boundary p p p +units lj + +# Fix active uses quaternions to orient the active force so you need +# ellipsoidal particles, even if you do not care about asymmetric shape +atom_style ellipsoid + + +# Set up box and particles: +variable L equal 10 +region total block -$L $L -$L $L -$L $L +create_box 1 total +variable N equal 1000 +create_atoms 1 random $N 123456 total + +# Set particle shape: +set group all shape 1 1 1 +set group all quat/random 18238 +variable V equal "4.0*PI" +variable rho equal "1.0 / v_V" +# Density so that mass = 1 +set group all density ${rho} + +# Simple repulsive LJ +pair_style lj/cut 1.1225 +pair_coeff * * 1.0 1.0 +pair_modify shift yes + +# Remove initial overlap: +minimize 1e-4 1e-6 1000 10000 +reset_timestep 0 + +# Fixes for time integration +fix step all nve/asphere +fix active all active 1.0 +fix temp all langevin 1.0 1.0 1.0 12341 angmom 3.333333333333 + +variable DUMP_OUT string "active_langevin.dump" +variable MSD_OUT string "active_langevin_msd.dat" + +# Compute temperature and orientation +compute T all temp/asphere +compute quat all property/atom quatw quati quatj quatk +compute msd all msd +compute vacf all vacf + +# Some output: +thermo_style custom time step pe ke etotal temp c_T +thermo 1000 + +dump traj all custom 100 ${DUMP_OUT} id type x y z & + c_quat[1] c_quat[2] c_quat[3] c_quat[4] + +fix msd all ave/time 1 10 50 c_msd[4] c_vacf[1] c_vacf[2] c_vacf[3] c_vacf[4] file ${MSD_OUT} + +timestep 0.001 +run 10000 diff --git a/examples/USER/misc/active/active_viscous.in b/examples/USER/misc/active/active_viscous.in new file mode 100644 index 0000000000..f7660187e1 --- /dev/null +++ b/examples/USER/misc/active/active_viscous.in @@ -0,0 +1,66 @@ +# Uses a combnation of fix active and fix langevin to achieve +# self-propelled particles. + +dimension 3 +boundary p p p +units lj + +# Fix active uses quaternions to orient the active force so you need +# ellipsoidal particles, even if you do not care about asymmetric shape +atom_style ellipsoid + + +# Set up box and particles: +variable L equal 10 +region total block -$L $L -$L $L -$L $L +create_box 1 total +variable N equal 1000 +create_atoms 1 random $N 123456 total + +# Set particle shape: +set group all shape 1 1 1 +set group all quat/random 18238 +variable V equal "4.0*PI" +variable rho equal "1.0 / v_V" +# Density so that mass = 1 +set group all density ${rho} + +# Simple repulsive LJ +pair_style lj/cut 1.1225 +pair_coeff * * 0.0 1.0 +pair_modify shift yes + +# Remove initial overlap: +minimize 1e-4 1e-6 1000 10000 +reset_timestep 0 + +# Fixes for time integration +fix step all nve/asphere + +# Should lead to constant velocity of 10 per particle +fix active all active 1.0 +fix visc all viscous 0.1 + +# Initial velocities: +velocity all create 1.0 1234 + +variable DUMP_OUT string "active_viscous.dump" +variable MSD_OUT string "active_viscous_msd.dat" + +# Compute temperature and orientation +compute T all temp/asphere +compute quat all property/atom quatw quati quatj quatk +compute msd all msd +compute vacf all vacf + +# Some output: +thermo_style custom time step pe ke etotal temp c_T +thermo 1000 + +dump traj all custom 100 ${DUMP_OUT} id type x y z & + c_quat[1] c_quat[2] c_quat[3] c_quat[4] + +fix msd all ave/time 1 10 50 c_msd[4] c_vacf[1] c_vacf[2] c_vacf[3] c_vacf[4] file ${MSD_OUT} + +timestep 0.001 +run 10000 diff --git a/src/KOKKOS/fix_momentum_kokkos.cpp b/src/KOKKOS/fix_momentum_kokkos.cpp index 3b615e82e9..2d4911bfda 100644 --- a/src/KOKKOS/fix_momentum_kokkos.cpp +++ b/src/KOKKOS/fix_momentum_kokkos.cpp @@ -121,13 +121,12 @@ void FixMomentumKokkos::end_of_step() auto xflag2 = xflag; auto yflag2 = yflag; auto zflag2 = zflag; - const Few &vcm_ref = vcm; Kokkos::parallel_for(nlocal, LAMMPS_LAMBDA(int i) { if (mask(i) & groupbit2) { - if (xflag2) v(i,0) -= vcm_ref[0]; - if (yflag2) v(i,1) -= vcm_ref[1]; - if (zflag2) v(i,2) -= vcm_ref[2]; + if (xflag2) v(i,0) -= vcm[0]; + if (yflag2) v(i,1) -= vcm[1]; + if (zflag2) v(i,2) -= vcm[2]; } }); atomKK->modified(execution_space, V_MASK); @@ -161,10 +160,6 @@ void FixMomentumKokkos::end_of_step() auto prd = Few(domain->prd); auto h = Few(domain->h); auto triclinic = domain->triclinic; - - const Few &xcm_ref = xcm; - const Few &omega_ref = omega; - Kokkos::parallel_for(nlocal, LAMMPS_LAMBDA(int i) { if (mask[i] & groupbit2) { Few x_i; @@ -172,12 +167,12 @@ void FixMomentumKokkos::end_of_step() x_i[1] = x(i,1); x_i[2] = x(i,2); auto unwrap = DomainKokkos::unmap(prd,h,triclinic,x_i,image(i)); - auto dx = unwrap[0] - xcm_ref[0]; - auto dy = unwrap[1] - xcm_ref[1]; - auto dz = unwrap[2] - xcm_ref[2]; - v(i,0) -= omega_ref[1]*dz - omega_ref[2]*dy; - v(i,1) -= omega_ref[2]*dx - omega_ref[0]*dz; - v(i,2) -= omega_ref[0]*dy - omega_ref[1]*dx; + auto dx = unwrap[0] - xcm[0]; + auto dy = unwrap[1] - xcm[1]; + auto dz = unwrap[2] - xcm[2]; + v(i,0) -= omega[1]*dz - omega[2]*dy; + v(i,1) -= omega[2]*dx - omega[0]*dz; + v(i,2) -= omega[0]*dy - omega[1]*dx; } }); atomKK->modified(execution_space, V_MASK); @@ -208,3 +203,4 @@ template class FixMomentumKokkos; template class FixMomentumKokkos; #endif } + diff --git a/src/USER-MISC/fix_active.cpp b/src/USER-MISC/fix_active.cpp new file mode 100644 index 0000000000..1bb0726748 --- /dev/null +++ b/src/USER-MISC/fix_active.cpp @@ -0,0 +1,163 @@ +/* ---------------------------------------------------------------------- + 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 +#include + +#include "fix_active.h" + +#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; + +/* + Might be used later, who knows... +static const char* cite_fix_active = + "fix active command:\n\n" + "@article{paquay-2018,\n" + " author = {Paquay, Stefan},\n" + " doi = {},\n" + " issn = {},\n" + " journal = {},\n" + " month = ,\n" + " number = ,\n" + " pages = ,\n" + " title = ,\n" + " volume = ,\n" + " year = \n" + "}\n\n"; +*/ + +static constexpr const bool debug_out = false; + +FixActive::FixActive( LAMMPS *lmp, int narg, char **argv ) + : Fix(lmp, narg, argv) +{ + // if( lmp->citeme) lmp->citeme->add(cite_fix_active); + if( narg < 4 ) error->all(FLERR, "Illegal fix active command"); + + AtomVecEllipsoid *av = static_cast(atom->avec); + if (!av) error->all(FLERR, "FixActive requires atom_style ellipsoid"); + + if( debug_out && comm->me == 0 ){ + fprintf(screen, "This is fix active, narg is %d\n", narg ); + fprintf(screen, "args:"); + for( int i = 0; i < narg; ++i ){ + fprintf(screen, " %s", argv[i]); + } + fprintf(screen, "\n"); + } + + // args: fix ID all active magnitude prop1 prop2 prop3 + // Optional args are + magnitude = force->numeric( FLERR, argv[3] ); +} + + +FixActive::~FixActive() +{} + + +int FixActive::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + + return mask; +} + + +double FixActive::memory_usage() +{ + double bytes = 0.0; + return bytes; +} + + + +void FixActive::post_force(int /* vflag */ ) +{ + // Then do the rest: + double **f = atom->f; + + AtomVecEllipsoid *av = static_cast(atom->avec); + + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + AtomVecEllipsoid::Bonus *bonus = av->bonus; + // Add the active force to the atom force: + for( int i = 0; i < nlocal; ++i ){ + if( mask[i] & groupbit ){ + double f_act[3] = { 0.0, 0.0, 1.0 }; + double f_rot[3]; + + int* ellipsoid = atom->ellipsoid; + double *quat = bonus[ellipsoid[i]].quat; + tagint *tag = atom->tag; + + double Q[3][3]; + MathExtra::quat_to_mat( quat, Q ); + MathExtra::matvec( Q, f_act, f_rot ); + + if (debug_out && comm->me == 0) { + // Magical reference particle: + if (tag[i] == 12) { + fprintf(screen, "rotation quaternion: (%f %f %f %f)\n", + quat[0], quat[1], quat[2], quat[3]); + fprintf(screen, "rotation matrix: / %f %f %f \\\n", + Q[0][0] ,Q[0][1], Q[0][2]); + fprintf(screen, " | %f %f %f |\n", + Q[1][0] ,Q[1][1], Q[1][2]); + fprintf(screen, " \\ %f %f %f /\n", + Q[2][0] ,Q[2][1], Q[2][2]); + + fprintf(screen, "Active force on atom %d: (%f %f %f)\n", + tag[i], f_rot[0], f_rot[1], f_rot[2]); + fprintf(screen, " Total force before: (%f %f %f)\n", + f[i][0], f[i][1], f[i][2]); + } + } + + f[i][0] += magnitude * f_rot[0]; + f[i][1] += magnitude * f_rot[1]; + f[i][2] += magnitude * f_rot[2]; + + } + } +} diff --git a/src/USER-MISC/fix_active.h b/src/USER-MISC/fix_active.h new file mode 100644 index 0000000000..02a91db14a --- /dev/null +++ b/src/USER-MISC/fix_active.h @@ -0,0 +1,55 @@ +/* -*- 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(active,FixActive) + +#else + +#ifndef LMP_FIX_ACTIVE_H +#define LMP_FIX_ACTIVE_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixActive : public Fix { + public: + FixActive(class LAMMPS *, int, char **); + virtual ~FixActive(); + virtual int setmask(); + virtual void post_force(int); + // virtual void post_force_respa(int, int, int); + + double memory_usage(); + +protected: + double magnitude; + int thermostat_orient; +}; + +} + +#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. + +*/