Renamed fix active to fix propel/self

"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

/* ----------------------------------------------------------------------

   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 <stdio.h>

#include <string.h>

#include "fix_propel_self.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;

static constexpr const bool debug_out = false;

FixPropelSelf::FixPropelSelf( 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 propel/self command");

  AtomVecEllipsoid *av = static_cast<AtomVecEllipsoid*>(atom->avec);

  if (!av) error->all(FLERR, "FixPropelSelf 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] );

}

FixPropelSelf::~FixPropelSelf()

{}

int FixPropelSelf::setmask()

{

  int mask = 0;

  mask |= POST_FORCE;

  return mask;

}

double FixPropelSelf::memory_usage()

{

  double bytes = 0.0;

  return bytes;

}

void FixPropelSelf::post_force(int /* vflag */ )

{

  // Then do the rest:

  double **f = atom->f;

  AtomVecEllipsoid *av = static_cast<AtomVecEllipsoid*>(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];

    }

  }

}

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

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

*/