Made changes to propel/self. Kept all features of previous version (and tested they stayed the same), but additionally added dipole option for direction of self-propulsion. Also updated examples.
This commit is contained in:
Sam Cameron
@ -8,53 +8,117 @@ Syntax
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID propel/self mode magnitude keyword values ...
|
||||
fix ID group-ID propel/self 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*
|
||||
* magnitude = magnitude of self-propulsion force
|
||||
* one (and only one) keyword/value pair must be appended to args
|
||||
* keyword = *dipole* or *velocity* or *quat*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*dipole* value = none = apply force along dipole direction
|
||||
*velocity* value = none = apply force along velocity direction
|
||||
*quat* values = direction vector *sx* and *sy* and *sz*
|
||||
*sx* = x component of force direction in ellipsoid frame
|
||||
*sy* = y component of force direction in ellipsoid frame
|
||||
*sz* = z component of force direction in ellipsoid frame
|
||||
|
||||
|
||||
*types* values = one or more atom types
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
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
|
||||
fix propel/self all 40.0 dipole
|
||||
fix propel/self all 10.0 velocity
|
||||
fix propel/self all 15.7 quat 1.0 0.0 0.0
|
||||
|
||||
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.
|
||||
Add a force to each atom in the group due to a self-propulsion force. The
|
||||
force is given by
|
||||
|
||||
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>`.
|
||||
.. math::
|
||||
|
||||
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>`.
|
||||
F_i = f_P e_i
|
||||
|
||||
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.
|
||||
where *i* is the particle the force is being applied to, :math:`f_P`
|
||||
is the magnitude of the force, and :math:`e_i` is the vector direction
|
||||
of the force. The specification of :math:`e_i` is based on which of the
|
||||
three keywords (*dipole* or *velocity* or *quat*) one selects.
|
||||
|
||||
For keyword *dipole*, :math:`e_i` is just equal to
|
||||
the dipole vectors of the atoms in the group. Therefore, if the dipoles
|
||||
are not unit vectors, the :math:`e_i` will not be unit vectors.
|
||||
|
||||
.. note::
|
||||
|
||||
If another command changes the magnitude of the dipole, this force will
|
||||
change accordingly (since :math:`|e_i|` will change, which is physically
|
||||
equivalent to re-scaling :math:`f_P` while keeping :math:`|e_i|` constant),
|
||||
and no warning will be provided by LAMMPS. This is almost never what you
|
||||
want, so ensure you aren't changing dipole magnitudes with another LAMMPS
|
||||
fix or pair style. Furthermore, self-propulsion forces (almost) always
|
||||
set :math:`e_i` to be a unit vector for all times, so it's best to set
|
||||
all the dipole magnitudes to 1.0 unless you have a good reason not to
|
||||
(see the :doc:`set <set>` command on how to do this).
|
||||
|
||||
For keyword *velocity*, :math:`e_i` points in the direction
|
||||
of the current velocity (a unit-vector). This can be interpreted as a
|
||||
velocity-dependent friction, as proposed by e.g. :ref:`(Erdmann) <Erdmann1>`.
|
||||
|
||||
For keyword *quat*, :math:`e_i` points in the direction of the unit
|
||||
vector defined by its arguments *sx*, *sy*, and *sz*, which are
|
||||
themselves defined within the coordinate frame of the atom's
|
||||
ellipsoid. For instance, for an ellipsoid with long axis along
|
||||
its x-direction, if one wanted the self-propulsion force to also
|
||||
be along this axis, set *sx* equal to 1 and *sy*, *sz* both equal
|
||||
to zero. For *quat*, :math:`e_i` will always be a unit vector,
|
||||
so multiplying all three arguments *sx*, *sy*, and *sz* by a
|
||||
positive scalar will not change the self-propulsion force
|
||||
(multiplying by a negative scalar will change the sign of the
|
||||
force).
|
||||
|
||||
Along with adding a force contribution, this fix can also
|
||||
contribute to the virial (pressure) of the system, defined as
|
||||
:math:`f_P \sum_i <e_i . r_i>/(d V)`, where :math:`r_i` is the
|
||||
*unwrapped* coordinate of particle i in the case of periodic
|
||||
boundary conditions. See :ref:`(Winkler) <Winkler1>` for a
|
||||
discussion of this active pressure contribution.
|
||||
|
||||
For keywords *dipole* and *quat*, this fix is by default
|
||||
included in pressure computations.
|
||||
|
||||
For keyword *velocity*, this fix is by default not included
|
||||
in pressure computations.
|
||||
|
||||
|
||||
.. note::
|
||||
|
||||
In contrast to equilibrium systems, pressure of active systems
|
||||
in general depends on the geometry of the container.
|
||||
The active pressure contribution as calculated in this fix
|
||||
is only valid for certain boundary conditions (spherical
|
||||
walls, rectangular walls, or periodic boundary conditions).
|
||||
For other geometries, the pressure must be measured via
|
||||
explicit calculation of the force per unit area on a wall,
|
||||
and so one must not calculate it using this fix.
|
||||
(Use :doc:`fix_modify <fix_modify>` as described below
|
||||
to turn off the virial contribution of this fix). Again,
|
||||
see :ref:`(Winkler) <Winkler1>` for discussion of why this
|
||||
is the case.
|
||||
|
||||
Furthermore, when dealing with active systems, the temperature
|
||||
is no longer well defined. Therefore, one should ensure that
|
||||
the *virial* flag is used in the
|
||||
:doc:`compute pressure <compute_pressure>` command (turning
|
||||
off temperature contributions).
|
||||
|
||||
|
||||
----------
|
||||
|
||||
Restart, fix_modify, output, run start/stop, minimize info
|
||||
@ -62,40 +126,48 @@ 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.
|
||||
The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
|
||||
fix to add the contribution due to the added forces on atoms to the
|
||||
system's virial as part of :doc:`thermodynamic output <thermo_style>`.
|
||||
The default is *virial yes* for keywords *dipole* and *quat*. The
|
||||
default is *virial no* for keyword *velocity*.
|
||||
|
||||
|
||||
No parameter of this fix can be used with the *start/stop* keywords of
|
||||
the :doc:`run <run>` command.
|
||||
|
||||
|
||||
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.
|
||||
With keyword *dipole*, this fix only works when the DIPOLE package is enabled.
|
||||
See the :doc:`Build package <Build_package>` doc page for more info.
|
||||
|
||||
This fix is part of the USER-BROWNIAN package. It is only enabled if
|
||||
LAMMPS was built with that package. See the :doc:`Build package <Build_package>`
|
||||
doc page for more info.
|
||||
|
||||
|
||||
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.
|
||||
:doc:`fix efield <fix_efield>` , :doc:`fix setforce <fix_setforce>`,
|
||||
:doc:`fix addforce <fix_addforce>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
types
|
||||
none
|
||||
|
||||
----------
|
||||
|
||||
|
||||
.. _Erdmann1:
|
||||
|
||||
**(Erdmann)** U. Erdmann , W. Ebeling, L. Schimansky-Geier, and F. Schweitzer,
|
||||
Eur. Phys. J. B 15, 105-113, 2000.
|
||||
|
||||
|
||||
.. _Winkler1:
|
||||
|
||||
**(Winkler)** Winkler, Wysocki, and Gompper, Soft Matter, 11, 6680 (2015).
|
||||
|
||||
70
examples/USER/brownian/spherical_ABP/in.2d_abp
Normal file
70
examples/USER/brownian/spherical_ABP/in.2d_abp
Normal file
@ -0,0 +1,70 @@
|
||||
# 2D overdamped active brownian particle dynamics (ABP)
|
||||
# with WCA potential
|
||||
|
||||
variable gamma_t equal 1.0
|
||||
variable gamma_r equal 1.0
|
||||
variable D_t equal 1.0
|
||||
variable D_r equal 3.0
|
||||
variable seed equal 1974019
|
||||
variable fp equal 4.0
|
||||
|
||||
variable params string ${gamma_t}_${gamma_r}_${D_t}_${D_r}_${fp}
|
||||
|
||||
|
||||
log log_WCA_${params}_2d.lammps.log
|
||||
units lj
|
||||
atom_style hybrid dipole sphere
|
||||
dimension 2
|
||||
newton off
|
||||
|
||||
|
||||
lattice sq 0.4
|
||||
region box block -16 16 -16 16 -0.2 0.2
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
mass * 1.0
|
||||
set type * dipole/random ${seed} 1.0
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
# more careful with neighbors since higher diffusion in abps
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
# WCA potential (purely repulsive)
|
||||
pair_style lj/cut 2.5
|
||||
pair_coeff * * 1.0 1.0 1.1224
|
||||
pair_modify shift yes
|
||||
|
||||
|
||||
|
||||
# overdamped brownian dynamics time-step
|
||||
fix step all brownian/sphere ${gamma_t} ${gamma_r} ${D_t} ${D_r} ${seed} dipole
|
||||
# self-propulsion force along the dipole direction
|
||||
fix activity all propel/self ${fp} dipole
|
||||
fix 2 all enforce2d
|
||||
|
||||
|
||||
compute press all pressure NULL virial
|
||||
|
||||
thermo_style custom step temp epair c_press
|
||||
|
||||
#equilibration
|
||||
timestep 0.0000000001
|
||||
thermo 50001
|
||||
run 50000
|
||||
reset_timestep 0
|
||||
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*] c_press
|
||||
|
||||
|
||||
timestep 0.00001
|
||||
thermo 10000
|
||||
|
||||
# main run
|
||||
run 200000
|
||||
74
examples/USER/brownian/spherical_ABP/in.3d_ideal_abp
Normal file
74
examples/USER/brownian/spherical_ABP/in.3d_ideal_abp
Normal file
@ -0,0 +1,74 @@
|
||||
# 3D overdamped active brownian dynamics with no interactions
|
||||
|
||||
variable gamma_t equal 1.0
|
||||
variable gamma_r equal 1.0
|
||||
variable D_t equal 1.0
|
||||
variable D_r equal 3.0
|
||||
variable seed equal 1974019
|
||||
variable fp equal 4.0
|
||||
|
||||
variable params string ${gamma_t}_${gamma_r}_${D_t}_${D_r}_${fp}
|
||||
|
||||
|
||||
log log_ideal_${params}_3d.lammps.log
|
||||
units lj
|
||||
atom_style hybrid dipole sphere
|
||||
dimension 3
|
||||
newton off
|
||||
|
||||
|
||||
lattice sc 0.4
|
||||
region box block -8 8 -8 8 -8 8
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
mass * 1.0
|
||||
set type * dipole/random ${seed} 1.0
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
pair_style none
|
||||
|
||||
# overdamped brownian dynamics time-step
|
||||
fix step all brownian/sphere ${gamma_t} ${gamma_r} ${D_t} ${D_r} ${seed} dipole
|
||||
# self-propulsion force along the dipole direction
|
||||
fix activity all propel/self ${fp} dipole
|
||||
|
||||
|
||||
compute press all pressure NULL virial
|
||||
|
||||
thermo_style custom step temp epair c_press
|
||||
|
||||
#equilibration
|
||||
timestep 0.0000000001
|
||||
thermo 50001
|
||||
run 50000
|
||||
reset_timestep 0
|
||||
|
||||
|
||||
# MSD to demonstrate expected diffusive behaviour for ideal active
|
||||
# brownian motion, which is
|
||||
#
|
||||
# MSD = (2*d*D_t + 2*fp**2/(gamma_t**2*(d-1)*D_r))*t
|
||||
#
|
||||
# with d being simulation dimension
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*] c_press
|
||||
|
||||
|
||||
timestep 0.00001
|
||||
thermo 10000
|
||||
|
||||
# main run
|
||||
run 120000
|
||||
|
||||
# if you want to check that rotational diffusion is behaving as expected,
|
||||
# uncomment next three lines for dump output and then plot <e(t).e(0)>,
|
||||
# which should decay exponentially with timescale (d-1)*D_r (with d
|
||||
# being simulation dimension)
|
||||
|
||||
#dump 1 all custom 2000 dump_ideal_${params}_3d.lammpstrj id type &
|
||||
# x y xu yu mux muy muz fx fy fz
|
||||
#dump_modify 1 first yes sort id
|
||||
|
||||
#run 120000
|
||||
@ -0,0 +1,163 @@
|
||||
units lj
|
||||
atom_style hybrid dipole sphere
|
||||
WARNING: Atom_style hybrid defines both pertype and peratom masses - both must be set, only peratom masses will be used (src/atom_vec_hybrid.cpp:157)
|
||||
dimension 2
|
||||
newton off
|
||||
|
||||
|
||||
lattice sq 0.4
|
||||
Lattice spacing in x,y,z = 1.5811388 1.5811388 1.5811388
|
||||
region box block -16 16 -16 16 -0.2 0.2
|
||||
create_box 1 box
|
||||
Created orthogonal box = (-25.298221 -25.298221 -0.31622777) to (25.298221 25.298221 0.31622777)
|
||||
2 by 2 by 1 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 1024 atoms
|
||||
create_atoms CPU = 0.001 seconds
|
||||
mass * 1.0
|
||||
set type * dipole/random ${seed} 1.0
|
||||
set type * dipole/random 1974019 1.0
|
||||
Setting atom values ...
|
||||
1024 settings made for dipole/random
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
# more careful with neighbors since higher diffusion in abps
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
# WCA potential (purely repulsive)
|
||||
pair_style lj/cut 2.5
|
||||
pair_coeff * * 1.0 1.0 1.1224
|
||||
pair_modify shift yes
|
||||
|
||||
|
||||
|
||||
# overdamped brownian dynamics time-step
|
||||
fix step all brownian/sphere ${gamma_t} ${gamma_r} ${D_t} ${D_r} ${seed} dipole
|
||||
fix step all brownian/sphere 1 ${gamma_r} ${D_t} ${D_r} ${seed} dipole
|
||||
fix step all brownian/sphere 1 1 ${D_t} ${D_r} ${seed} dipole
|
||||
fix step all brownian/sphere 1 1 1 ${D_r} ${seed} dipole
|
||||
fix step all brownian/sphere 1 1 1 3 ${seed} dipole
|
||||
fix step all brownian/sphere 1 1 1 3 1974019 dipole
|
||||
# self-propulsion force along the dipole direction
|
||||
fix activity all propel/self ${fp} dipole
|
||||
fix activity all propel/self 4 dipole
|
||||
fix 2 all enforce2d
|
||||
|
||||
|
||||
compute press all pressure NULL virial
|
||||
|
||||
thermo_style custom step temp epair c_press
|
||||
|
||||
#equilibration
|
||||
timestep 0.0000000001
|
||||
thermo 50001
|
||||
run 50000
|
||||
Neighbor list info ...
|
||||
update every 1 steps, delay 1 steps, check yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 2.1224
|
||||
ghost atom cutoff = 2.1224
|
||||
binsize = 1.0612, bins = 48 48 1
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair lj/cut, perpetual
|
||||
attributes: half, newton off
|
||||
pair build: half/bin/newtoff
|
||||
stencil: half/bin/2d/newtoff
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 5.052 | 5.052 | 5.052 Mbytes
|
||||
Step Temp E_pair c_press
|
||||
0 1 0 -0.53979198
|
||||
50000 1.0371295e+10 0 -0.542818
|
||||
Loop time of 2.25396 on 4 procs for 50000 steps with 1024 atoms
|
||||
|
||||
Performance: 0.192 tau/day, 22183.200 timesteps/s
|
||||
99.8% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.1034 | 0.10382 | 0.10438 | 0.1 | 4.61
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.26698 | 0.26833 | 0.26924 | 0.2 | 11.90
|
||||
Output | 2.2284e-05 | 2.4926e-05 | 3.2332e-05 | 0.0 | 0.00
|
||||
Modify | 1.7222 | 1.7237 | 1.727 | 0.1 | 76.48
|
||||
Other | | 0.1581 | | | 7.01
|
||||
|
||||
Nlocal: 256.000 ave 256 max 256 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 105.000 ave 105 max 105 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 544.000 ave 544 max 544 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 2176
|
||||
Ave neighs/atom = 2.1250000
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
reset_timestep 0
|
||||
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*] c_press
|
||||
|
||||
|
||||
timestep 0.00001
|
||||
thermo 10000
|
||||
|
||||
# main run
|
||||
run 200000
|
||||
Per MPI rank memory allocation (min/avg/max) = 5.427 | 5.427 | 5.427 Mbytes
|
||||
Step Temp E_pair c_msd[1] c_msd[2] c_msd[3] c_msd[4] c_press
|
||||
0 1.0371295e+10 0 0 0 0 0 -0.542818
|
||||
10000 107356.09 0.079828495 0.19584264 0.19679822 0 0.39264086 0.00078740793
|
||||
20000 101692.44 0.11317533 0.40847364 0.42097802 0 0.82945167 0.74111888
|
||||
30000 105763.55 0.10261852 0.68303669 0.66125751 0 1.3442942 0.71112533
|
||||
40000 105127.29 0.12371743 0.97990144 0.94005552 0 1.919957 1.0574552
|
||||
50000 101579.58 0.12771813 1.3059069 1.2364468 0 2.5423537 1.059263
|
||||
60000 104914.36 0.12055843 1.6215593 1.525488 0 3.1470473 0.79873537
|
||||
70000 106629.18 0.1278745 1.9639958 1.8682794 0 3.8322752 0.91950208
|
||||
80000 103286.54 0.13927689 2.3201565 2.2373383 0 4.5574948 1.1875034
|
||||
90000 106451.61 0.093479681 2.6287902 2.5753347 0 5.2041249 1.0861163
|
||||
100000 102199.72 0.13269425 2.9127976 2.9369237 0 5.8497214 1.4841998
|
||||
110000 105229.73 0.10594209 3.1798718 3.3495317 0 6.5294035 1.5444784
|
||||
120000 106262.36 0.11902575 3.6267452 3.7188125 0 7.3455578 1.3366518
|
||||
130000 109388.12 0.10562576 3.929973 4.0226942 0 7.9526672 1.324534
|
||||
140000 107697.35 0.13028752 4.231893 4.3780995 0 8.6099925 1.7406167
|
||||
150000 103928.72 0.12278994 4.5826286 4.7578662 0 9.3404948 1.3024003
|
||||
160000 103370.23 0.11391216 4.8767011 5.1181189 0 9.99482 1.4325241
|
||||
170000 103821.53 0.11163256 5.153318 5.3785963 0 10.531914 1.4569115
|
||||
180000 106824.99 0.13225083 5.4080929 5.7399804 0 11.148073 1.334984
|
||||
190000 101794.6 0.10632938 5.7384925 6.080955 0 11.819448 0.81285422
|
||||
200000 102128.67 0.13703498 6.0414673 6.4511058 0 12.492573 0.42904128
|
||||
Loop time of 9.60419 on 4 procs for 200000 steps with 1024 atoms
|
||||
|
||||
Performance: 17992.140 tau/day, 20824.236 timesteps/s
|
||||
100.0% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.71756 | 0.75121 | 0.79008 | 3.0 | 7.82
|
||||
Neigh | 0.018158 | 0.018773 | 0.019357 | 0.3 | 0.20
|
||||
Comm | 1.0469 | 1.0597 | 1.0738 | 1.2 | 11.03
|
||||
Output | 0.00051435 | 0.00057078 | 0.00070838 | 0.0 | 0.01
|
||||
Modify | 6.8012 | 6.9883 | 7.1513 | 4.9 | 72.76
|
||||
Other | | 0.7857 | | | 8.18
|
||||
|
||||
Nlocal: 256.000 ave 265 max 240 min
|
||||
Histogram: 1 0 0 0 0 0 1 0 0 2
|
||||
Nghost: 88.5000 ave 91 max 87 min
|
||||
Histogram: 1 0 2 0 0 0 0 0 0 1
|
||||
Neighs: 678.500 ave 713 max 597 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 1 2
|
||||
|
||||
Total # of neighbors = 2714
|
||||
Ave neighs/atom = 2.6503906
|
||||
Neighbor list builds = 241
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:11
|
||||
@ -0,0 +1,150 @@
|
||||
units lj
|
||||
atom_style hybrid dipole sphere
|
||||
WARNING: Atom_style hybrid defines both pertype and peratom masses - both must be set, only peratom masses will be used (src/atom_vec_hybrid.cpp:157)
|
||||
dimension 3
|
||||
newton off
|
||||
|
||||
|
||||
lattice sc 0.4
|
||||
Lattice spacing in x,y,z = 1.3572088 1.3572088 1.3572088
|
||||
region box block -8 8 -8 8 -8 8
|
||||
create_box 1 box
|
||||
Created orthogonal box = (-10.857670 -10.857670 -10.857670) to (10.857670 10.857670 10.857670)
|
||||
2 by 1 by 2 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 4096 atoms
|
||||
create_atoms CPU = 0.002 seconds
|
||||
mass * 1.0
|
||||
set type * dipole/random ${seed} 1.0
|
||||
set type * dipole/random 1974019 1.0
|
||||
Setting atom values ...
|
||||
4096 settings made for dipole/random
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
pair_style none
|
||||
|
||||
# overdamped brownian dynamics time-step
|
||||
fix step all brownian/sphere ${gamma_t} ${gamma_r} ${D_t} ${D_r} ${seed} dipole
|
||||
fix step all brownian/sphere 1 ${gamma_r} ${D_t} ${D_r} ${seed} dipole
|
||||
fix step all brownian/sphere 1 1 ${D_t} ${D_r} ${seed} dipole
|
||||
fix step all brownian/sphere 1 1 1 ${D_r} ${seed} dipole
|
||||
fix step all brownian/sphere 1 1 1 3 ${seed} dipole
|
||||
fix step all brownian/sphere 1 1 1 3 1974019 dipole
|
||||
# self-propulsion force along the dipole direction
|
||||
fix activity all propel/self ${fp} dipole
|
||||
fix activity all propel/self 4 dipole
|
||||
|
||||
|
||||
compute press all pressure NULL virial
|
||||
|
||||
thermo_style custom step temp epair c_press
|
||||
|
||||
#equilibration
|
||||
timestep 0.0000000001
|
||||
thermo 50001
|
||||
run 50000
|
||||
WARNING: No pairwise cutoff or binsize set. Atom sorting therefore disabled. (src/atom.cpp:2118)
|
||||
WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (src/comm_brick.cpp:167)
|
||||
Per MPI rank memory allocation (min/avg/max) = 4.319 | 4.319 | 4.319 Mbytes
|
||||
Step Temp E_pair c_press
|
||||
0 1 0 0.068021726
|
||||
50000 1.046425e+10 0 0.067505633
|
||||
Loop time of 7.83903 on 4 procs for 50000 steps with 4096 atoms
|
||||
|
||||
Performance: 0.055 tau/day, 6378.340 timesteps/s
|
||||
97.9% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.50238 | 0.51891 | 0.53617 | 1.7 | 6.62
|
||||
Output | 2.6343e-05 | 3.6075e-05 | 4.6997e-05 | 0.0 | 0.00
|
||||
Modify | 6.9536 | 6.9732 | 7.0105 | 0.8 | 88.95
|
||||
Other | | 0.3469 | | | 4.43
|
||||
|
||||
Nlocal: 1024.00 ave 1024 max 1024 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 353.000 ave 353 max 353 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 0
|
||||
Ave neighs/atom = 0.0000000
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
reset_timestep 0
|
||||
|
||||
|
||||
# MSD to demonstrate expected diffusive behaviour for ideal active
|
||||
# brownian motion, which is
|
||||
#
|
||||
# MSD = (2*d*D_t + 2*fp**2/(gamma_t**2*(d-1)*D_r))*t
|
||||
#
|
||||
# with d being simulation dimension
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*] c_press
|
||||
|
||||
|
||||
timestep 0.00001
|
||||
thermo 10000
|
||||
|
||||
# main run
|
||||
run 120000
|
||||
WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms may get lost. (src/comm_brick.cpp:167)
|
||||
Per MPI rank memory allocation (min/avg/max) = 4.694 | 4.694 | 4.694 Mbytes
|
||||
Step Temp E_pair c_msd[1] c_msd[2] c_msd[3] c_msd[4] c_press
|
||||
0 1.046425e+10 0 0 0 0 0 0.067505633
|
||||
10000 106340.56 0 0.2469318 0.23662295 0.2441413 0.72769605 0.19939966
|
||||
20000 104620.27 0 0.55429324 0.5231436 0.54976641 1.6272032 0.26601423
|
||||
30000 106130.45 0 0.87562668 0.84813496 0.89321299 2.6169746 0.30836996
|
||||
40000 105773.31 0 1.2262635 1.1899278 1.2626926 3.6788838 0.35392219
|
||||
50000 103804.88 0 1.5851624 1.5645815 1.6434185 4.7931624 0.33326997
|
||||
60000 105746.45 0 1.9928016 1.9347072 1.9837329 5.9112417 0.2550878
|
||||
70000 104500.3 0 2.3269429 2.2774077 2.3368821 6.9412326 0.25218225
|
||||
80000 105381.46 0 2.7114959 2.6937299 2.7171132 8.122339 0.36940892
|
||||
90000 104542.79 0 3.0828648 3.084417 3.0783207 9.2456025 0.36106481
|
||||
100000 104246.75 0 3.4635513 3.5105066 3.4545226 10.42858 0.3712313
|
||||
110000 103099.55 0 3.8471061 3.9389997 3.8220676 11.608173 0.38466185
|
||||
120000 103098.45 0 4.2014598 4.3456831 4.1888659 12.736009 0.36710217
|
||||
Loop time of 22.8893 on 4 procs for 120000 steps with 4096 atoms
|
||||
|
||||
Performance: 4529.619 tau/day, 5242.615 timesteps/s
|
||||
100.0% CPU use with 4 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Neigh | 0.0049489 | 0.0050479 | 0.0050978 | 0.1 | 0.02
|
||||
Comm | 0.082752 | 0.084491 | 0.085332 | 0.4 | 0.37
|
||||
Output | 0.00054352 | 0.0006034 | 0.00064793 | 0.0 | 0.00
|
||||
Modify | 21.069 | 21.521 | 21.964 | 7.0 | 94.02
|
||||
Other | | 1.278 | | | 5.58
|
||||
|
||||
Nlocal: 1024.00 ave 1050 max 1010 min
|
||||
Histogram: 2 0 0 1 0 0 0 0 0 1
|
||||
Nghost: 0.00000 ave 0 max 0 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 4 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 0
|
||||
Ave neighs/atom = 0.0000000
|
||||
Neighbor list builds = 2169
|
||||
Dangerous builds = 0
|
||||
|
||||
# if you want to check that rotational diffusion is behaving as expected,
|
||||
# uncomment next three lines for dump output and then plot <e(t).e(0)>,
|
||||
# which should decay exponentially with timescale (d-1)*D_r (with d
|
||||
# being simulation dimension)
|
||||
|
||||
#dump 1 all custom 2000 dump_ideal_${params}_3d.lammpstrj id type # x y xu yu mux muy muz fx fy fz
|
||||
#dump_modify 1 first yes sort id
|
||||
|
||||
#run 120000
|
||||
Total wall time: 0:00:30
|
||||
@ -1,54 +0,0 @@
|
||||
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
|
||||
|
||||
|
||||
@ -1,37 +0,0 @@
|
||||
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
|
||||
@ -1,40 +0,0 @@
|
||||
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
|
||||
312
src/USER-BROWNIAN/fix_propel_self.cpp
Normal file
312
src/USER-BROWNIAN/fix_propel_self.cpp
Normal file
@ -0,0 +1,312 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
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!
|
||||
|
||||
Current maintainer: Sam Cameron @ University of Bristol
|
||||
|
||||
----------------------------------------------------------------------- */
|
||||
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
#include <string.h>
|
||||
#include "fix_propel_self.h"
|
||||
#include "math_extra.h"
|
||||
#include "atom.h"
|
||||
#include "atom_vec_ellipsoid.h"
|
||||
#include "force.h"
|
||||
#include "update.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
enum{DIPOLE,VELOCITY,QUAT};
|
||||
|
||||
#define TOL 1e-14
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixPropelSelf::FixPropelSelf(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
|
||||
virial_flag = 1;
|
||||
|
||||
if (narg < 5)
|
||||
error->all(FLERR,"Illegal fix propel/self command");
|
||||
|
||||
|
||||
magnitude = utils::numeric(FLERR,arg[3],false,lmp);
|
||||
|
||||
if (strcmp(arg[4],"velocity") == 0) {
|
||||
mode = VELOCITY;
|
||||
thermo_virial = 0;
|
||||
if (narg != 5) {
|
||||
error->all(FLERR,"Illegal fix propel/self command");
|
||||
}
|
||||
} else if (strcmp(arg[4],"dipole") == 0) {
|
||||
mode = DIPOLE;
|
||||
thermo_virial = 1;
|
||||
if (narg != 5) {
|
||||
error->all(FLERR,"Illegal fix propel/self command");
|
||||
}
|
||||
} else if (strcmp(arg[4],"quat") == 0) {
|
||||
mode = QUAT;
|
||||
thermo_virial = 1;
|
||||
if (narg != 8) {
|
||||
error->all(FLERR,"Illegal fix propel/self command");
|
||||
} else {
|
||||
sx = utils::numeric(FLERR,arg[5],false,lmp);
|
||||
sy = utils::numeric(FLERR,arg[6],false,lmp);
|
||||
sz = utils::numeric(FLERR,arg[7],false,lmp);
|
||||
double qnorm = sqrt(sx*sx + sy*sy + sz*sz);
|
||||
sx = sx/qnorm;
|
||||
sy = sy/qnorm;
|
||||
sz = sz/qnorm;
|
||||
}
|
||||
} else {
|
||||
error->all(FLERR,"Illegal fix propel/self command");
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixPropelSelf::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= POST_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixPropelSelf::~FixPropelSelf()
|
||||
{
|
||||
|
||||
}
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixPropelSelf::init()
|
||||
{
|
||||
if (mode == DIPOLE && !atom->mu_flag)
|
||||
error->all(FLERR,"Fix propel/self requires atom attribute mu "
|
||||
"with option dipole.");
|
||||
|
||||
if (mode == QUAT) {
|
||||
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
|
||||
if (!avec)
|
||||
error->all(FLERR,"Fix propel/self requires "
|
||||
"atom style ellipsoid with option quat.");
|
||||
|
||||
// check that all particles are finite-size ellipsoids
|
||||
// no point particles allowed, spherical is OK
|
||||
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
if (ellipsoid[i] < 0)
|
||||
error->one(FLERR,"Fix propel/self requires extended particles "
|
||||
"with option quat.");
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
void FixPropelSelf::setup(int vflag)
|
||||
{
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
void FixPropelSelf::post_force(int vflag)
|
||||
{
|
||||
if (mode == DIPOLE)
|
||||
post_force_dipole(vflag);
|
||||
else if (mode == VELOCITY)
|
||||
post_force_velocity(vflag);
|
||||
else if (mode == QUAT)
|
||||
post_force_quaternion(vflag);
|
||||
|
||||
}
|
||||
|
||||
|
||||
void FixPropelSelf::post_force_dipole(int vflag)
|
||||
{
|
||||
double **f = atom->f;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
double **x = atom->x;
|
||||
double **mu = atom->mu;
|
||||
double fx,fy,fz;
|
||||
|
||||
|
||||
// energy and virial setup
|
||||
double vi[6];
|
||||
if (vflag) v_setup(vflag);
|
||||
else evflag = 0;
|
||||
|
||||
// if domain has PBC, need to unwrap for virial
|
||||
double unwrap[3];
|
||||
imageint *image = atom->image;
|
||||
|
||||
// Add the active force to the atom force:
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
fx = magnitude*mu[i][0];
|
||||
fy = magnitude*mu[i][1];
|
||||
fz = magnitude*mu[i][2];
|
||||
f[i][0] += fx;
|
||||
f[i][1] += fy;
|
||||
f[i][2] += fz;
|
||||
|
||||
if (evflag) {
|
||||
domain->unmap(x[i],image[i],unwrap);
|
||||
vi[0] = fx*unwrap[0];
|
||||
vi[1] = fy*unwrap[1];
|
||||
vi[2] = fz*unwrap[2];
|
||||
vi[3] = fx*unwrap[1];
|
||||
vi[4] = fx*unwrap[2];
|
||||
vi[5] = fy*unwrap[2];
|
||||
v_tally(i, vi);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void FixPropelSelf::post_force_velocity(int vflag)
|
||||
{
|
||||
double **f = atom->f;
|
||||
double **v = atom->v;
|
||||
double **x = atom->x;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int *type = atom->type;
|
||||
double nv2,fnorm,fx,fy,fz;
|
||||
|
||||
|
||||
// energy and virial setup
|
||||
double vi[6];
|
||||
if (vflag) v_setup(vflag);
|
||||
else evflag = 0;
|
||||
|
||||
// if domain has PBC, need to unwrap for virial
|
||||
double unwrap[3];
|
||||
imageint *image = atom->image;
|
||||
|
||||
// Add the active force to the atom force:
|
||||
for(int i = 0; i < nlocal; ++i) {
|
||||
if( mask[i] & groupbit ){
|
||||
|
||||
nv2 = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
|
||||
fnorm = 0.0;
|
||||
|
||||
if (nv2 > TOL) {
|
||||
|
||||
// Without this check you can run into numerical
|
||||
// issues because fnorm will blow up.
|
||||
|
||||
fnorm = magnitude / sqrt(nv2);
|
||||
}
|
||||
fx = fnorm * v[i][0];
|
||||
fy = fnorm * v[i][1];
|
||||
fz = fnorm * v[i][2];
|
||||
|
||||
f[i][0] += fx;
|
||||
f[i][1] += fy;
|
||||
f[i][2] += fz;
|
||||
|
||||
if (evflag) {
|
||||
domain->unmap(x[i],image[i],unwrap);
|
||||
vi[0] = fx*unwrap[0];
|
||||
vi[1] = fy*unwrap[1];
|
||||
vi[2] = fz*unwrap[2];
|
||||
vi[3] = fx*unwrap[1];
|
||||
vi[4] = fx*unwrap[2];
|
||||
vi[5] = fy*unwrap[2];
|
||||
v_tally(i, vi);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void FixPropelSelf::post_force_quaternion(int vflag)
|
||||
{
|
||||
double **f = atom->f;
|
||||
double **x = atom->x;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int *type = atom->type;
|
||||
int* ellipsoid = atom->ellipsoid;
|
||||
|
||||
// ellipsoidal properties
|
||||
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
|
||||
double f_act[3] = { sx, sy, sz };
|
||||
double f_rot[3];
|
||||
double *quat;
|
||||
double Q[3][3];
|
||||
double fx,fy,fz;
|
||||
|
||||
|
||||
// energy and virial setup
|
||||
double vi[6];
|
||||
if (vflag) v_setup(vflag);
|
||||
else evflag = 0;
|
||||
|
||||
// if domain has PBC, need to unwrap for virial
|
||||
double unwrap[3];
|
||||
imageint *image = atom->image;
|
||||
|
||||
// Add the active force to the atom force:
|
||||
for( int i = 0; i < nlocal; ++i ){
|
||||
if( mask[i] & groupbit ){
|
||||
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
MathExtra::quat_to_mat( quat, Q );
|
||||
MathExtra::matvec( Q, f_act, f_rot );
|
||||
|
||||
fx = magnitude*f_rot[0];
|
||||
fy = magnitude*f_rot[1];
|
||||
fz = magnitude*f_rot[2];
|
||||
|
||||
f[i][0] += fx;
|
||||
f[i][1] += fy;
|
||||
f[i][2] += fz;
|
||||
|
||||
if (evflag) {
|
||||
domain->unmap(x[i],image[i],unwrap);
|
||||
vi[0] = fx*unwrap[0];
|
||||
vi[1] = fy*unwrap[1];
|
||||
vi[2] = fz*unwrap[2];
|
||||
vi[3] = fx*unwrap[1];
|
||||
vi[4] = fx*unwrap[2];
|
||||
vi[5] = fy*unwrap[2];
|
||||
v_tally(i, vi);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -21,51 +21,50 @@ FixStyle(propel/self,FixPropelSelf)
|
||||
#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);
|
||||
void init();
|
||||
void post_force(int);
|
||||
void setup(int);
|
||||
int setmask();
|
||||
|
||||
double memory_usage();
|
||||
|
||||
protected:
|
||||
enum operation_modes {
|
||||
VELOCITY = 0,
|
||||
QUATERNION = 1
|
||||
};
|
||||
|
||||
private:
|
||||
private:
|
||||
double magnitude;
|
||||
double sx,sy,sz;
|
||||
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.
|
||||
void post_force_dipole(int);
|
||||
void post_force_velocity(int);
|
||||
void post_force_quaternion(int);
|
||||
|
||||
class AtomVecEllipsoid *avec;
|
||||
|
||||
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
|
||||
E: Illegal fix propel/self 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.
|
||||
Wrong number/type of input arguments.
|
||||
|
||||
E: Fix propel/self requires atom attribute mu with option dipole.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix propel/self requires atom style ellipsoid with option quat.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
Fix propel/self requires extended particles with option quat.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
*/
|
||||
@ -68,7 +68,6 @@ fix nvk, Efrem Braun (UC Berkeley), efrem.braun at gmail.com, https://github.com
|
||||
fix orient/eco Adrian A. Schratt and Volker Mohles (Ruhr-Uni Bochum), volker.mohles at rub.de, 6 Jun 2020
|
||||
fix pafi, Thomas Swinburne (CNRS), swinburne at cinam.univ-mrs.fr, 1st Sep 2020
|
||||
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
|
||||
|
||||
@ -1,258 +0,0 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://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 "atom.h"
|
||||
#include "atom_vec_ellipsoid.h"
|
||||
#include "error.h"
|
||||
#include "math_const.h"
|
||||
#include "math_extra.h"
|
||||
|
||||
#include <cctype>
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
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(nullptr)
|
||||
{
|
||||
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 = utils::numeric( FLERR, argv[4] ,false,lmp);
|
||||
|
||||
// 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 = utils::inumeric(FLERR,argv[iarg],false,lmp);
|
||||
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;
|
||||
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
int *type = atom->type;
|
||||
int* ellipsoid = atom->ellipsoid;
|
||||
|
||||
AtomVecEllipsoid *av = static_cast<AtomVecEllipsoid*>(atom->style_match("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;
|
||||
}
|
||||
Reference in New Issue
Block a user