forked from lijiext/lammps
Added a fix to add active force to particles.
This commit is contained in:
Stefan Paquay
Pierre de Buyl
parent
5139f3af33
commit
c93ca5b4a4
|
|
@ -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
|
||||||
|
|
@ -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
|
||||||
|
|
@ -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
|
||||||
|
|
@ -121,13 +121,12 @@ void FixMomentumKokkos<DeviceType>::end_of_step()
|
||||||
auto xflag2 = xflag;
|
auto xflag2 = xflag;
|
||||||
auto yflag2 = yflag;
|
auto yflag2 = yflag;
|
||||||
auto zflag2 = zflag;
|
auto zflag2 = zflag;
|
||||||
const Few<double,3> &vcm_ref = vcm;
|
|
||||||
|
|
||||||
Kokkos::parallel_for(nlocal, LAMMPS_LAMBDA(int i) {
|
Kokkos::parallel_for(nlocal, LAMMPS_LAMBDA(int i) {
|
||||||
if (mask(i) & groupbit2) {
|
if (mask(i) & groupbit2) {
|
||||||
if (xflag2) v(i,0) -= vcm_ref[0];
|
if (xflag2) v(i,0) -= vcm[0];
|
||||||
if (yflag2) v(i,1) -= vcm_ref[1];
|
if (yflag2) v(i,1) -= vcm[1];
|
||||||
if (zflag2) v(i,2) -= vcm_ref[2];
|
if (zflag2) v(i,2) -= vcm[2];
|
||||||
}
|
}
|
||||||
});
|
});
|
||||||
atomKK->modified(execution_space, V_MASK);
|
atomKK->modified(execution_space, V_MASK);
|
||||||
|
|
@ -161,10 +160,6 @@ void FixMomentumKokkos<DeviceType>::end_of_step()
|
||||||
auto prd = Few<double,3>(domain->prd);
|
auto prd = Few<double,3>(domain->prd);
|
||||||
auto h = Few<double,6>(domain->h);
|
auto h = Few<double,6>(domain->h);
|
||||||
auto triclinic = domain->triclinic;
|
auto triclinic = domain->triclinic;
|
||||||
|
|
||||||
const Few<double,3> &xcm_ref = xcm;
|
|
||||||
const Few<double,3> &omega_ref = omega;
|
|
||||||
|
|
||||||
Kokkos::parallel_for(nlocal, LAMMPS_LAMBDA(int i) {
|
Kokkos::parallel_for(nlocal, LAMMPS_LAMBDA(int i) {
|
||||||
if (mask[i] & groupbit2) {
|
if (mask[i] & groupbit2) {
|
||||||
Few<double,3> x_i;
|
Few<double,3> x_i;
|
||||||
|
|
@ -172,12 +167,12 @@ void FixMomentumKokkos<DeviceType>::end_of_step()
|
||||||
x_i[1] = x(i,1);
|
x_i[1] = x(i,1);
|
||||||
x_i[2] = x(i,2);
|
x_i[2] = x(i,2);
|
||||||
auto unwrap = DomainKokkos::unmap(prd,h,triclinic,x_i,image(i));
|
auto unwrap = DomainKokkos::unmap(prd,h,triclinic,x_i,image(i));
|
||||||
auto dx = unwrap[0] - xcm_ref[0];
|
auto dx = unwrap[0] - xcm[0];
|
||||||
auto dy = unwrap[1] - xcm_ref[1];
|
auto dy = unwrap[1] - xcm[1];
|
||||||
auto dz = unwrap[2] - xcm_ref[2];
|
auto dz = unwrap[2] - xcm[2];
|
||||||
v(i,0) -= omega_ref[1]*dz - omega_ref[2]*dy;
|
v(i,0) -= omega[1]*dz - omega[2]*dy;
|
||||||
v(i,1) -= omega_ref[2]*dx - omega_ref[0]*dz;
|
v(i,1) -= omega[2]*dx - omega[0]*dz;
|
||||||
v(i,2) -= omega_ref[0]*dy - omega_ref[1]*dx;
|
v(i,2) -= omega[0]*dy - omega[1]*dx;
|
||||||
}
|
}
|
||||||
});
|
});
|
||||||
atomKK->modified(execution_space, V_MASK);
|
atomKK->modified(execution_space, V_MASK);
|
||||||
|
|
@ -208,3 +203,4 @@ template class FixMomentumKokkos<LMPDeviceType>;
|
||||||
template class FixMomentumKokkos<LMPHostType>;
|
template class FixMomentumKokkos<LMPHostType>;
|
||||||
#endif
|
#endif
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
||||||
|
|
@ -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 <stdio.h>
|
||||||
|
#include <string.h>
|
||||||
|
|
||||||
|
#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<AtomVecEllipsoid*>(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<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];
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
@ -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.
|
||||||
|
|
||||||
|
*/
|
||||||
Loading…
Reference in New Issue