Corrected implementation of ellipsoidal dynamics, made ashared base class for the time-integrators, templated the time-integrators (and so reversed changes that this PR had previously made to random_mars src files), combined docs of all three integrators.
This commit is contained in:
Sam Cameron
@ -40,8 +40,8 @@ OPT.
|
||||
* :doc:`aveforce <fix_aveforce>`
|
||||
* :doc:`balance <fix_balance>`
|
||||
* :doc:`brownian <fix_brownian>`
|
||||
* :doc:`brownian/asphere <fix_brownian_asphere>`
|
||||
* :doc:`brownian/sphere <fix_brownian_sphere>`
|
||||
* :doc:`brownian/asphere <fix_brownian>`
|
||||
* :doc:`brownian/sphere <fix_brownian>`
|
||||
* :doc:`bocs <fix_bocs>`
|
||||
* :doc:`bond/break <fix_bond_break>`
|
||||
* :doc:`bond/create <fix_bond_create>`
|
||||
|
||||
@ -183,8 +183,8 @@ accelerated styles exist.
|
||||
* :doc:`aveforce <fix_aveforce>` - add an averaged force to each atom
|
||||
* :doc:`balance <fix_balance>` - perform dynamic load-balancing
|
||||
* :doc:`brownian <fix_brownian>` - overdamped translational brownian motion
|
||||
* :doc:`brownian/asphere <fix_brownian_asphere>` - overdamped translational and rotational brownian motion for ellipsoids
|
||||
* :doc:`brownian/sphere <fix_brownian_sphere>` - overdamped translational and rotational brownian motion for spheres
|
||||
* :doc:`brownian/asphere <fix_brownian>` - overdamped translational and rotational brownian motion for ellipsoids
|
||||
* :doc:`brownian/sphere <fix_brownian>` - overdamped translational and rotational brownian motion for spheres
|
||||
* :doc:`bocs <fix_bocs>` - NPT style time integration with pressure correction
|
||||
* :doc:`bond/break <fix_bond_break>` - break bonds on the fly
|
||||
* :doc:`bond/create <fix_bond_create>` - create bonds on the fly
|
||||
|
||||
@ -1,21 +1,30 @@
|
||||
.. index:: fix brownian
|
||||
.. index:: fix brownian/sphere
|
||||
.. index:: fix brownian/asphere
|
||||
|
||||
fix brownian command
|
||||
====================
|
||||
===========================
|
||||
|
||||
fix brownian/sphere command
|
||||
===========================
|
||||
|
||||
fix brownian/sphere command
|
||||
===========================
|
||||
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID brownian gamma_t diff_t seed keyword args
|
||||
fix ID group-ID style_name temp seed keyword args
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* brownian/sphere = style name of this fix command
|
||||
* gamma_t = translational friction coefficient
|
||||
* diff_t = translational diffusion coefficient
|
||||
* zero or more keyword/value pairs may be appended
|
||||
* keyword = *rng*
|
||||
* style_name = *brownian* or *brownian/sphere* or *brownian/asphere*
|
||||
* temp = temperature
|
||||
* seed = random number generator seed
|
||||
* one or more keyword/value pairs may be appended
|
||||
* keyword = *rng* or *dipole* or *gamma_r_eigen* or *gamma_t_eigen* or *gamma_r* or *gamma_t*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
@ -23,39 +32,79 @@ Syntax
|
||||
*uniform* = use uniform random number generator
|
||||
*gaussian* = use gaussian random number generator
|
||||
*none* = turn off noise
|
||||
*dipole* value = *mux* and *muy* and *muz* for *brownian/asphere*
|
||||
*mux*, *muy*, and *muz* = update orientation of dipole having direction (*mux*,*muy*,*muz*) in body frame of rigid body
|
||||
*gamma_r_eigen* values = *gr1* and *gr2* and *gr3* for *brownian/asphere*
|
||||
*gr1*, *gr2*, and *gr3* = diagonal entries of body frame rotational friction tensor
|
||||
*gamma_r* values = *gr* for *brownian/sphere*
|
||||
*gr* = magnitude of the (isotropic) rotational friction tensor
|
||||
*gamma_t_eigen* values = *gt1* and *gt2* and *gt3* for *brownian/asphere*
|
||||
*gt1*, *gt2*, and *gt3* = diagonal entries of body frame translational friction tensor
|
||||
*gamma_t* values = *gt* for *brownian* and *brownian/sphere*
|
||||
*gt* = magnitude of the (isotropic) translational friction tensor
|
||||
|
||||
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix 1 all brownian 1.0 3.0 1294019
|
||||
fix 1 all brownian 1.0 3.0 19581092 rng none
|
||||
fix 1 all brownian 1.0 3.0 19581092 rng uniform
|
||||
fix 1 all brownian 1.0 3.0 19581092 rng gaussian
|
||||
fix 1 all brownian 1.0 12908410 gamma_t 1.0
|
||||
fix 1 all brownian 1.0 12908410 gamma_t 3.0 rng gaussian
|
||||
fix 1 all brownian/sphere 1.0 1294019 gamma_t 3.0 gamma_r 1.0
|
||||
fix 1 all brownian/sphere 1.0 19581092 gamma_t 1.0 gamma_r 0.3 rng none
|
||||
fix 1 all brownian/asphere 1.0 1294019 gamma_t_eigen 1.0 2.0 3.0 gamma_r_eigen 4.0 7.0 8.0 rng gaussian
|
||||
fix 1 all brownian/asphere 1.0 1294019 gamma_t_eigen 1.0 2.0 3.0 gamma_r_eigen 4.0 7.0 8.0 dipole 1.0 0.0 0.0
|
||||
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Perform Brownian Dynamics integration to update position and velocity
|
||||
of atoms in the group each timestep. Brownian Dynamics uses Newton's laws of
|
||||
Perform Brownian Dynamics integration to update position, velocity,
|
||||
dipole orientation (for spheres) and quaternion orientation (for ellipsoids,
|
||||
with optional dipole update as well) in the group each timestep.
|
||||
Brownian Dynamics uses Newton's laws of
|
||||
motion in the limit that inertial forces are negligible compared to
|
||||
viscous forces. The stochastic equation of motion is
|
||||
viscous forces. The stochastic equation of motion for the centre of mass
|
||||
positions are
|
||||
|
||||
.. math::
|
||||
|
||||
dr = \frac{F}{\gamma_t}dt+\sqrt{2D_t}dW_t, \\
|
||||
d\mathbf{r} = \mathbf{\gamma}_t^{-1}\mathbf{F}dt+\sqrt{2k_BT}\mathbf{\gamma}_t^{-1/2}d\mathbf{W}_t,
|
||||
|
||||
where :math:`dW_t` is a Wiener processes (see e.g. :ref:`(Gardiner) <GardinerC1>`).
|
||||
in the lab-frame (i.e. :math:`\mathbf{\gamma}_t` is not diagonal, but only depends on
|
||||
orientation and so the noise is still additive).
|
||||
|
||||
The rotational motion for the spherical and ellipsoidal particles is not as simple an
|
||||
expression, but is chosen to replicate the Boltzmann distribution for the case of
|
||||
conservative torques (see :ref:`(Ilie) <Ilie1>` or :ref:`(Delong) <Delong1>`).
|
||||
|
||||
For the style *brownian*, only the positions of the particles are updated. This is
|
||||
therefore suitable for point particle simulations.
|
||||
|
||||
For the style *brownian/sphere*, the positions of the particles are updated, and a dipole
|
||||
slaved to the spherical orientation is also updated. This style therefore requires the
|
||||
hybrid atom style :doc:`atom_style dipole <atom_style>` and
|
||||
:doc:`atom_style sphere <atom_style>`.
|
||||
|
||||
For the style *brownian/asphere*, the centre of mass positions and the quaternions of
|
||||
ellipsoidal particles are updated. This fix style is suitable for equations of motion
|
||||
where the rotational and translational friction tensors are diagonalisable in a certain
|
||||
(body) reference frame.
|
||||
|
||||
|
||||
---------
|
||||
|
||||
.. note::
|
||||
This integrator is designed for generic non-equilibrium
|
||||
simulations with additive noise. There are two important cases which
|
||||
(conceptually) reduce the number of free parameters in this fix.
|
||||
(a) In equilibrium simulations
|
||||
(where fluctuation dissipation theorems are obeyed), one can define
|
||||
the thermal energy :math:`k_bT=D_t\gamma_t`.
|
||||
This integrator does not by default assume a relationship between the
|
||||
rotational and translational friction tensors, though such a relationship
|
||||
should exist in the case of no-slip boundary conditions between the particles and
|
||||
the surrounding (implicit) solvent. E.g. in the case of spherical particles,
|
||||
the condition :math:`\gamma_t=3\gamma_r/\sigma^2` must be explicitly
|
||||
accounted for by setting *gamma_t* to 3x and *gamma_r* to x (where
|
||||
:math:`\sigma` is the spherical diameter). A similar (though more complex)
|
||||
relationship holds for ellipsoids and rod-like particles.
|
||||
|
||||
---------
|
||||
|
||||
@ -63,38 +112,52 @@ where :math:`dW_t` is a Wiener processes (see e.g. :ref:`(Gardiner) <GardinerC1>
|
||||
Temperature computation using the :doc:`compute temp <compute_temp>`
|
||||
will not correctly compute temperature of these overdamped dynamics
|
||||
since we are explicitly neglecting inertial effects.
|
||||
See e.g. chapter 6 of :ref:`(Doi) <Doi1>` for more details on this.
|
||||
Temperature is instead defined in terms of the note above (for
|
||||
equilibrium systems).
|
||||
Furthermore, this time integrator does not add the stochastic terms or
|
||||
viscous terms to the force and/or torques. Rather, they are just added
|
||||
in to the equations of motion to update the degrees of freedom.
|
||||
|
||||
---------
|
||||
|
||||
.. note::
|
||||
The diffusion coefficient :math:`D_t` is measured
|
||||
in units of (length*length)/time, where time and length
|
||||
are in the units specified on the :doc:`units <units>` page.
|
||||
Similarly, :math:`\gamma_t` is measured in
|
||||
units of mass/time.
|
||||
|
||||
---------
|
||||
|
||||
If the *rng* keyword is used with the *uniform* value, then the noise
|
||||
is generated from a uniform distribution (see
|
||||
:ref:`(Dunweg) <Dunweg6>` for why this works). This is the same method
|
||||
:ref:`(Dunweg) <Dunweg7>` for why this works). This is the same method
|
||||
of noise generation as used in :doc:`fix_langevin <fix_langevin>`.
|
||||
|
||||
If the *rng* keyword is used with the *gaussian* value, then the noise
|
||||
is generated from a gaussian distribution. Typically this added
|
||||
complexity is unnecessary, and one should be fine using the *uniform*
|
||||
value for reasons argued in :ref:`(Dunweg) <Dunweg6>`.
|
||||
value for reasons argued in :ref:`(Dunweg) <Dunweg7>`.
|
||||
|
||||
If the *rng* keyword is used with the *none* value, then the noise
|
||||
terms are set to zero.
|
||||
|
||||
|
||||
The *gamma_t* keyword sets the (isotropic) translational viscous damping.
|
||||
Required for (and only compatible with) *brownian* and *brownian/sphere*.
|
||||
The units of *gamma_t* are mass/time.
|
||||
|
||||
The *gamma_r* keyword sets the (isotropic) rotational viscous damping.
|
||||
Required for (and only compatible with) *brownian/sphere*.
|
||||
The units of *gamma_r* are mass*length**2/time.
|
||||
|
||||
The *gamma_r_eigen*, and *gamma_t_eigen* keywords are the eigenvalues of
|
||||
the rotational and viscous damping tensors (having the same units as
|
||||
their isotropic counterparts). Required for (and only compatible with)
|
||||
*brownian/asphere*. For a 2D system, the first two values of *gamma_r_eigen*
|
||||
must be inf (only rotation in xy plane), and the third value of *gamma_t_eigen*
|
||||
must be inf (only diffusion in xy plane).
|
||||
|
||||
If the *dipole* keyword is used, then the dipole moments of the particles
|
||||
are updated as described above. Only compatible with *brownian/asphere*
|
||||
(as *brownian/sphere* updates dipoles automatically).
|
||||
|
||||
----------
|
||||
|
||||
.. include:: accel_styles.rst
|
||||
.. note::
|
||||
For style *brownian/asphere*, the components *gamma_t_eigen* =(x,x,x) and
|
||||
*gamma_r_eigen* = (y,y,y), the dynamics will replicate those of the
|
||||
*brownian/sphere* style with *gamma_t* = x and *gamma_r* = y.
|
||||
|
||||
----------
|
||||
|
||||
@ -105,10 +168,6 @@ No information about this fix is written to :doc:`binary restart files <restart>
|
||||
No global or per-atom quantities are stored
|
||||
by this fix for access by various :doc:`output commands <Howto_output>`.
|
||||
|
||||
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 no*.
|
||||
|
||||
No parameter of this fix can be used with the *start/stop* keywords of
|
||||
the :doc:`run <run>` command. This fix is not invoked during
|
||||
@ -117,6 +176,13 @@ the :doc:`run <run>` command. This fix is not invoked during
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
The style *brownian/sphere* fix requires that atoms store torque and angular velocity (omega)
|
||||
as defined by the :doc:`atom_style sphere <atom_style>` command.
|
||||
The style *brownian/asphere* fix requires that atoms store torque and quaternions
|
||||
as defined by the :doc:`atom_style ellipsoid <atom_style>` command.
|
||||
If the *dipole* keyword is used, they must also store a dipole moment
|
||||
as defined by the :doc:`atom_style dipole <atom_style>` command.
|
||||
|
||||
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.
|
||||
@ -125,26 +191,28 @@ doc page for more info.
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`fix propel/self <fix_propel_self>`,
|
||||
:doc:`fix langevin <fix_langevin>`, :doc:`fix nve/sphere <fix_nve_sphere>`,
|
||||
:doc:`fix brownian/sphere <fix_brownian_sphere>`,
|
||||
:doc:`fix brownian/asphere <fix_brownian_asphere>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
The default for *rng* is *uniform*.
|
||||
The default for *rng* is *uniform*. The default for the rotational and translational friction
|
||||
tensors are the identity tensor.
|
||||
|
||||
----------
|
||||
|
||||
.. _GardinerC1:
|
||||
|
||||
**(Gardiner)** Gardiner, A Handbook for the Natural and Social Sciences 4th Ed. (2009).
|
||||
.. _Ilie1:
|
||||
|
||||
.. _Doi1:
|
||||
**(Ilie)** Ilie, Briels, den Otter, Journal of Chemical Physics, 142, 114103 (2015).
|
||||
|
||||
**(Doi)** Doi, Soft Matter Physics (2013).
|
||||
|
||||
.. _Dunweg6:
|
||||
.. _Delong1:
|
||||
|
||||
**(Delong)** Delong, Usabiaga, Donev, Journal of Chemical Physics. 143, 144107 (2015)
|
||||
|
||||
.. _Dunweg7:
|
||||
|
||||
**(Dunweg)** Dunweg and Paul, Int J of Modern Physics C, 2, 817-27 (1991).
|
||||
|
||||
|
||||
@ -1,150 +0,0 @@
|
||||
.. index:: fix brownian/asphere
|
||||
|
||||
fix brownian/asphere command
|
||||
============================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID brownian/asphere gamma_t gamma_r diff_t diff_r seed keyword args
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* brownian/asphere = style name of this fix command
|
||||
* gamma_t = translational friction coefficient
|
||||
* gamma_r = rotational friction coefficient
|
||||
* diff_t = translational diffusion coefficient
|
||||
* diff_r = rotational diffusion coefficient
|
||||
* zero or more keyword/value pairs may be appended
|
||||
* keyword = *rng* or *dipole*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*rng* value = *uniform* or *gaussian* or *none*
|
||||
*uniform* = use uniform random number generator
|
||||
*gaussian* = use gaussian random number generator
|
||||
*none* = turn off noise
|
||||
*dipole* value = none = update orientation of dipoles during integration
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix 1 all brownian/asphere 1.0 1.0 1.0 1.0 1294019
|
||||
fix 1 all brownian/asphere 1.0 1.0 1.0 1.0 19581092 rng none dipole
|
||||
fix 1 all brownian/asphere 1.0 1.0 1.0 1.0 19581092 rng uniform
|
||||
fix 1 all brownian/asphere 1.0 1.0 1.0 1.0 19581092 dipole rng gaussian
|
||||
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Perform Brownian Dynamics integration to update position, velocity,
|
||||
angular velocity, particle orientation, and dipole moment for
|
||||
finite-size elipsoidal particles in the group each timestep.
|
||||
Brownian Dynamics uses Newton's laws of
|
||||
motion in the limit that inertial forces are negligible compared to
|
||||
viscous forces. The stochastic equations of motion are
|
||||
|
||||
.. math::
|
||||
|
||||
dr = \frac{F}{\gamma_t}dt+\sqrt{2D_t}dW_t, \\
|
||||
d\Omega = \frac{T}{\gamma_r}dt + \sqrt{2D_r}dW_r,
|
||||
|
||||
where :math:`d\Omega` is an infinitesimal rotation vector (see e.g.
|
||||
Chapter 4 of :ref:`(Goldstein) <GoldsteinCM2>`), :math:`dW_t` and
|
||||
:math:`dW_r` are Wiener processes (see e.g. :ref:`(Gardiner) <GardinerC3>`).
|
||||
The quaternions :math:`q` of the ellipsoid are updated each timestep from
|
||||
the angular velocity vector.
|
||||
|
||||
See :doc:`fix brownian/sphere <fix_brownian_sphere>` for discussion on the
|
||||
values of :math:`\gamma_t`, :math:`\gamma_r`, :math:`D_t`,
|
||||
:math:`D_r`, and temperature when simulating equilibrium systems.
|
||||
|
||||
|
||||
If the *rng* keyword is used with the *uniform* value, then the noise
|
||||
is generated from a uniform distribution (see
|
||||
:ref:`(Dunweg) <Dunweg8>` for why this works). This is the same method
|
||||
of noise generation as used in :doc:`fix_langevin <fix_langevin>`.
|
||||
|
||||
If the *rng* keyword is used with the *gaussian* value, then the noise
|
||||
is generated from a gaussian distribution. Typically this added
|
||||
complexity is unnecessary, and one should be fine using the *uniform*
|
||||
value for reasons argued in :ref:`(Dunweg) <Dunweg8>`.
|
||||
|
||||
If the *rng* keyword is used with the *none* value, then the noise
|
||||
terms are set to zero.
|
||||
|
||||
If the *dipole* keyword is used, then the dipole moments of the particles
|
||||
are updated by setting them along the x axis of the ellipsoidal frames of
|
||||
reference.
|
||||
|
||||
----------
|
||||
|
||||
.. include:: accel_styles.rst
|
||||
|
||||
----------
|
||||
|
||||
Restart, fix_modify, output, run start/stop, minimize info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
No information about this fix is written to :doc:`binary restart files <restart>`.
|
||||
No global or per-atom quantities are stored
|
||||
by this fix for access by various :doc:`output commands <Howto_output>`.
|
||||
|
||||
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 no*.
|
||||
|
||||
No parameter of this fix can be used with the *start/stop* keywords of
|
||||
the :doc:`run <run>` command. This fix is not invoked during
|
||||
:doc:`energy minimization <minimize>`.
|
||||
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This fix requires that atoms store torque and angular velocity (omega)
|
||||
as defined by the :doc:`atom_style sphere <atom_style>` command, as well
|
||||
as atoms which have a definite orientation as defined by the
|
||||
:doc:`atom_style ellipsoid <atom_style>` command.
|
||||
Optionally, they can also store a dipole moment as defined by the
|
||||
:doc:`atom_style dipole <atom_style>` command.
|
||||
|
||||
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.
|
||||
|
||||
All particles in the group must be finite-size ellipsoids. They cannot
|
||||
be point particles.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`fix brownian <fix_brownian>`, :doc:`fix brownian/sphere <fix_brownian_sphere>`,
|
||||
:doc:`fix propel/self <fix_propel_self>`, :doc:`fix langevin <fix_langevin>`,
|
||||
:doc:`fix nve/asphere <fix_nve_asphere>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
The default for *rng* is *uniform*.
|
||||
|
||||
----------
|
||||
|
||||
.. _GoldsteinCM2:
|
||||
|
||||
**(Goldstein)** Goldstein, Poole, and Safko, Classical Mechanics, 3rd Ed. (2001).
|
||||
|
||||
.. _GardinerC3:
|
||||
|
||||
**(Gardiner)** Gardiner, A Handbook for the Natural and Social Sciences 4th Ed. (2009).
|
||||
|
||||
.. _Dunweg8:
|
||||
|
||||
**(Dunweg)** Dunweg and Paul, Int J of Modern Physics C, 2, 817-27 (1991).
|
||||
|
||||
|
||||
@ -1,186 +0,0 @@
|
||||
.. index:: fix brownian/sphere
|
||||
|
||||
fix brownian/sphere command
|
||||
===========================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID brownian/sphere gamma_t gamma_r diff_t diff_r seed keyword args
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* brownian/sphere = style name of this fix command
|
||||
* gamma_t = translational friction coefficient
|
||||
* gamma_r = rotational friction coefficient
|
||||
* diff_t = translational diffusion coefficient
|
||||
* diff_r = rotational diffusion coefficient
|
||||
* zero or more keyword/value pairs may be appended
|
||||
* keyword = *rng* or *dipole*
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
*rng* value = *uniform* or *gaussian* or *none*
|
||||
*uniform* = use uniform random number generator
|
||||
*gaussian* = use gaussian random number generator
|
||||
*none* = turn off noise
|
||||
*dipole* value = none = update orientation of dipoles during integration
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
fix 1 all brownian/sphere 1.0 1.0 1.0 1.0 1294019
|
||||
fix 1 all brownian/sphere 1.0 1.0 1.0 1.0 19581092 rng none dipole
|
||||
fix 1 all brownian/sphere 1.0 1.0 1.0 1.0 19581092 rng uniform
|
||||
fix 1 all brownian/sphere 1.0 1.0 1.0 1.0 19581092 dipole rng gaussian
|
||||
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
Perform Brownian Dynamics integration to update position, velocity,
|
||||
angular velocity, and dipole moment for finite-size spherical particles
|
||||
in the group each timestep. Brownian Dynamics uses Newton's laws of
|
||||
motion in the limit that inertial forces are negligible compared to
|
||||
viscous forces. The stochastic equations of motion are
|
||||
|
||||
.. math::
|
||||
|
||||
dr = \frac{F}{\gamma_t}dt+\sqrt{2D_t}dW_t, \\
|
||||
d\Omega = \frac{T}{\gamma_r}dt + \sqrt{2D_r}dW_r,
|
||||
|
||||
where :math:`d\Omega` is an infinitesimal rotation vector (see e.g.
|
||||
Chapter 4 of :ref:`(Goldstein) <GoldsteinCM1>`), :math:`dW_t` and
|
||||
:math:`dW_r` are Wiener processes (see e.g. :ref:`(Gardiner) <GardinerC2>`).
|
||||
The dipole vectors :math:`e_i` are updated using the rotation matrix
|
||||
|
||||
.. math::
|
||||
|
||||
e_i(t+dt) = e^{\theta_X} e_i(t),\\
|
||||
|
||||
where :math:`\omega=d\Omega/dt` is the angular velocity,
|
||||
:math:`\Delta\theta=|\omega|dt` is the rotation angle about
|
||||
the :math:`\omega` axis, and
|
||||
:math:`(\theta_X)_{ij}=-\epsilon_{ijk}d\Omega_k` is the
|
||||
infinitesimal rotation matrix (see e.g. :ref:`(Callegari) <Callegari1>`,
|
||||
section 7.4).
|
||||
|
||||
.. note::
|
||||
This integrator is designed for generic non-equilibrium
|
||||
simulations with additive noise. There are two important cases which
|
||||
(conceptually) reduce the number of free parameters in this fix.
|
||||
(a) In equilibrium simulations
|
||||
(where fluctuation dissipation theorems are obeyed), one can define
|
||||
the thermal energy :math:`k_bT=D_t\gamma_t=D_r\gamma_r`.
|
||||
(b) When a no-slip boundary condition is expected between the spheres and
|
||||
the surrounding medium, the condition
|
||||
:math:`\gamma_t=3\gamma_r/\sigma^2` must be explicitly
|
||||
accounted for (e.g. by setting *gamma_t* to 3 and *gamma_r* to 1) where
|
||||
:math:`sigma` is the particle diameter.
|
||||
If both (a) and (b) are true, then one must ensure this explicitly via
|
||||
the above relationships.
|
||||
|
||||
---------
|
||||
|
||||
See note on the unphysical result of using :doc:`compute temp <compute_temp>`
|
||||
with this fix in :doc:`fix brownian <fix_brownian>`.
|
||||
|
||||
---------
|
||||
|
||||
.. note::
|
||||
The diffusion coefficient :math:`D_t` and the translational
|
||||
drag coefficient :math:`\gamma_t` are discussed in
|
||||
:doc:`fix brownian <fix_brownian>`. The diffusion coefficient
|
||||
:math:`D_r` is measured in units of 1/time, where time is in the
|
||||
units specified on the :doc:`units <units>` page. Similarly,
|
||||
and :math:`\gamma_r` is measured in
|
||||
units of mass/time and (mass*length*length)/(time).
|
||||
|
||||
---------
|
||||
|
||||
If the *rng* keyword is used with the *uniform* value, then the noise
|
||||
is generated from a uniform distribution (see
|
||||
:ref:`(Dunweg) <Dunweg7>` for why this works). This is the same method
|
||||
of noise generation as used in :doc:`fix_langevin <fix_langevin>`.
|
||||
|
||||
If the *rng* keyword is used with the *gaussian* value, then the noise
|
||||
is generated from a gaussian distribution. Typically this added
|
||||
complexity is unnecessary, and one should be fine using the *uniform*
|
||||
value for reasons argued in :ref:`(Dunweg) <Dunweg7>`.
|
||||
|
||||
If the *rng* keyword is used with the *none* value, then the noise
|
||||
terms are set to zero.
|
||||
|
||||
If the *dipole* keyword is used, then the dipole moments of the particles
|
||||
are updated as described above.
|
||||
|
||||
----------
|
||||
|
||||
.. include:: accel_styles.rst
|
||||
|
||||
----------
|
||||
|
||||
Restart, fix_modify, output, run start/stop, minimize info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
No information about this fix is written to :doc:`binary restart files <restart>`.
|
||||
No global or per-atom quantities are stored
|
||||
by this fix for access by various :doc:`output commands <Howto_output>`.
|
||||
|
||||
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 no*.
|
||||
|
||||
No parameter of this fix can be used with the *start/stop* keywords of
|
||||
the :doc:`run <run>` command. This fix is not invoked during
|
||||
:doc:`energy minimization <minimize>`.
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This fix requires that atoms store torque and angular velocity (omega)
|
||||
as defined by the :doc:`atom_style sphere <atom_style>` command.
|
||||
If the *dipole* keyword is used, they must also store a dipole moment
|
||||
as defined by the :doc:`atom_style dipole <atom_style>` command.
|
||||
|
||||
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 brownian <fix_brownian>`, :doc:`fix brownian/asphere <fix_brownian_asphere>`,
|
||||
:doc:`fix propel/self <fix_propel_self>`,
|
||||
:doc:`fix langevin <fix_langevin>`, :doc:`fix nve/sphere <fix_nve_sphere>`,
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
The default for *rng* is *uniform*.
|
||||
|
||||
----------
|
||||
|
||||
.. _GoldsteinCM1:
|
||||
|
||||
**(Goldstein)** Goldstein, Poole, and Safko, Classical Mechanics, 3rd Ed. (2001).
|
||||
|
||||
.. _GardinerC2:
|
||||
|
||||
**(Gardiner)** Gardiner, A Handbook for the Natural and Social Sciences 4th Ed. (2009).
|
||||
|
||||
.. _Callegari1:
|
||||
|
||||
**(Callegari)** Callegari and Volpe, *Numerical Simulations of Active Brownian
|
||||
Particles*, Flowing Matter, 211-238 (2019).
|
||||
|
||||
.. _Dunweg7:
|
||||
|
||||
**(Dunweg)** Dunweg and Paul, Int J of Modern Physics C, 2, 817-27 (1991).
|
||||
|
||||
|
||||
@ -1,12 +1,11 @@
|
||||
# 2d overdamped brownian dynamics with self-propulsion
|
||||
# force in direction of velocity.
|
||||
##### 2d overdamped brownian dynamics with self-propulsion force in direction of velocity. #####
|
||||
|
||||
variable gamma_t equal 1.0
|
||||
variable D_t equal 1.0
|
||||
variable temp equal 1.0
|
||||
variable seed equal 1974019
|
||||
variable fp equal 4.0
|
||||
|
||||
variable params string ${gamma_t}_${D_t}_${fp}
|
||||
variable params string ${gamma_t}_${temp}_${fp}
|
||||
|
||||
|
||||
log log_${params}_2d.lammps.log
|
||||
@ -30,7 +29,7 @@ neigh_modify every 1 delay 1 check yes
|
||||
|
||||
pair_style none
|
||||
|
||||
fix step all brownian ${gamma_t} ${D_t} ${seed}
|
||||
fix step all brownian ${temp} ${seed} gamma_t ${gamma_t}
|
||||
fix vel all propel/self velocity ${fp}
|
||||
fix 2 all enforce2d
|
||||
fix_modify vel virial yes
|
||||
@ -23,16 +23,12 @@ neigh_modify every 1 delay 1 check yes
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
|
||||
#compute d all property/atom mux muy muz
|
||||
|
||||
fix step all brownian ${gamma_t} ${D_t} ${seed}
|
||||
fix step all brownian 1 ${D_t} ${seed}
|
||||
fix step all brownian 1 1 ${seed}
|
||||
fix step all brownian 1 1 1974019
|
||||
fix vel all propel/self ${fp} velocity
|
||||
fix vel all propel/self 4 velocity
|
||||
fix step all brownian ${temp} ${seed} gamma_t ${gamma_t}
|
||||
fix step all brownian 1 ${seed} gamma_t ${gamma_t}
|
||||
fix step all brownian 1 1974019 gamma_t ${gamma_t}
|
||||
fix step all brownian 1 1974019 gamma_t 1
|
||||
fix vel all propel/self velocity ${fp}
|
||||
fix vel all propel/self velocity 4
|
||||
fix 2 all enforce2d
|
||||
fix_modify vel virial yes
|
||||
|
||||
@ -49,10 +45,10 @@ WARNING: Communication cutoff is 0.0. No ghost atoms will be generated. Atoms ma
|
||||
Per MPI rank memory allocation (min/avg/max) = 2.289 | 2.289 | 2.289 Mbytes
|
||||
Step Temp E_pair c_press
|
||||
0 1 0 -0.18336111
|
||||
50000 1.9663098e+10 0 -0.75033044
|
||||
Loop time of 2.45902 on 1 procs for 50000 steps with 1024 atoms
|
||||
50000 2.0329082e+10 0 -0.30450754
|
||||
Loop time of 1.76045 on 1 procs for 50000 steps with 1024 atoms
|
||||
|
||||
Performance: 0.176 tau/day, 20333.276 timesteps/s
|
||||
Performance: 0.245 tau/day, 28401.873 timesteps/s
|
||||
100.0% CPU use with 1 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
@ -60,10 +56,10 @@ 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.0078395 | 0.0078395 | 0.0078395 | 0.0 | 0.32
|
||||
Output | 2.2947e-05 | 2.2947e-05 | 2.2947e-05 | 0.0 | 0.00
|
||||
Modify | 2.332 | 2.332 | 2.332 | 0.0 | 94.83
|
||||
Other | | 0.1192 | | | 4.85
|
||||
Comm | 0.008736 | 0.008736 | 0.008736 | 0.0 | 0.50
|
||||
Output | 2.3562e-05 | 2.3562e-05 | 2.3562e-05 | 0.0 | 0.00
|
||||
Modify | 1.6284 | 1.6284 | 1.6284 | 0.0 | 92.50
|
||||
Other | | 0.1233 | | | 7.00
|
||||
|
||||
Nlocal: 1024.00 ave 1024 max 1024 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
@ -96,141 +92,141 @@ 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) = 2.664 | 2.664 | 2.664 Mbytes
|
||||
Step Temp E_pair c_msd[1] c_msd[2] c_msd[3] c_msd[4] c_press
|
||||
0 1.9663098e+10 0 0 0 0 0 -0.75033044
|
||||
1000 199346.07 0 0.01933096 0.020555579 0 0.039886539 -0.45701943
|
||||
2000 198310 0 0.040165459 0.041283119 0 0.081448577 0.34264096
|
||||
3000 204115.93 0 0.057654441 0.060411193 0 0.11806563 -0.41710363
|
||||
4000 196903.18 0 0.075874071 0.08470884 0 0.16058291 -0.18627889
|
||||
5000 201382.13 0 0.097484871 0.1049401 0 0.20242497 -0.12876225
|
||||
6000 195317.96 0 0.11872475 0.12469358 0 0.24341834 -0.3084651
|
||||
7000 192139.34 0 0.14148561 0.14363452 0 0.28512013 0.0032867364
|
||||
8000 201737.84 0 0.16055109 0.16405717 0 0.32460826 0.36435453
|
||||
9000 199390.58 0 0.17897382 0.18795928 0 0.3669331 0.025659298
|
||||
10000 207807.41 0 0.19680417 0.20733821 0 0.40414239 -0.35368379
|
||||
11000 201936.29 0 0.21666426 0.22702132 0 0.44368558 0.044253449
|
||||
12000 196863.68 0 0.2394452 0.24672678 0 0.48617198 0.059027892
|
||||
13000 199368.34 0 0.26647368 0.2700584 0 0.53653208 0.27090461
|
||||
14000 201246.39 0 0.28799289 0.29823282 0 0.58622571 0.59883778
|
||||
15000 195355.47 0 0.29975278 0.32348787 0 0.62324065 -0.70763643
|
||||
16000 198372.41 0 0.32191014 0.34434864 0 0.66625878 -0.36543908
|
||||
17000 193442.08 0 0.33927003 0.36811239 0 0.70738242 0.28541473
|
||||
18000 197441 0 0.36067818 0.38982011 0 0.75049829 -0.45670227
|
||||
19000 208769.5 0 0.37965583 0.41015661 0 0.78981244 -0.47803396
|
||||
20000 198311.2 0 0.3968078 0.42701175 0 0.82381955 -0.18642397
|
||||
21000 201365.22 0 0.4151043 0.44909345 0 0.86419775 0.86839756
|
||||
22000 198253.24 0 0.4396583 0.46388261 0 0.90354091 -0.19592545
|
||||
23000 204598.51 0 0.45382292 0.49253671 0 0.94635963 -0.24169987
|
||||
24000 211421.88 0 0.46086338 0.51831304 0 0.97917642 0.49751915
|
||||
25000 198690.71 0 0.47110913 0.53640271 0 1.0075118 -0.24475563
|
||||
26000 203981.49 0 0.49265476 0.55310571 0 1.0457605 0.20237224
|
||||
27000 201128.99 0 0.52865208 0.57516064 0 1.1038127 -0.40826104
|
||||
28000 198529.77 0 0.54479087 0.59876678 0 1.1435576 0.41576857
|
||||
29000 205024.27 0 0.56195744 0.61217109 0 1.1741285 -0.79146635
|
||||
30000 201565.62 0 0.58276132 0.63743585 0 1.2201972 -0.065832917
|
||||
31000 197893.43 0 0.61132665 0.66126375 0 1.2725904 0.47907079
|
||||
32000 201395.11 0 0.61904956 0.67520462 0 1.2942542 -0.7408472
|
||||
33000 202064.22 0 0.64760511 0.69087605 0 1.3384812 -0.14601514
|
||||
34000 191227.75 0 0.65698736 0.73857849 0 1.3955659 0.2177548
|
||||
35000 199474.65 0 0.66491543 0.76604575 0 1.4309612 -0.64627039
|
||||
36000 195252.77 0 0.67565581 0.79911139 0 1.4747672 0.0293298
|
||||
37000 198167.14 0 0.68899202 0.81268008 0 1.5016721 -0.0055245918
|
||||
38000 202995.35 0 0.70224976 0.82436547 0 1.5266152 -0.2826768
|
||||
39000 197129.03 0 0.71270072 0.8579444 0 1.5706451 0.063623666
|
||||
40000 199153.69 0 0.73777312 0.88820969 0 1.6259828 -0.26740551
|
||||
41000 205347.31 0 0.75613153 0.9006214 0 1.6567529 0.82415354
|
||||
42000 199423.73 0 0.76864739 0.92457092 0 1.6932183 -0.16636304
|
||||
43000 198052 0 0.79841199 0.93832523 0 1.7367372 -0.016241224
|
||||
44000 205205.39 0 0.81727188 0.96823569 0 1.7855076 -0.10405924
|
||||
45000 199116.46 0 0.83208052 0.99352694 0 1.8256075 0.040835286
|
||||
46000 198759.48 0 0.84291542 1.0038949 0 1.8468104 0.46109436
|
||||
47000 189676.9 0 0.86430719 1.0131299 0 1.8774371 -0.67947102
|
||||
48000 199801.4 0 0.90662572 1.0286589 0 1.9352846 -0.5293946
|
||||
49000 199730.89 0 0.93530132 1.0568273 0 1.9921286 -0.27562311
|
||||
50000 203272.56 0 0.96366375 1.0790026 0 2.0426664 -0.10629234
|
||||
51000 196992.09 0 0.97818106 1.1030549 0 2.0812359 0.31719382
|
||||
52000 204063.62 0 1.0068773 1.1239506 0 2.130828 -7.1998264e-05
|
||||
53000 204349.41 0 1.0277098 1.1546477 0 2.1823575 0.16897786
|
||||
54000 203207.3 0 1.0415673 1.1881409 0 2.2297082 -0.25033492
|
||||
55000 194841.07 0 1.0600954 1.2179033 0 2.2779987 0.0062433629
|
||||
56000 198601.58 0 1.071562 1.2363958 0 2.3079578 -0.13124124
|
||||
57000 196654.6 0 1.0997086 1.2486573 0 2.3483659 -0.46425912
|
||||
58000 197781.48 0 1.112526 1.2706195 0 2.3831455 0.26007922
|
||||
59000 199853.9 0 1.1295132 1.2978402 0 2.4273533 -0.030877018
|
||||
60000 201698.14 0 1.1690892 1.3228782 0 2.4919675 -0.23190043
|
||||
61000 199260.52 0 1.1870513 1.3431945 0 2.5302458 -0.040373885
|
||||
62000 196909.15 0 1.2007194 1.3683525 0 2.5690719 0.21318495
|
||||
63000 203927.79 0 1.2442809 1.3964232 0 2.6407041 0.10156282
|
||||
64000 205322.27 0 1.2721207 1.4159422 0 2.6880629 0.33393307
|
||||
65000 199142.6 0 1.2989685 1.4330729 0 2.7320414 -0.65644005
|
||||
66000 205023.06 0 1.2948208 1.4255094 0 2.7203302 0.22290721
|
||||
67000 209750.01 0 1.3127511 1.4369628 0 2.7497139 -0.40241279
|
||||
68000 200205.19 0 1.3355277 1.4541568 0 2.7896846 0.6665415
|
||||
69000 198653.01 0 1.3500764 1.4962697 0 2.8463461 -0.28396983
|
||||
70000 207485.71 0 1.3825583 1.5095552 0 2.8921135 0.34774599
|
||||
71000 203918.68 0 1.4090995 1.5246756 0 2.9337751 0.071958672
|
||||
72000 199038.47 0 1.4246969 1.5612784 0 2.9859753 -0.42129173
|
||||
73000 197380.7 0 1.445835 1.6025372 0 3.0483722 -0.099854135
|
||||
74000 205006.49 0 1.4703253 1.606784 0 3.0771093 0.137244
|
||||
75000 196040.42 0 1.4897388 1.6297195 0 3.1194583 -0.46715263
|
||||
76000 200022.09 0 1.5178017 1.6560075 0 3.1738092 -0.21504553
|
||||
77000 200766.4 0 1.5184584 1.6673791 0 3.1858374 0.54447858
|
||||
78000 199636.76 0 1.5344452 1.6845098 0 3.2189549 0.021903761
|
||||
79000 204188.86 0 1.5567356 1.7205197 0 3.2772553 0.15356898
|
||||
80000 199862.31 0 1.5669731 1.7265239 0 3.2934969 -0.26342032
|
||||
81000 198557.57 0 1.5735674 1.7468943 0 3.3204617 0.22068216
|
||||
82000 203675.4 0 1.5898596 1.7646027 0 3.3544624 0.51221445
|
||||
83000 203412.56 0 1.6066548 1.7813624 0 3.3880172 0.33512263
|
||||
84000 200896.36 0 1.6112635 1.812552 0 3.4238154 1.1657793
|
||||
85000 198987.17 0 1.652135 1.8336748 0 3.4858098 -0.20419704
|
||||
86000 203027.04 0 1.6728041 1.864244 0 3.5370481 0.61746313
|
||||
87000 203913.8 0 1.6783324 1.8762967 0 3.554629 0.28428076
|
||||
88000 205061.23 0 1.6781682 1.879458 0 3.5576262 0.089285353
|
||||
89000 198210.89 0 1.7017682 1.9029345 0 3.6047027 -0.23977904
|
||||
90000 196170.91 0 1.7291253 1.9258436 0 3.6549689 0.15806438
|
||||
91000 202846.79 0 1.7592339 1.9431644 0 3.7023983 0.17723099
|
||||
92000 198976.91 0 1.762318 1.9742003 0 3.7365183 0.25668658
|
||||
93000 194578.12 0 1.7880716 2.0009816 0 3.7890532 0.20231476
|
||||
94000 200319.82 0 1.7854494 2.020855 0 3.8063044 0.14729427
|
||||
95000 202430.33 0 1.7925005 2.0480213 0 3.8405218 -0.28291245
|
||||
96000 200145.58 0 1.8113602 2.0621043 0 3.8734644 0.05547277
|
||||
97000 194617.53 0 1.8286924 2.0729365 0 3.9016289 -0.59829377
|
||||
98000 200829.37 0 1.8485835 2.077731 0 3.9263145 0.78242718
|
||||
99000 198197.55 0 1.8537119 2.0873455 0 3.9410573 -0.52970118
|
||||
100000 204149.63 0 1.8897394 2.0942927 0 3.9840321 0.58118967
|
||||
101000 198654.59 0 1.9126448 2.1380708 0 4.0507156 -0.61766977
|
||||
102000 201198.05 0 1.9433521 2.143049 0 4.0864011 0.15570108
|
||||
103000 200383.43 0 1.9669578 2.1361518 0 4.1031095 -0.10556532
|
||||
104000 204363.26 0 1.9827272 2.1564095 0 4.1391367 -0.060748593
|
||||
105000 198267.87 0 1.9904164 2.1804141 0 4.1708305 0.40995747
|
||||
106000 202082.12 0 1.993709 2.1925288 0 4.1862378 -0.52458386
|
||||
107000 200209.49 0 2.015427 2.219161 0 4.234588 -0.67350679
|
||||
108000 203112.58 0 2.0467303 2.2405372 0 4.2872675 0.25033168
|
||||
109000 203844.13 0 2.056314 2.2623929 0 4.3187068 0.3384149
|
||||
110000 201740.33 0 2.0706519 2.285547 0 4.3561989 0.35582365
|
||||
111000 202493.05 0 2.0725826 2.308685 0 4.3812676 -0.26487212
|
||||
112000 205404.09 0 2.0969613 2.3252767 0 4.422238 0.58346057
|
||||
113000 201223.54 0 2.0876103 2.3316941 0 4.4193044 0.37747414
|
||||
114000 208649.11 0 2.1095116 2.3488008 0 4.4583124 0.068129648
|
||||
115000 202708.5 0 2.1233837 2.3701129 0 4.4934966 -0.23734073
|
||||
116000 202482.42 0 2.1221907 2.4262516 0 4.5484424 -0.087142119
|
||||
117000 200384.65 0 2.1487723 2.4619437 0 4.6107161 -0.13673271
|
||||
118000 196885.36 0 2.158057 2.4818335 0 4.6398905 0.31912412
|
||||
119000 194064.42 0 2.1821315 2.5032336 0 4.6853651 0.17615749
|
||||
120000 195203.09 0 2.1939678 2.539838 0 4.7338058 0.5106086
|
||||
Loop time of 6.05038 on 1 procs for 120000 steps with 1024 atoms
|
||||
0 2.0329082e+10 0 0 0 0 0 -0.30450754
|
||||
1000 194002.24 0 0.019541352 0.020412616 0 0.039953968 0.17960144
|
||||
2000 201170.68 0 0.037773433 0.040755027 0 0.078528461 0.030184783
|
||||
3000 200206.8 0 0.058949662 0.060726886 0 0.11967655 0.080592855
|
||||
4000 205243 0 0.084357492 0.084713582 0 0.16907107 -0.9284229
|
||||
5000 201741.41 0 0.1024163 0.10453135 0 0.20694765 0.31035655
|
||||
6000 203143.61 0 0.12337925 0.12397645 0 0.2473557 -0.88772072
|
||||
7000 202627.98 0 0.14537908 0.14173625 0 0.28711534 0.17510592
|
||||
8000 200688.93 0 0.15993736 0.16450382 0 0.32444118 0.20926033
|
||||
9000 200360.25 0 0.18277138 0.18597743 0 0.36874881 0.14728415
|
||||
10000 203691.46 0 0.20588386 0.20467547 0 0.41055933 -0.041058877
|
||||
11000 203028.99 0 0.2282616 0.21792715 0 0.44618875 0.14648625
|
||||
12000 195738.49 0 0.24584925 0.22877733 0 0.47462658 0.11242329
|
||||
13000 198977.93 0 0.26808832 0.25467835 0 0.52276667 -0.82500263
|
||||
14000 193395.28 0 0.29074728 0.27404624 0 0.56479352 0.29664298
|
||||
15000 200408 0 0.30241331 0.29063643 0 0.59304975 0.22506685
|
||||
16000 194570.8 0 0.32914382 0.30896611 0 0.63810992 -0.011781021
|
||||
17000 201699.18 0 0.34741469 0.33153019 0 0.67894488 -0.073195501
|
||||
18000 198665.14 0 0.35708615 0.35032716 0 0.7074133 0.33684482
|
||||
19000 201011.15 0 0.37328322 0.36518777 0 0.73847099 -0.20884796
|
||||
20000 201387.99 0 0.39484415 0.3783921 0 0.77323625 -0.30200013
|
||||
21000 202952.38 0 0.41838333 0.3958867 0 0.81427003 -0.98801165
|
||||
22000 205150.98 0 0.44071908 0.42001514 0 0.86073422 -0.65548109
|
||||
23000 200164.9 0 0.46156041 0.44411636 0 0.90567677 0.24186828
|
||||
24000 196027.89 0 0.4834982 0.46586728 0 0.94936548 -0.43593925
|
||||
25000 198438.08 0 0.50328321 0.4865808 0 0.98986402 -0.25516412
|
||||
26000 200569.69 0 0.52535125 0.52102617 0 1.0463774 -0.41486705
|
||||
27000 200598.31 0 0.52958053 0.53142647 0 1.061007 0.011109738
|
||||
28000 205592.24 0 0.54694416 0.5464941 0 1.0934383 -0.18134038
|
||||
29000 199085.55 0 0.57222754 0.57188852 0 1.1441161 -0.46788121
|
||||
30000 201398.39 0 0.59813146 0.58520822 0 1.1833397 -0.15341079
|
||||
31000 197674.44 0 0.61191002 0.60572867 0 1.2176387 0.25394145
|
||||
32000 205188.42 0 0.63435682 0.63745001 0 1.2718068 0.26362999
|
||||
33000 202316.72 0 0.66509392 0.65462302 0 1.3197169 -0.37062797
|
||||
34000 200673.31 0 0.67923697 0.69170315 0 1.3709401 -0.36342828
|
||||
35000 199626.58 0 0.71214435 0.71976006 0 1.4319044 0.14520308
|
||||
36000 206317.51 0 0.73098443 0.75421862 0 1.4852031 -0.049274455
|
||||
37000 200878.1 0 0.75404669 0.77976535 0 1.533812 -0.020890744
|
||||
38000 204399.3 0 0.77905716 0.80048467 0 1.5795418 0.41855233
|
||||
39000 194807.41 0 0.80901097 0.83413027 0 1.6431412 0.061887194
|
||||
40000 197736.05 0 0.82530638 0.85259329 0 1.6778997 -0.21390647
|
||||
41000 196234.62 0 0.84399769 0.87823172 0 1.7222294 -0.069449607
|
||||
42000 202916.17 0 0.84692321 0.90553215 0 1.7524554 0.014080553
|
||||
43000 199860.49 0 0.87526292 0.92743523 0 1.8026981 0.14037107
|
||||
44000 195023.53 0 0.88991314 0.95932182 0 1.849235 -0.33208549
|
||||
45000 195356.87 0 0.91402747 0.98447964 0 1.8985071 -0.54494657
|
||||
46000 208076.31 0 0.92947502 1.0151107 0 1.9445857 -0.62483899
|
||||
47000 196855.73 0 0.94724593 1.0592619 0 2.0065078 -0.20412571
|
||||
48000 201398.95 0 0.98527482 1.0801565 0 2.0654313 -0.12234755
|
||||
49000 196643.46 0 0.9999414 1.115902 0 2.1158434 -0.014252649
|
||||
50000 195863.57 0 1.0150887 1.1141794 0 2.1292681 0.75335013
|
||||
51000 206236.4 0 1.060301 1.1367091 0 2.1970101 -0.22091636
|
||||
52000 202457.92 0 1.0661594 1.1431793 0 2.2093387 0.20119657
|
||||
53000 200285.87 0 1.0911263 1.1652507 0 2.256377 0.64981185
|
||||
54000 192647.27 0 1.1240369 1.1911613 0 2.3151982 -0.31991262
|
||||
55000 203771.16 0 1.1340164 1.2059129 0 2.3399293 -0.69282837
|
||||
56000 208567.7 0 1.1644643 1.2358751 0 2.4003394 -0.31413258
|
||||
57000 197231.77 0 1.1817765 1.2580471 0 2.4398236 0.0033296714
|
||||
58000 203466.84 0 1.2256115 1.297949 0 2.5235605 0.5978746
|
||||
59000 204071.73 0 1.2635225 1.3217703 0 2.5852928 0.080439812
|
||||
60000 204497.78 0 1.280816 1.3529427 0 2.6337587 -0.13596951
|
||||
61000 201581.53 0 1.3062524 1.3840896 0 2.6903421 0.2402402
|
||||
62000 201808.73 0 1.3158468 1.40441 0 2.7202568 -0.49735816
|
||||
63000 195485.7 0 1.3246003 1.4115588 0 2.7361591 -0.034708225
|
||||
64000 203289.43 0 1.3515942 1.4270935 0 2.7786877 -0.47288139
|
||||
65000 198360.38 0 1.3655219 1.4410644 0 2.8065863 -0.60256944
|
||||
66000 207037.46 0 1.4024293 1.447203 0 2.8496323 0.19886918
|
||||
67000 197118.92 0 1.4400601 1.4746068 0 2.914667 0.53432407
|
||||
68000 202323.37 0 1.4737023 1.4757925 0 2.9494948 -0.50123703
|
||||
69000 202679.62 0 1.4996722 1.4982292 0 2.9979014 -0.36930796
|
||||
70000 197067.46 0 1.5202881 1.5079024 0 3.0281905 0.20627183
|
||||
71000 201481.55 0 1.5292493 1.5243374 0 3.0535868 0.1664502
|
||||
72000 197565.37 0 1.5649605 1.5415579 0 3.1065184 0.1538837
|
||||
73000 198996.68 0 1.5823587 1.5687541 0 3.1511128 0.33347453
|
||||
74000 196056.36 0 1.5876207 1.5880262 0 3.1756469 0.51039007
|
||||
75000 196383.82 0 1.6238368 1.6126393 0 3.2364761 0.58359751
|
||||
76000 198050.9 0 1.6332373 1.6444992 0 3.2777365 -0.30732735
|
||||
77000 198867.35 0 1.6713412 1.6649495 0 3.3362907 -0.77019295
|
||||
78000 197809.01 0 1.7072475 1.6978103 0 3.4050578 0.062786188
|
||||
79000 198410.74 0 1.7408508 1.7258001 0 3.4666509 0.13736398
|
||||
80000 203417.5 0 1.757256 1.748869 0 3.506125 -0.020989175
|
||||
81000 195062.76 0 1.7608203 1.7602204 0 3.5210407 0.013727312
|
||||
82000 200632.52 0 1.7755665 1.7921425 0 3.567709 -0.052325244
|
||||
83000 191576.45 0 1.7959579 1.8241176 0 3.6200755 -0.20352178
|
||||
84000 200502.95 0 1.8249302 1.8426066 0 3.6675369 -0.38424345
|
||||
85000 201308.33 0 1.8545228 1.8747626 0 3.7292854 0.68968471
|
||||
86000 203836.21 0 1.8649352 1.9007982 0 3.7657335 -0.56168571
|
||||
87000 208433.36 0 1.8740109 1.9350895 0 3.8091003 0.14026532
|
||||
88000 198009.97 0 1.9059077 1.9537513 0 3.859659 0.017064255
|
||||
89000 200694.98 0 1.9312351 1.9704954 0 3.9017306 -0.19090346
|
||||
90000 200436.01 0 1.9567279 2.0115534 0 3.9682814 0.48791684
|
||||
91000 201810.57 0 1.9784333 2.0255054 0 4.0039386 -0.21421409
|
||||
92000 208365.35 0 2.0050623 2.049849 0 4.0549113 0.30465483
|
||||
93000 197213.15 0 2.017481 2.0689265 0 4.0864075 0.22419309
|
||||
94000 201428.76 0 2.0426301 2.0972527 0 4.1398828 -0.24230583
|
||||
95000 198472.77 0 2.0677457 2.111042 0 4.1787877 0.13071029
|
||||
96000 201070.22 0 2.1019013 2.1226741 0 4.2245754 0.15711614
|
||||
97000 197899.22 0 2.1019041 2.1575764 0 4.2594805 -0.18586251
|
||||
98000 196679.32 0 2.1162285 2.170629 0 4.2868575 -0.34056846
|
||||
99000 197724.36 0 2.1326497 2.1886854 0 4.321335 0.58967646
|
||||
100000 205354.32 0 2.1612412 2.2209085 0 4.3821498 0.044930992
|
||||
101000 201516.54 0 2.1954542 2.2365787 0 4.4320329 0.065295032
|
||||
102000 204976.1 0 2.2029662 2.2694803 0 4.4724465 0.11481479
|
||||
103000 194476.8 0 2.2257884 2.2955314 0 4.5213197 -0.11610949
|
||||
104000 203736.97 0 2.2268211 2.2907985 0 4.5176196 0.015187894
|
||||
105000 196842.81 0 2.2620115 2.3381811 0 4.6001926 0.32154947
|
||||
106000 202018.95 0 2.2836727 2.3560537 0 4.6397264 0.33256726
|
||||
107000 198280.55 0 2.2802504 2.3840825 0 4.6643329 0.50865119
|
||||
108000 198415.05 0 2.299079 2.402819 0 4.7018979 0.028090763
|
||||
109000 214429.79 0 2.326329 2.4071902 0 4.7335192 0.33519083
|
||||
110000 206059.4 0 2.3452258 2.431054 0 4.7762798 0.012116372
|
||||
111000 198624.14 0 2.3484889 2.4352672 0 4.7837561 -0.35529818
|
||||
112000 201890.76 0 2.3810208 2.4718016 0 4.8528225 -0.55980456
|
||||
113000 199303.2 0 2.4011131 2.5070887 0 4.9082018 0.21744731
|
||||
114000 207834.3 0 2.4341043 2.5221162 0 4.9562204 0.11696186
|
||||
115000 202055.54 0 2.4714851 2.5576047 0 5.0290898 -0.42272945
|
||||
116000 203298.5 0 2.5043401 2.5626771 0 5.0670173 0.17688701
|
||||
117000 194028.49 0 2.5273873 2.5945816 0 5.1219689 -0.70088882
|
||||
118000 201726.79 0 2.5475583 2.6532013 0 5.2007596 0.22458895
|
||||
119000 200037.8 0 2.599892 2.6701733 0 5.2700653 0.036916228
|
||||
120000 194969.52 0 2.623779 2.6921819 0 5.3159609 -0.6067451
|
||||
Loop time of 4.32114 on 1 procs for 120000 steps with 1024 atoms
|
||||
|
||||
Performance: 17136.102 tau/day, 19833.451 timesteps/s
|
||||
Performance: 23993.668 tau/day, 27770.449 timesteps/s
|
||||
100.0% CPU use with 1 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.00031301 | 0.00031301 | 0.00031301 | 0.0 | 0.01
|
||||
Comm | 0.0072472 | 0.0072472 | 0.0072472 | 0.0 | 0.12
|
||||
Output | 0.0036492 | 0.0036492 | 0.0036492 | 0.0 | 0.06
|
||||
Modify | 5.7173 | 5.7173 | 5.7173 | 0.0 | 94.50
|
||||
Other | | 0.3218 | | | 5.32
|
||||
Neigh | 0.00029646 | 0.00029646 | 0.00029646 | 0.0 | 0.01
|
||||
Comm | 0.0079174 | 0.0079174 | 0.0079174 | 0.0 | 0.18
|
||||
Output | 0.0044738 | 0.0044738 | 0.0044738 | 0.0 | 0.10
|
||||
Modify | 3.9818 | 3.9818 | 3.9818 | 0.0 | 92.15
|
||||
Other | | 0.3267 | | | 7.56
|
||||
|
||||
Nlocal: 1024.00 ave 1024 max 1024 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
@ -241,8 +237,8 @@ Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 0
|
||||
Ave neighs/atom = 0.0000000
|
||||
Neighbor list builds = 154
|
||||
Neighbor list builds = 153
|
||||
Dangerous builds = 0
|
||||
|
||||
|
||||
Total wall time: 0:00:08
|
||||
Total wall time: 0:00:06
|
||||
|
||||
@ -1,2 +0,0 @@
|
||||
The input file in3d.brownian demonstrates how to run a 3d simulation
|
||||
of ellipsoidal particles undergoing overdamped brownian motion.
|
||||
72
examples/USER/brownian/asphere/in2d.ellipsoid
Normal file
72
examples/USER/brownian/asphere/in2d.ellipsoid
Normal file
@ -0,0 +1,72 @@
|
||||
##### overdamped dynamics of non-interacting ellipsoids in 2D #####
|
||||
|
||||
|
||||
variable rng string gaussian
|
||||
variable seed string 198098
|
||||
|
||||
|
||||
variable temp string 1.0
|
||||
|
||||
variable gamma_r_1 string inf
|
||||
variable gamma_r_2 string inf
|
||||
variable gamma_r_3 string 0.1
|
||||
|
||||
variable gamma_t_1 string 5.0
|
||||
variable gamma_t_2 string 7.0
|
||||
variable gamma_t_3 string inf
|
||||
|
||||
variable params string ${rng}_${temp}_${gamma_r_1}_${gamma_r_2}_${gamma_r_3}_${gamma_t_1}_${gamma_t_2}_${gamma_t_3}
|
||||
|
||||
|
||||
log log_${params}_2d.lammps.log
|
||||
units lj
|
||||
atom_style hybrid dipole ellipsoid
|
||||
dimension 2
|
||||
newton off
|
||||
|
||||
|
||||
lattice sq 0.4
|
||||
region box block -30 30 -30 30 -0.2 0.2
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
mass * 1.0
|
||||
set type * dipole/random ${seed} 1.0
|
||||
set type * shape 3.0 1.0 1.0
|
||||
set type * quat/random ${seed}
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
|
||||
fix 1 all brownian/asphere ${temp} ${seed} rng ${rng} &
|
||||
gamma_r_eigen ${gamma_r_1} ${gamma_r_2} ${gamma_r_3} &
|
||||
gamma_t_eigen ${gamma_t_1} ${gamma_t_2} ${gamma_t_3} &
|
||||
dipole 1.0 0.0 0.0
|
||||
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*]
|
||||
|
||||
|
||||
dump 1 all custom 1000 dump_${params}_2d.lammpstrj id type &
|
||||
x y z xu yu zu mux muy muz fx fy fz
|
||||
dump_modify 1 first yes sort id
|
||||
|
||||
timestep 0.00001
|
||||
thermo 50
|
||||
|
||||
# main run
|
||||
run 30000
|
||||
@ -1,63 +0,0 @@
|
||||
# 3d overdamped brownian dynamics for ellipsoids
|
||||
# with dipole moment also being updated
|
||||
|
||||
# choose variables to obey thermal equilibrium
|
||||
# (fluctuation dissipation theorem) and no-slip
|
||||
# boundary condition
|
||||
variable rng string gaussian
|
||||
variable gamma_t equal 1.0
|
||||
variable gamma_r equal 0.33
|
||||
variable D_t equal 1.0
|
||||
variable D_r equal 3.0
|
||||
variable seed equal 1974019
|
||||
|
||||
variable params string ${rng}_${gamma_t}_${gamma_r}_${D_t}_${D_r}
|
||||
|
||||
|
||||
log log_${params}_3d.lammps.log
|
||||
units lj
|
||||
atom_style hybrid dipole sphere ellipsoid
|
||||
dimension 3
|
||||
newton off
|
||||
|
||||
|
||||
lattice sc 0.4
|
||||
region box block -4 4 -4 4 -4 4
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
mass * 1.0
|
||||
set type * dipole/random ${seed} 1.0
|
||||
set type * shape 1 1 1
|
||||
set type * quat/random ${seed}
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
fix 1 all brownian/asphere ${gamma_t} ${gamma_r} ${D_t} ${D_r} &
|
||||
${seed} rng ${rng} dipole
|
||||
|
||||
|
||||
compute press all pressure NULL virial
|
||||
|
||||
thermo_style custom step temp epair c_press
|
||||
|
||||
#equilibration
|
||||
timestep 0.0000000001
|
||||
thermo 50000
|
||||
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 120000
|
||||
73
examples/USER/brownian/asphere/in3d.ellipsoid
Normal file
73
examples/USER/brownian/asphere/in3d.ellipsoid
Normal file
@ -0,0 +1,73 @@
|
||||
##### overdamped dynamics of non-interacting ellipsoids in 3D #####
|
||||
|
||||
|
||||
variable rng string uniform
|
||||
variable seed string 198098
|
||||
|
||||
|
||||
variable temp string 1.0
|
||||
|
||||
variable gamma_r_1 string 2.0
|
||||
variable gamma_r_2 string 0.25
|
||||
variable gamma_r_3 string 0.1
|
||||
|
||||
variable gamma_t_1 string 5.0
|
||||
variable gamma_t_2 string 7.0
|
||||
variable gamma_t_3 string 9.0
|
||||
|
||||
variable params string ${rng}_${temp}_${gamma_r_1}_${gamma_r_2}_${gamma_r_3}_${gamma_t_1}_${gamma_t_2}_${gamma_t_3}
|
||||
|
||||
|
||||
log log_${params}_3d.lammps.log
|
||||
units lj
|
||||
atom_style hybrid dipole ellipsoid
|
||||
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
|
||||
set type * shape 3.0 1.0 1.0
|
||||
set type * quat/random ${seed}
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
|
||||
|
||||
fix 1 all brownian/asphere ${temp} ${seed} rng ${rng} &
|
||||
gamma_r_eigen ${gamma_r_1} ${gamma_r_2} ${gamma_r_3} &
|
||||
gamma_t_eigen ${gamma_t_1} ${gamma_t_2} ${gamma_t_3} &
|
||||
dipole 1.0 0.0 0.0
|
||||
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*]
|
||||
|
||||
|
||||
dump 1 all custom 1000 dump_${params}_3d.lammpstrj id type &
|
||||
x y z xu yu zu mux muy muz fx fy fz
|
||||
dump_modify 1 first yes sort id
|
||||
|
||||
timestep 0.00001
|
||||
thermo 50
|
||||
|
||||
# main run
|
||||
run 30000
|
||||
@ -1,141 +0,0 @@
|
||||
units lj
|
||||
atom_style hybrid dipole sphere ellipsoid
|
||||
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)
|
||||
WARNING: Peratom rmass is in multiple sub-styles - must be used consistently (src/atom_vec_hybrid.cpp:219)
|
||||
dimension 3
|
||||
newton off
|
||||
|
||||
|
||||
lattice sc 0.4
|
||||
Lattice spacing in x,y,z = 1.3572088 1.3572088 1.3572088
|
||||
region box block -4 4 -4 4 -4 4
|
||||
create_box 1 box
|
||||
Created orthogonal box = (-5.4288352 -5.4288352 -5.4288352) to (5.4288352 5.4288352 5.4288352)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 512 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 ...
|
||||
512 settings made for dipole/random
|
||||
set type * shape 1 1 1
|
||||
Setting atom values ...
|
||||
512 settings made for shape
|
||||
set type * quat/random ${seed}
|
||||
set type * quat/random 1974019
|
||||
Setting atom values ...
|
||||
512 settings made for quat/random
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
fix 1 all brownian/asphere ${gamma_t} ${gamma_r} ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 1 ${gamma_r} ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 1 0.33 ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 1 0.33 1 ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 1 0.33 1 3 ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 1 0.33 1 3 1974019 rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 1 0.33 1 3 1974019 rng gaussian dipole
|
||||
|
||||
|
||||
compute press all pressure NULL virial
|
||||
|
||||
thermo_style custom step temp epair c_press
|
||||
|
||||
#equilibration
|
||||
timestep 0.0000000001
|
||||
thermo 50000
|
||||
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) = 5.379 | 5.379 | 5.379 Mbytes
|
||||
Step Temp E_pair c_press
|
||||
0 1 0 0
|
||||
50000 1.9891104e+10 0 0
|
||||
Loop time of 5.53749 on 1 procs for 50000 steps with 512 atoms
|
||||
|
||||
Performance: 0.078 tau/day, 9029.362 timesteps/s
|
||||
100.0% CPU use with 1 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.11711 | 0.11711 | 0.11711 | 0.0 | 2.11
|
||||
Output | 2.034e-05 | 2.034e-05 | 2.034e-05 | 0.0 | 0.00
|
||||
Modify | 5.3401 | 5.3401 | 5.3401 | 0.0 | 96.44
|
||||
Other | | 0.08027 | | | 1.45
|
||||
|
||||
Nlocal: 512.000 ave 512 max 512 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 217.000 ave 217 max 217 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 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
|
||||
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) = 5.754 | 5.754 | 5.754 Mbytes
|
||||
Step Temp E_pair c_msd[1] c_msd[2] c_msd[3] c_msd[4] c_press
|
||||
0 1.9891104e+10 0 0 0 0 0 0
|
||||
10000 201972.17 0 0.19918647 0.20079752 0.20495007 0.60493407 0
|
||||
20000 197255.49 0 0.40800225 0.37910274 0.38545643 1.1725614 0
|
||||
30000 195533.67 0 0.60991554 0.5898132 0.56621596 1.7659447 0
|
||||
40000 192777.99 0 0.7761198 0.86101776 0.76531344 2.402451 0
|
||||
50000 195241.43 0 1.0053256 1.0477568 0.96681401 3.0198964 0
|
||||
60000 201887.49 0 1.2298892 1.1979933 1.1950141 3.6228965 0
|
||||
70000 200187.14 0 1.4407329 1.356743 1.3992739 4.1967498 0
|
||||
80000 202737.24 0 1.6305637 1.5663775 1.5724692 4.7694104 0
|
||||
90000 185530.51 0 1.7937597 1.7795995 1.7222809 5.2956401 0
|
||||
100000 204405.47 0 2.0149709 1.9738573 1.9423625 5.9311907 0
|
||||
110000 194892.4 0 2.1974948 2.1560014 2.1453303 6.4988264 0
|
||||
120000 198462.51 0 2.3761388 2.334739 2.287964 6.9988418 0
|
||||
Loop time of 13.0463 on 1 procs for 120000 steps with 512 atoms
|
||||
|
||||
Performance: 7947.110 tau/day, 9198.044 timesteps/s
|
||||
100.0% CPU use with 1 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.0011789 | 0.0011789 | 0.0011789 | 0.0 | 0.01
|
||||
Comm | 0.030109 | 0.030109 | 0.030109 | 0.0 | 0.23
|
||||
Output | 0.00030614 | 0.00030614 | 0.00030614 | 0.0 | 0.00
|
||||
Modify | 12.834 | 12.834 | 12.834 | 0.0 | 98.38
|
||||
Other | | 0.1803 | | | 1.38
|
||||
|
||||
Nlocal: 512.000 ave 512 max 512 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 0
|
||||
Ave neighs/atom = 0.0000000
|
||||
Neighbor list builds = 1748
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:18
|
||||
@ -1,153 +0,0 @@
|
||||
units lj
|
||||
atom_style hybrid dipole sphere ellipsoid
|
||||
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)
|
||||
WARNING: Peratom rmass is in multiple sub-styles - must be used consistently (src/atom_vec_hybrid.cpp:219)
|
||||
dimension 3
|
||||
newton off
|
||||
|
||||
|
||||
lattice sc 0.4
|
||||
Lattice spacing in x,y,z = 1.3572088 1.3572088 1.3572088
|
||||
region box block -4 4 -4 4 -4 4
|
||||
create_box 1 box
|
||||
Created orthogonal box = (-5.4288352 -5.4288352 -5.4288352) to (5.4288352 5.4288352 5.4288352)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 512 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 ...
|
||||
512 settings made for dipole/random
|
||||
set type * shape 1 1 1
|
||||
Setting atom values ...
|
||||
512 settings made for shape
|
||||
set type * quat/random ${seed}
|
||||
set type * quat/random 1974019
|
||||
Setting atom values ...
|
||||
512 settings made for quat/random
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
fix 1 all brownian/asphere ${gamma_t} ${gamma_r} ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 4 ${gamma_r} ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 4 1 ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 4 1 7 ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 4 1 7 13 ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 4 1 7 13 1974019 rng ${rng} dipole
|
||||
fix 1 all brownian/asphere 4 1 7 13 1974019 rng gaussian dipole
|
||||
|
||||
|
||||
compute press all pressure NULL virial
|
||||
|
||||
thermo_style custom step temp epair c_press
|
||||
|
||||
#equilibration
|
||||
timestep 0.0000000001
|
||||
thermo 50000
|
||||
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) = 5.379 | 5.379 | 5.379 Mbytes
|
||||
Step Temp E_pair c_press
|
||||
0 1 0 0
|
||||
50000 1.3923773e+11 0 0
|
||||
Loop time of 5.61345 on 1 procs for 50000 steps with 512 atoms
|
||||
|
||||
Performance: 0.077 tau/day, 8907.171 timesteps/s
|
||||
100.0% CPU use with 1 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.11822 | 0.11822 | 0.11822 | 0.0 | 2.11
|
||||
Output | 2.0199e-05 | 2.0199e-05 | 2.0199e-05 | 0.0 | 0.00
|
||||
Modify | 5.4135 | 5.4135 | 5.4135 | 0.0 | 96.44
|
||||
Other | | 0.0817 | | | 1.46
|
||||
|
||||
Nlocal: 512.000 ave 512 max 512 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 217.000 ave 217 max 217 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 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
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*] c_press
|
||||
|
||||
|
||||
# write trajectory and thermo in a log-scale frequency
|
||||
#dump 1 all custom 1000 dump_${params}_3d.lammpstrj id type # x y xu yu mux muy muz fx fy fz
|
||||
#dump_modify 1 first yes sort id
|
||||
|
||||
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) = 5.754 | 5.754 | 5.754 Mbytes
|
||||
Step Temp E_pair c_msd[1] c_msd[2] c_msd[3] c_msd[4] c_press
|
||||
0 1.3923773e+11 0 0 0 0 0 0
|
||||
10000 1413805.2 0 1.3943053 1.4055827 1.4346505 4.2345385 0
|
||||
20000 1380788.4 0 2.8560158 2.6537192 2.698195 8.2079299 0
|
||||
30000 1368735.7 0 4.2694087 4.1286924 3.9635117 12.361613 0
|
||||
40000 1349446 0 5.4328386 6.0271243 5.3571941 16.817157 0
|
||||
50000 1366690 0 7.0372792 7.3342977 6.7676981 21.139275 0
|
||||
60000 1413212.4 0 8.6092241 8.3859529 8.3650987 25.360276 0
|
||||
70000 1401310 0 10.085131 9.4972009 9.7949174 29.377249 0
|
||||
80000 1419160.7 0 11.413946 10.964643 11.007284 33.385873 0
|
||||
90000 1298713.5 0 12.556318 12.457196 12.055966 37.06948 0
|
||||
100000 1430838.3 0 14.104796 13.817001 13.596538 41.518335 0
|
||||
110000 1364246.8 0 15.382464 15.09201 15.017312 45.491785 0
|
||||
120000 1389237.6 0 16.632972 16.343173 16.015748 48.991892 0
|
||||
Loop time of 13.2439 on 1 procs for 120000 steps with 512 atoms
|
||||
|
||||
Performance: 7828.487 tau/day, 9060.749 timesteps/s
|
||||
100.0% CPU use with 1 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.00075178 | 0.00075178 | 0.00075178 | 0.0 | 0.01
|
||||
Comm | 0.0282 | 0.0282 | 0.0282 | 0.0 | 0.21
|
||||
Output | 0.00032692 | 0.00032692 | 0.00032692 | 0.0 | 0.00
|
||||
Modify | 13.017 | 13.017 | 13.017 | 0.0 | 98.29
|
||||
Other | | 0.1976 | | | 1.49
|
||||
|
||||
Nlocal: 512.000 ave 512 max 512 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 0
|
||||
Ave neighs/atom = 0.0000000
|
||||
Neighbor list builds = 1110
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:18
|
||||
62
examples/USER/brownian/point/in2d.point
Normal file
62
examples/USER/brownian/point/in2d.point
Normal file
@ -0,0 +1,62 @@
|
||||
##### dynamics of non-interacting point particles in 2D #####
|
||||
|
||||
variable rng string gaussian
|
||||
variable seed string 198098
|
||||
|
||||
|
||||
variable temp string 5.0
|
||||
|
||||
variable gamma_t string 1.0
|
||||
|
||||
|
||||
|
||||
variable params string ${rng}_${temp}_${gamma_t}
|
||||
|
||||
|
||||
log log_${params}_2d.lammps.log
|
||||
units lj
|
||||
atom_style atomic
|
||||
dimension 2
|
||||
newton off
|
||||
|
||||
|
||||
lattice sq 0.4
|
||||
region box block -30 30 -30 30 -0.2 0.2
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
mass * 1.0
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
fix 1 all brownian ${temp} ${seed} rng ${rng} &
|
||||
gamma_t ${gamma_t}
|
||||
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*]
|
||||
|
||||
|
||||
dump 1 all custom 1000 dump_${params}_3d.lammpstrj id type &
|
||||
x y z xu yu zu fx fy fz
|
||||
dump_modify 1 first yes sort id
|
||||
|
||||
timestep 0.00001
|
||||
thermo 50
|
||||
|
||||
# main run
|
||||
run 30000
|
||||
|
||||
61
examples/USER/brownian/point/in3d.point
Normal file
61
examples/USER/brownian/point/in3d.point
Normal file
@ -0,0 +1,61 @@
|
||||
##### overdamped dynamics of non-interacting point particles in 3D #####
|
||||
|
||||
variable rng string gaussian
|
||||
variable seed string 198098
|
||||
|
||||
|
||||
variable temp string 5.0
|
||||
|
||||
variable gamma_t string 1.0
|
||||
|
||||
|
||||
|
||||
variable params string ${rng}_${temp}_${gamma_t}
|
||||
|
||||
log log_${params}_3d.lammps.log
|
||||
units lj
|
||||
atom_style atomic
|
||||
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
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
fix 1 all brownian ${temp} ${seed} rng ${rng} &
|
||||
gamma_t ${gamma_t}
|
||||
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*]
|
||||
|
||||
|
||||
dump 1 all custom 1000 dump_${params}_3d.lammpstrj id type &
|
||||
x y z xu yu zu fx fy fz
|
||||
dump_modify 1 first yes sort id
|
||||
|
||||
timestep 0.00001
|
||||
thermo 50
|
||||
|
||||
# main run
|
||||
run 30000
|
||||
|
||||
@ -1,26 +0,0 @@
|
||||
The input file in2ddipole.brownian demonstrates how to run a
|
||||
2d simulation of particles undergoing overdamped brownian motion
|
||||
in both translational and rotational degrees of freedom. Dipole
|
||||
updating is turned on, so the package DIPOLE must be built
|
||||
for this example to work.
|
||||
|
||||
The input file in3d_virial_on.brownian demonstrates how to run
|
||||
a similar simulation but in 3d. In this case, the virial
|
||||
contribution of the brownian dynamics (the sum
|
||||
sum_i <r_i dot sqrt{2D_t}W>/(3*volume) where W is
|
||||
a random variable with mean 0 and variance dt) is
|
||||
calculated via the fix_modify command. For long
|
||||
enough times, this will be equal to rho*D_t*gamma_t
|
||||
(the ideal gas term in equilibrium systems). Note that
|
||||
no dipole updating is performed.
|
||||
|
||||
To confirm rotational diffusion is working correctly,
|
||||
run the above simulations with dump files on and
|
||||
measure \sum_i<e_i(0) dot e_i(t)> where e_i is the
|
||||
dipole vector of particle i, and one should
|
||||
find that this decays as an exponential with
|
||||
timescale 1/((d-1)*D_r).
|
||||
|
||||
Note that both of the simulations above are not long
|
||||
enough to get good statistics on e.g. ideal gas
|
||||
pressure, rotational diffusion, or translational diffusion.
|
||||
66
examples/USER/brownian/sphere/in2d.sphere
Normal file
66
examples/USER/brownian/sphere/in2d.sphere
Normal file
@ -0,0 +1,66 @@
|
||||
##### overdamped dynamics of a sphere (with dipole attached to it) in 2D #####
|
||||
|
||||
variable rng string uniform
|
||||
variable seed string 198098
|
||||
|
||||
|
||||
variable temp string 1.0
|
||||
|
||||
variable gamma_t string 5.0
|
||||
|
||||
variable gamma_r string 0.7
|
||||
|
||||
|
||||
|
||||
variable params string ${rng}_${temp}_${gamma_r}_${gamma_t}
|
||||
|
||||
|
||||
log log_${params}_2d.lammps.log
|
||||
units lj
|
||||
atom_style hybrid dipole sphere
|
||||
dimension 2
|
||||
newton off
|
||||
|
||||
|
||||
lattice sq 0.4
|
||||
region box block -30 30 -30 30 -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
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
|
||||
|
||||
fix 1 all brownian/sphere ${temp} ${seed} rng ${rng} &
|
||||
gamma_r ${gamma_r} gamma_t ${gamma_t}
|
||||
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*]
|
||||
|
||||
|
||||
dump 1 all custom 1000 dump_${params}_2d.lammpstrj id type &
|
||||
x y z xu yu zu mux muy muz fx fy fz
|
||||
dump_modify 1 first yes sort id
|
||||
|
||||
timestep 0.00001
|
||||
thermo 50
|
||||
|
||||
# main run
|
||||
run 30000
|
||||
@ -1,59 +0,0 @@
|
||||
# 2d overdamped brownian dynamics of a sphere
|
||||
# with dipole also being updated
|
||||
|
||||
variable rng string gaussian
|
||||
variable gamma_t equal 4.0
|
||||
variable gamma_r equal 1.0
|
||||
variable D_t equal 7.0
|
||||
variable D_r equal 13.0
|
||||
variable seed equal 1974019
|
||||
|
||||
|
||||
variable params string ${rng}_${gamma_t}_${gamma_r}_${D_t}_${D_r}
|
||||
|
||||
|
||||
log log_${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
|
||||
|
||||
pair_style none
|
||||
|
||||
fix 1 all brownian/sphere ${gamma_t} ${gamma_r} ${D_t} &
|
||||
${D_r} ${seed} rng ${rng} 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 1000
|
||||
|
||||
# main run
|
||||
run 120000
|
||||
|
||||
|
||||
65
examples/USER/brownian/sphere/in3d.sphere
Normal file
65
examples/USER/brownian/sphere/in3d.sphere
Normal file
@ -0,0 +1,65 @@
|
||||
##### overdamped dynamics of a sphere (with dipole attached to it) in 3D#####
|
||||
|
||||
variable rng string uniform
|
||||
variable seed string 198098
|
||||
|
||||
|
||||
variable temp string 1.0
|
||||
|
||||
variable gamma_t string 5.0
|
||||
|
||||
variable gamma_r string 0.7
|
||||
|
||||
|
||||
|
||||
variable params string ${rng}_${temp}_${gamma_r}_${gamma_t}
|
||||
|
||||
|
||||
log log_${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
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
|
||||
fix 1 all brownian/sphere ${temp} ${seed} rng ${rng} &
|
||||
gamma_r ${gamma_r} gamma_t ${gamma_t}
|
||||
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*]
|
||||
|
||||
|
||||
dump 1 all custom 1000 dump_${params}_3d.lammpstrj id type &
|
||||
x y z xu yu zu mux muy muz fx fy fz
|
||||
dump_modify 1 first yes sort id
|
||||
|
||||
timestep 0.00001
|
||||
thermo 50
|
||||
|
||||
# main run
|
||||
run 30000
|
||||
@ -1,59 +0,0 @@
|
||||
# 3d overdamped brownian dynamics of sphere, with
|
||||
# virial contribution (i.e. ideal gas pressure)
|
||||
# included
|
||||
|
||||
|
||||
variable rng string uniform
|
||||
variable gamma_t equal 10.0
|
||||
variable gamma_r equal 1.0
|
||||
variable D_t equal 5.0
|
||||
variable D_r equal 15.0
|
||||
variable seed equal 553910
|
||||
|
||||
|
||||
variable params string ${rng}_${gamma_t}_${gamma_r}_${D_t}_${D_r}
|
||||
|
||||
|
||||
log log_${params}_3d.lammps.log
|
||||
units lj
|
||||
atom_style sphere
|
||||
dimension 3
|
||||
|
||||
|
||||
lattice sc 0.4
|
||||
region box block -8 8 -8 8 -8 8
|
||||
create_box 1 box
|
||||
create_atoms 1 box
|
||||
#mass * 1.0
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
fix 1 all brownian/sphere ${gamma_t} ${gamma_r} ${D_t} &
|
||||
${D_r} ${seed} rng ${rng}
|
||||
fix_modify 1 virial yes
|
||||
|
||||
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 1000
|
||||
|
||||
# main run
|
||||
run 120000
|
||||
|
||||
|
||||
@ -1,260 +0,0 @@
|
||||
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)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 1024 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 ...
|
||||
1024 settings made for dipole/random
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
|
||||
#compute d all property/atom mux muy muz
|
||||
|
||||
fix 1 all brownian/sphere ${gamma_t} ${gamma_r} ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/sphere 4 ${gamma_r} ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/sphere 4 1 ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/sphere 4 1 7 ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/sphere 4 1 7 13 ${seed} rng ${rng} dipole
|
||||
fix 1 all brownian/sphere 4 1 7 13 1974019 rng ${rng} dipole
|
||||
fix 1 all brownian/sphere 4 1 7 13 1974019 rng gaussian 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
|
||||
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.289 | 4.289 | 4.289 Mbytes
|
||||
Step Temp E_pair c_press
|
||||
0 1 0 0
|
||||
50000 7.3671759e+10 0 0
|
||||
Loop time of 11.2532 on 1 procs for 50000 steps with 1024 atoms
|
||||
|
||||
Performance: 0.038 tau/day, 4443.198 timesteps/s
|
||||
100.0% CPU use with 1 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.031456 | 0.031456 | 0.031456 | 0.0 | 0.28
|
||||
Output | 2.0958e-05 | 2.0958e-05 | 2.0958e-05 | 0.0 | 0.00
|
||||
Modify | 11.067 | 11.067 | 11.067 | 0.0 | 98.35
|
||||
Other | | 0.1546 | | | 1.37
|
||||
|
||||
Nlocal: 1024.00 ave 1024 max 1024 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 65.0000 ave 65 max 65 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 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
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*] c_press
|
||||
|
||||
|
||||
# write trajectory and thermo in a log-scale frequency
|
||||
# uncomment next three lines for dump output
|
||||
dump 1 all custom 2000 dump_${params}_2d.lammpstrj id type x y xu yu mux muy muz fx fy fz
|
||||
dump 1 all custom 2000 dump_gaussian_4_1_7_13_2d.lammpstrj id type x y xu yu mux muy muz fx fy fz
|
||||
dump_modify 1 first yes sort id
|
||||
|
||||
timestep 0.00001
|
||||
thermo 1000
|
||||
|
||||
# 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) = 6.113 | 6.113 | 6.113 Mbytes
|
||||
Step Temp E_pair c_msd[1] c_msd[2] c_msd[3] c_msd[4] c_press
|
||||
0 7.3671759e+10 0 0 0 0 0 0
|
||||
1000 710744.42 0 0.1332856 0.13577642 0 0.26906203 0
|
||||
2000 705025.5 0 0.27078701 0.27355073 0 0.54433774 0
|
||||
3000 750845.29 0 0.41036556 0.39710766 0 0.80747321 0
|
||||
4000 752542.47 0 0.53097803 0.53654063 0 1.0675187 0
|
||||
5000 731812.26 0 0.66927129 0.65872533 0 1.3279966 0
|
||||
6000 716119.24 0 0.79097615 0.76825364 0 1.5592298 0
|
||||
7000 726294.28 0 0.93322099 0.91220222 0 1.8454232 0
|
||||
8000 673771.54 0 1.0789689 1.031843 0 2.1108119 0
|
||||
9000 757173.95 0 1.2328252 1.2024778 0 2.435303 0
|
||||
10000 693711.95 0 1.3972184 1.3045575 0 2.7017759 0
|
||||
11000 774897.65 0 1.584036 1.4168797 0 3.0009157 0
|
||||
12000 764910.09 0 1.7344079 1.5751172 0 3.3095251 0
|
||||
13000 752859.5 0 1.8538323 1.7336107 0 3.5874429 0
|
||||
14000 746333.26 0 2.0392353 1.9027299 0 3.9419652 0
|
||||
15000 735704.24 0 2.1714465 2.000952 0 4.1723985 0
|
||||
16000 722768.01 0 2.258542 2.1084289 0 4.3669709 0
|
||||
17000 741408.63 0 2.4194247 2.2469958 0 4.6664206 0
|
||||
18000 713609.29 0 2.5366244 2.3325203 0 4.8691447 0
|
||||
19000 770066.65 0 2.6662758 2.4726928 0 5.1389686 0
|
||||
20000 739828.53 0 2.792873 2.6188232 0 5.4116962 0
|
||||
21000 719478.96 0 2.9109491 2.7872008 0 5.6981499 0
|
||||
22000 737448.77 0 3.0162218 2.9213733 0 5.9375951 0
|
||||
23000 714695.3 0 3.1309486 3.0325142 0 6.1634628 0
|
||||
24000 712386.78 0 3.31348 3.1323675 0 6.4458475 0
|
||||
25000 747298.11 0 3.5148965 3.2978465 0 6.812743 0
|
||||
26000 737240.19 0 3.7117818 3.4245004 0 7.1362822 0
|
||||
27000 720878.36 0 3.8342853 3.5681829 0 7.4024682 0
|
||||
28000 776824.55 0 3.9904931 3.6853691 0 7.6758622 0
|
||||
29000 714696.4 0 4.0733919 3.8807741 0 7.954166 0
|
||||
30000 725827.3 0 4.2213715 4.0552488 0 8.2766204 0
|
||||
31000 737978.89 0 4.2918759 4.1795647 0 8.4714407 0
|
||||
32000 766620.93 0 4.4620184 4.2765194 0 8.7385378 0
|
||||
33000 714739.71 0 4.5942749 4.426881 0 9.0211559 0
|
||||
34000 694821.06 0 4.686348 4.5057433 0 9.1920912 0
|
||||
35000 721948.38 0 4.7526371 4.5689873 0 9.3216245 0
|
||||
36000 713621.16 0 4.8645036 4.6355302 0 9.5000339 0
|
||||
37000 704079.85 0 4.9652601 4.8235856 0 9.7888457 0
|
||||
38000 706914.29 0 5.1045998 4.9585169 0 10.063117 0
|
||||
39000 736940.33 0 5.2172626 5.1168713 0 10.334134 0
|
||||
40000 738302.73 0 5.3720228 5.3207578 0 10.692781 0
|
||||
41000 746518.85 0 5.5376524 5.5192708 0 11.056923 0
|
||||
42000 719970.82 0 5.7610696 5.6335312 0 11.394601 0
|
||||
43000 736875.87 0 5.8769052 5.7504356 0 11.627341 0
|
||||
44000 765218.38 0 6.0308866 5.9650683 0 11.995955 0
|
||||
45000 739382.69 0 6.1577029 6.0736304 0 12.231333 0
|
||||
46000 730967.87 0 6.3846223 6.2090389 0 12.593661 0
|
||||
47000 758881.95 0 6.5604053 6.3883315 0 12.948737 0
|
||||
48000 751135.53 0 6.6160248 6.5199017 0 13.135926 0
|
||||
49000 721058.97 0 6.6124519 6.6938253 0 13.306277 0
|
||||
50000 686287.15 0 6.8527601 6.8518895 0 13.70465 0
|
||||
51000 735175.7 0 6.9928329 6.8663476 0 13.85918 0
|
||||
52000 760085.67 0 7.1009497 7.0345784 0 14.135528 0
|
||||
53000 741211.89 0 7.2287282 7.2318431 0 14.460571 0
|
||||
54000 723668.97 0 7.3732634 7.3847301 0 14.757994 0
|
||||
55000 751846.88 0 7.5566634 7.4927365 0 15.0494 0
|
||||
56000 762648.27 0 7.7651818 7.5379796 0 15.303161 0
|
||||
57000 710175.2 0 7.9051989 7.6644542 0 15.569653 0
|
||||
58000 740068.94 0 8.0296854 7.8307471 0 15.860433 0
|
||||
59000 728119.43 0 8.1562433 8.039333 0 16.195576 0
|
||||
60000 728768.48 0 8.316192 8.1643843 0 16.480576 0
|
||||
61000 707538.83 0 8.3727565 8.3895027 0 16.762259 0
|
||||
62000 713931.41 0 8.4657927 8.5517529 0 17.017546 0
|
||||
63000 742604.95 0 8.6550684 8.7189851 0 17.374053 0
|
||||
64000 738981.54 0 8.7331005 8.9208616 0 17.653962 0
|
||||
65000 701685.87 0 8.898571 8.977933 0 17.876504 0
|
||||
66000 716151.39 0 9.0229699 9.1292403 0 18.15221 0
|
||||
67000 732641.29 0 9.151611 9.3249179 0 18.476529 0
|
||||
68000 740870.85 0 9.233452 9.4667607 0 18.700213 0
|
||||
69000 733159.17 0 9.1743719 9.5648134 0 18.739185 0
|
||||
70000 693562.75 0 9.3414089 9.7425241 0 19.083933 0
|
||||
71000 748694.53 0 9.5036102 9.8686571 0 19.372267 0
|
||||
72000 721046.02 0 9.7116931 10.0195 0 19.731193 0
|
||||
73000 736631.66 0 9.7095662 10.205306 0 19.914872 0
|
||||
74000 751264.08 0 9.9413234 10.327844 0 20.269167 0
|
||||
75000 729223.27 0 10.211903 10.412516 0 20.624419 0
|
||||
76000 747811 0 10.332266 10.544251 0 20.876517 0
|
||||
77000 717505.01 0 10.521195 10.737774 0 21.258969 0
|
||||
78000 712288.89 0 10.712514 10.874932 0 21.587446 0
|
||||
79000 728912.34 0 10.755227 10.968257 0 21.723484 0
|
||||
80000 743505.67 0 11.026063 11.070234 0 22.096298 0
|
||||
81000 732218.14 0 11.321417 11.31773 0 22.639147 0
|
||||
82000 777186.61 0 11.480074 11.401465 0 22.881539 0
|
||||
83000 734805.95 0 11.7524 11.62552 0 23.37792 0
|
||||
84000 753703.38 0 11.863318 11.833925 0 23.697243 0
|
||||
85000 755730.75 0 12.068564 11.937143 0 24.005707 0
|
||||
86000 725021.19 0 12.258195 12.034586 0 24.292781 0
|
||||
87000 731844.67 0 12.341844 12.331792 0 24.673636 0
|
||||
88000 707368.15 0 12.431452 12.547314 0 24.978766 0
|
||||
89000 784756.28 0 12.405699 12.686021 0 25.09172 0
|
||||
90000 760278.97 0 12.542669 12.851925 0 25.394593 0
|
||||
91000 753765.97 0 12.703534 13.08092 0 25.784454 0
|
||||
92000 705869.22 0 12.971838 13.210092 0 26.18193 0
|
||||
93000 741699.59 0 13.110745 13.383115 0 26.49386 0
|
||||
94000 717372.41 0 13.252109 13.466388 0 26.718497 0
|
||||
95000 721820.49 0 13.373146 13.660605 0 27.033751 0
|
||||
96000 777409.97 0 13.596822 13.886961 0 27.483782 0
|
||||
97000 741480.91 0 13.752545 14.204624 0 27.957169 0
|
||||
98000 725755.72 0 13.918302 14.219292 0 28.137594 0
|
||||
99000 729710.24 0 14.065417 14.213142 0 28.278559 0
|
||||
100000 700366.74 0 14.278054 14.324135 0 28.602189 0
|
||||
101000 722737.18 0 14.398464 14.46904 0 28.867505 0
|
||||
102000 758930.85 0 14.408098 14.637964 0 29.046061 0
|
||||
103000 725233.66 0 14.526482 14.719688 0 29.24617 0
|
||||
104000 752416.52 0 14.641393 14.952087 0 29.59348 0
|
||||
105000 769047.5 0 14.780788 15.182945 0 29.963733 0
|
||||
106000 734849.98 0 14.840982 15.316112 0 30.157094 0
|
||||
107000 720210.74 0 15.010029 15.372792 0 30.382821 0
|
||||
108000 741216.55 0 15.113143 15.488575 0 30.601718 0
|
||||
109000 714158.43 0 15.057499 15.592865 0 30.650364 0
|
||||
110000 743262.25 0 15.381281 15.798368 0 31.179649 0
|
||||
111000 728836.92 0 15.488226 16.034645 0 31.522871 0
|
||||
112000 753490.67 0 15.679979 16.350819 0 32.030799 0
|
||||
113000 730699.5 0 15.82813 16.680136 0 32.508266 0
|
||||
114000 705526.22 0 16.099699 16.920095 0 33.019794 0
|
||||
115000 752408.36 0 16.393466 17.053935 0 33.447401 0
|
||||
116000 713697.16 0 16.634939 17.35349 0 33.98843 0
|
||||
117000 690041.5 0 16.727379 17.457936 0 34.185315 0
|
||||
118000 745594.37 0 16.838659 17.585392 0 34.424052 0
|
||||
119000 714487.42 0 17.064214 17.743478 0 34.807692 0
|
||||
120000 716557.14 0 17.105558 17.845925 0 34.951483 0
|
||||
Loop time of 28.1507 on 1 procs for 120000 steps with 1024 atoms
|
||||
|
||||
Performance: 3683.037 tau/day, 4262.774 timesteps/s
|
||||
99.9% CPU use with 1 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.0016237 | 0.0016237 | 0.0016237 | 0.0 | 0.01
|
||||
Comm | 0.027064 | 0.027064 | 0.027064 | 0.0 | 0.10
|
||||
Output | 0.1712 | 0.1712 | 0.1712 | 0.0 | 0.61
|
||||
Modify | 27.545 | 27.545 | 27.545 | 0.0 | 97.85
|
||||
Other | | 0.4062 | | | 1.44
|
||||
|
||||
Nlocal: 1024.00 ave 1024 max 1024 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 0
|
||||
Ave neighs/atom = 0.0000000
|
||||
Neighbor list builds = 1049
|
||||
Dangerous builds = 0
|
||||
|
||||
|
||||
Total wall time: 0:00:39
|
||||
@ -1,249 +0,0 @@
|
||||
units lj
|
||||
atom_style sphere
|
||||
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)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 4096 atoms
|
||||
create_atoms CPU = 0.001 seconds
|
||||
#mass * 1.0
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
fix 1 all bd/sphere ${gamma_t} ${gamma_r} ${D_t} ${D_r} ${seed} rng ${rng}
|
||||
fix 1 all bd/sphere 10 ${gamma_r} ${D_t} ${D_r} ${seed} rng ${rng}
|
||||
fix 1 all bd/sphere 10 1 ${D_t} ${D_r} ${seed} rng ${rng}
|
||||
fix 1 all bd/sphere 10 1 5 ${D_r} ${seed} rng ${rng}
|
||||
fix 1 all bd/sphere 10 1 5 15 ${seed} rng ${rng}
|
||||
fix 1 all bd/sphere 10 1 5 15 553910 rng ${rng}
|
||||
fix 1 all bd/sphere 10 1 5 15 553910 rng uniform
|
||||
fix_modify 1 virial yes
|
||||
|
||||
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) = 3.581 | 3.581 | 3.581 Mbytes
|
||||
Step Temp E_pair c_press
|
||||
0 1 0 -32553.271
|
||||
50000 5.2351457e+10 0 -15026.42
|
||||
Loop time of 16.0183 on 1 procs for 50000 steps with 4096 atoms
|
||||
|
||||
Performance: 0.027 tau/day, 3121.426 timesteps/s
|
||||
100.0% CPU use with 1 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.074984 | 0.074984 | 0.074984 | 0.0 | 0.47
|
||||
Output | 3.0041e-05 | 3.0041e-05 | 3.0041e-05 | 0.0 | 0.00
|
||||
Modify | 15.366 | 15.366 | 15.366 | 0.0 | 95.93
|
||||
Other | | 0.5772 | | | 3.60
|
||||
|
||||
Nlocal: 4096.00 ave 4096 max 4096 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 817.000 ave 817 max 817 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 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
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*] c_press
|
||||
|
||||
|
||||
# write trajectory and thermo in a log-scale frequency
|
||||
# uncomment the next three lines for dump file
|
||||
#dump 1 all custom 10000 dump_${params}_3d.lammpstrj id type # x y xu yu fx fy fz
|
||||
#dump_modify 1 first yes sort id
|
||||
|
||||
timestep 0.00001
|
||||
thermo 1000
|
||||
|
||||
# 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) = 3.956 | 3.956 | 3.956 Mbytes
|
||||
Step Temp E_pair c_msd[1] c_msd[2] c_msd[3] c_msd[4] c_press
|
||||
0 5.2351457e+10 0 0 0 0 0 -233.83012
|
||||
1000 522590.88 0 0.10259269 0.10487519 0.10391005 0.31137793 -13.730626
|
||||
2000 521776.08 0 0.20662097 0.2118237 0.20371673 0.62216139 -141.90848
|
||||
3000 529408.9 0 0.30711557 0.3105565 0.30293728 0.92060935 28.335187
|
||||
4000 518069.42 0 0.39918073 0.40643835 0.40168024 1.2072993 -12.262269
|
||||
5000 520487.75 0 0.50222348 0.50854461 0.50304351 1.5138116 303.56775
|
||||
6000 527970.5 0 0.60443861 0.60363711 0.59884932 1.806925 218.55799
|
||||
7000 519424.68 0 0.69529071 0.70025532 0.70347351 2.0990195 -363.51723
|
||||
8000 521925.72 0 0.80545039 0.81520807 0.79930708 2.4199655 510.47732
|
||||
9000 518397.28 0 0.91236827 0.9050426 0.89431985 2.7117307 537.74151
|
||||
10000 525696.47 0 0.99494325 0.99433832 0.98574884 2.9750304 411.56896
|
||||
11000 524505.88 0 1.0863677 1.1232636 1.0810609 3.2906921 32.652569
|
||||
12000 525350.08 0 1.1816136 1.2412674 1.189996 3.612877 45.304061
|
||||
13000 524598.12 0 1.2841566 1.3381915 1.2906249 3.912973 -376.50739
|
||||
14000 525438.09 0 1.3964076 1.4420032 1.3915611 4.2299719 -31.273339
|
||||
15000 523531.08 0 1.4829127 1.5436611 1.4892074 4.5157811 -357.26366
|
||||
16000 532366.43 0 1.5871965 1.6432286 1.6020146 4.8324397 435.10667
|
||||
17000 523629.21 0 1.6968632 1.7574539 1.6967316 5.1510487 280.59816
|
||||
18000 521432.94 0 1.7973983 1.8462477 1.8087416 5.4523876 134.69384
|
||||
19000 530778.78 0 1.888653 1.9328611 1.9268155 5.7483296 -64.346115
|
||||
20000 525623.21 0 1.9718697 2.0048196 2.0068177 5.983507 73.884028
|
||||
21000 524862.26 0 2.0414226 2.0967449 2.0878261 6.2259936 163.95173
|
||||
22000 529123.04 0 2.143889 2.1850282 2.1972491 6.5261662 234.87793
|
||||
23000 519629.76 0 2.2297952 2.2661611 2.2805197 6.776476 57.65702
|
||||
24000 523746.18 0 2.3293654 2.3658544 2.3842082 7.079428 46.291089
|
||||
25000 515851.07 0 2.4177426 2.4559274 2.5029961 7.3766661 -16.353546
|
||||
26000 515592.59 0 2.5352044 2.5438918 2.5657776 7.6448738 -65.259248
|
||||
27000 527764.96 0 2.6498176 2.6423592 2.6715377 7.9637146 324.48367
|
||||
28000 527177.18 0 2.7405087 2.7431537 2.7702133 8.2538757 -47.906223
|
||||
29000 521079.66 0 2.8294642 2.8410515 2.8546167 8.5251323 -135.36173
|
||||
30000 526092.8 0 2.961898 2.9361476 2.9319586 8.8300042 206.28635
|
||||
31000 529595.93 0 3.0638234 3.0128995 3.0434232 9.1201461 -35.946596
|
||||
32000 521173.83 0 3.159344 3.1397814 3.1455147 9.4446401 328.73193
|
||||
33000 522989.77 0 3.2567127 3.2348491 3.2444068 9.7359686 362.33435
|
||||
34000 522930.45 0 3.3489514 3.3433962 3.3391822 10.03153 178.73773
|
||||
35000 530664.64 0 3.4445399 3.4843743 3.4205275 10.349442 146.02898
|
||||
36000 521443.26 0 3.520475 3.6171115 3.5179112 10.655498 244.07172
|
||||
37000 529669.54 0 3.6076942 3.7233869 3.6169718 10.948053 79.46077
|
||||
38000 527521.93 0 3.7093699 3.8275855 3.7451165 11.282072 -91.131074
|
||||
39000 523237.65 0 3.8088585 3.9559041 3.8748751 11.639638 -70.265271
|
||||
40000 522798.75 0 3.89319 4.0420406 3.9196106 11.854841 91.631855
|
||||
41000 524228.77 0 3.9863379 4.1602581 4.0393136 12.18591 113.26515
|
||||
42000 524546.08 0 4.0969366 4.2738786 4.1218115 12.492627 313.83957
|
||||
43000 515948.22 0 4.2161418 4.3700009 4.2178962 12.804039 92.993253
|
||||
44000 527319.42 0 4.3267431 4.4688062 4.3274029 13.122952 -118.91153
|
||||
45000 526222.74 0 4.4289127 4.5531203 4.4604075 13.442441 -474.98113
|
||||
46000 524617.54 0 4.5422968 4.6316235 4.5445894 13.71851 -303.23458
|
||||
47000 524067.07 0 4.6284944 4.6989745 4.5884153 13.915884 -343.41107
|
||||
48000 535659.25 0 4.7221068 4.7505302 4.7212887 14.193926 149.75872
|
||||
49000 527425.65 0 4.8098179 4.8761354 4.8362998 14.522253 -272.18936
|
||||
50000 516588.98 0 4.9077783 4.9950242 4.9826783 14.885481 -52.925733
|
||||
51000 522291.86 0 5.0022325 5.0659091 5.0525647 15.120706 -52.985883
|
||||
52000 530056.07 0 5.123751 5.1470734 5.1603805 15.431205 -231.29334
|
||||
53000 525624.96 0 5.2412036 5.2543541 5.2717587 15.767316 294.19678
|
||||
54000 524407.27 0 5.3530927 5.359043 5.342451 16.054587 -55.282671
|
||||
55000 519989.37 0 5.4458173 5.4781678 5.413865 16.33785 -135.28414
|
||||
56000 524987.68 0 5.5282178 5.5825979 5.5127172 16.623533 84.654044
|
||||
57000 525610.66 0 5.6202713 5.6925409 5.5964466 16.909259 407.16526
|
||||
58000 523097.93 0 5.7390671 5.7830376 5.6921201 17.214225 -271.05243
|
||||
59000 525357.74 0 5.8037507 5.8654514 5.7920744 17.461276 -213.07815
|
||||
60000 522185.96 0 5.9067925 5.9909357 5.903525 17.801253 330.09637
|
||||
61000 527694.33 0 6.0576096 6.0907943 6.0202838 18.168688 310.22884
|
||||
62000 522936.57 0 6.1422565 6.1888677 6.124212 18.455336 -24.591697
|
||||
63000 526642.55 0 6.2261061 6.2832608 6.2023277 18.711695 -130.32823
|
||||
64000 514826.49 0 6.3032614 6.3628539 6.3184362 18.984552 195.43036
|
||||
65000 529338.23 0 6.4152502 6.4482512 6.4241507 19.287652 -77.817059
|
||||
66000 527425.07 0 6.5122194 6.5052947 6.5518104 19.569324 8.0143425
|
||||
67000 525258.85 0 6.6292334 6.6417881 6.6438938 19.914915 152.1905
|
||||
68000 527067.89 0 6.7569918 6.6951219 6.7357134 20.187827 -337.54654
|
||||
69000 522395.33 0 6.8401085 6.7836633 6.8243746 20.448146 34.755329
|
||||
70000 522200.34 0 6.9549662 6.8773008 6.9537894 20.786056 352.30283
|
||||
71000 527550.11 0 7.0441399 6.9837133 7.0860791 21.113932 171.42143
|
||||
72000 529564.22 0 7.1744229 7.0581688 7.1507041 21.383296 269.04569
|
||||
73000 524175.53 0 7.248828 7.1693905 7.2974868 21.715705 -84.993601
|
||||
74000 528604.1 0 7.3358903 7.2983204 7.390523 22.024734 -6.8495991
|
||||
75000 519136.15 0 7.445248 7.4224403 7.4690469 22.336735 250.23521
|
||||
76000 518470.58 0 7.5747768 7.5354904 7.5612252 22.671492 274.53956
|
||||
77000 518459.51 0 7.6660567 7.6557812 7.7020797 23.023918 430.07722
|
||||
78000 516646.23 0 7.7596088 7.7130901 7.7896654 23.262364 536.21298
|
||||
79000 517511.57 0 7.9021784 7.82098 7.8840028 23.607161 193.80468
|
||||
80000 520654.84 0 7.991107 7.8932454 8.0049184 23.889271 -6.0066355
|
||||
81000 522331.74 0 8.0616391 7.9990606 8.0954141 24.156114 47.271954
|
||||
82000 521540.3 0 8.1268306 8.0763829 8.1861513 24.389365 36.265762
|
||||
83000 529154.07 0 8.2537975 8.1854268 8.2745822 24.713807 113.94154
|
||||
84000 525977.29 0 8.2999431 8.26591 8.2970498 24.862903 -72.444836
|
||||
85000 520505.41 0 8.3872349 8.3551374 8.4251458 25.167518 -346.6703
|
||||
86000 516334.67 0 8.4947842 8.4504822 8.5238751 25.469142 -3.8332672
|
||||
87000 524745.34 0 8.5405071 8.5348571 8.6344387 25.709803 -149.22382
|
||||
88000 521384.87 0 8.5798155 8.6260785 8.7133565 25.919251 22.009599
|
||||
89000 519680.71 0 8.6657392 8.6857605 8.8319304 26.18343 -180.94487
|
||||
90000 522482.4 0 8.7310324 8.7844751 8.9636862 26.479194 -102.50361
|
||||
91000 527450.69 0 8.8250902 8.8445728 9.0456076 26.715271 325.68024
|
||||
92000 519060.77 0 8.9054775 8.9412744 9.1259074 26.972659 130.73028
|
||||
93000 518170.14 0 9.0378587 9.0315875 9.2346184 27.304065 -111.68498
|
||||
94000 531639.53 0 9.1268874 9.0952299 9.3437063 27.565824 -26.423328
|
||||
95000 519776.1 0 9.2332576 9.1903006 9.4516001 27.875158 -113.46356
|
||||
96000 527082.72 0 9.2993234 9.2533436 9.534409 28.087076 -46.720611
|
||||
97000 526965.16 0 9.3485677 9.3197496 9.6631942 28.331512 159.18386
|
||||
98000 529790.04 0 9.4759113 9.5177629 9.783846 28.77752 343.32872
|
||||
99000 520964.82 0 9.6488022 9.6573355 9.882432 29.18857 -257.85576
|
||||
100000 520863.58 0 9.7452168 9.7659565 9.9878877 29.499061 -108.52324
|
||||
101000 526760.86 0 9.8073751 9.8508213 10.127532 29.785728 120.22249
|
||||
102000 519249.89 0 9.9315855 9.9767409 10.221605 30.129931 -176.50647
|
||||
103000 525003.9 0 10.023982 10.062451 10.315002 30.401435 422.1401
|
||||
104000 519112.73 0 10.086117 10.136879 10.415222 30.638218 -147.66505
|
||||
105000 517898.24 0 10.180438 10.240942 10.525346 30.946725 158.63054
|
||||
106000 528046.35 0 10.258163 10.41411 10.601663 31.273936 -128.04669
|
||||
107000 527105.8 0 10.364015 10.485036 10.795177 31.644228 183.82165
|
||||
108000 522024.28 0 10.450786 10.544065 10.902961 31.897812 -39.692553
|
||||
109000 519497.83 0 10.556206 10.61633 11.046222 32.218758 173.75988
|
||||
110000 521070.8 0 10.64856 10.745382 11.071387 32.465329 -128.45389
|
||||
111000 525657.64 0 10.830485 10.852954 11.188295 32.871734 -214.60249
|
||||
112000 519426.33 0 10.987769 10.97575 11.26012 33.223639 292.35901
|
||||
113000 526472.45 0 11.029941 11.086376 11.418772 33.535089 189.69245
|
||||
114000 520070.28 0 11.107229 11.183295 11.454757 33.745281 -40.433571
|
||||
115000 525812.59 0 11.153992 11.305378 11.537521 33.99689 -106.38733
|
||||
116000 524464.26 0 11.256779 11.413082 11.633295 34.303156 -159.59643
|
||||
117000 519838.94 0 11.344413 11.480104 11.706366 34.530883 -2.0346135
|
||||
118000 524075.83 0 11.441416 11.597533 11.783056 34.822005 -350.05313
|
||||
119000 533816.71 0 11.52766 11.681986 11.917713 35.127359 -521.58975
|
||||
120000 524509.31 0 11.63865 11.792667 12.080379 35.511696 273.09647
|
||||
Loop time of 39.6359 on 1 procs for 120000 steps with 4096 atoms
|
||||
|
||||
Performance: 2615.809 tau/day, 3027.557 timesteps/s
|
||||
100.0% CPU use with 1 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.0059817 | 0.0059817 | 0.0059817 | 0.0 | 0.02
|
||||
Comm | 0.039473 | 0.039473 | 0.039473 | 0.0 | 0.10
|
||||
Output | 0.0063653 | 0.0063653 | 0.0063653 | 0.0 | 0.02
|
||||
Modify | 38.085 | 38.085 | 38.085 | 0.0 | 96.09
|
||||
Other | | 1.499 | | | 3.78
|
||||
|
||||
Nlocal: 4096.00 ave 4096 max 4096 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 0
|
||||
Ave neighs/atom = 0.0000000
|
||||
Neighbor list builds = 992
|
||||
Dangerous builds = 0
|
||||
|
||||
|
||||
Total wall time: 0:00:55
|
||||
@ -1,14 +1,13 @@
|
||||
# 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 gamma_t string 1.0
|
||||
variable gamma_r string 1.0
|
||||
variable temp string 1.0
|
||||
variable seed equal 1974019
|
||||
variable fp equal 4.0
|
||||
variable fp string 4.0
|
||||
|
||||
variable params string ${gamma_t}_${gamma_r}_${D_t}_${D_r}_${fp}
|
||||
variable params string ${temp}_${gamma_t}_${gamma_r}_${fp}
|
||||
|
||||
|
||||
log log_WCA_${params}_2d.lammps.log
|
||||
@ -39,7 +38,7 @@ 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 ${temp} ${seed} gamma_t ${gamma_t} gamma_r ${gamma_r}
|
||||
# self-propulsion force along the dipole direction
|
||||
fix activity all propel/self dipole ${fp}
|
||||
fix 2 all enforce2d
|
||||
@ -1,13 +1,12 @@
|
||||
# 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 gamma_t string 3.0
|
||||
variable gamma_r string 1.0
|
||||
variable temp string 1.0
|
||||
variable seed equal 1974019
|
||||
variable fp equal 4.0
|
||||
variable fp string 4.0
|
||||
|
||||
variable params string ${gamma_t}_${gamma_r}_${D_t}_${D_r}_${fp}
|
||||
variable params string ${temp}_${gamma_t}_${gamma_r}_${fp}
|
||||
|
||||
|
||||
log log_ideal_${params}_3d.lammps.log
|
||||
@ -28,7 +27,7 @@ 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 ${temp} ${seed} gamma_t ${gamma_t} gamma_r ${gamma_r}
|
||||
# self-propulsion force along the dipole direction
|
||||
fix activity all propel/self dipole ${fp}
|
||||
|
||||
@ -47,7 +46,8 @@ 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
|
||||
# MSD = (2*d*kb*T/gamma_t + 2*fp**2*gamma_r/(kb*T*gamma_t**2*(d-1)))*t
|
||||
# + 2*fp**2*gamma_r**2/(gamma_t**2*(d-1)**2*(kb*T)**2)*(e^(-(d-1)*t*kb*T/gamma_r)-1)
|
||||
#
|
||||
# with d being simulation dimension
|
||||
compute msd all msd
|
||||
@ -1,65 +0,0 @@
|
||||
# 3d overdamped brownian dynamics
|
||||
|
||||
variable gamma_t equal 1.0
|
||||
variable D_t equal 1.0
|
||||
variable seed equal 1974019
|
||||
|
||||
variable params string ${gamma_t}_${D_t}
|
||||
|
||||
|
||||
log log_${params}.lammps.log
|
||||
units lj
|
||||
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
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
# simple overdamped brownian dynamics time-stepper
|
||||
fix step all brownian ${gamma_t} ${D_t} ${seed}
|
||||
|
||||
# turn on the virial contribution from the noise
|
||||
# (this is the ideal gas pressure, but it is really noisy
|
||||
# for small systems)
|
||||
fix_modify step virial yes
|
||||
|
||||
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
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*] c_press
|
||||
|
||||
|
||||
timestep 0.00001
|
||||
thermo 1000
|
||||
|
||||
# main run
|
||||
run 120000
|
||||
|
||||
|
||||
@ -1,246 +0,0 @@
|
||||
units lj
|
||||
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)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
create_atoms 1 box
|
||||
Created 4096 atoms
|
||||
create_atoms CPU = 0.002 seconds
|
||||
mass * 1.0
|
||||
velocity all create 1.0 1 loop geom
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
# simple overdamped brownian dynamics time-stepper
|
||||
fix step all brownian ${gamma_t} ${D_t} ${seed}
|
||||
fix step all brownian 1 ${D_t} ${seed}
|
||||
fix step all brownian 1 1 ${seed}
|
||||
fix step all brownian 1 1 1974019
|
||||
|
||||
# turn on the virial contribution from the noise
|
||||
# (this is the ideal gas pressure, but it is really noisy
|
||||
# for small systems)
|
||||
fix_modify step virial yes
|
||||
|
||||
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) = 2.319 | 2.319 | 2.319 Mbytes
|
||||
Step Temp E_pair c_press
|
||||
0 1 0 1500.0667
|
||||
50000 2.0204192e+10 0 809.28898
|
||||
Loop time of 8.36695 on 1 procs for 50000 steps with 4096 atoms
|
||||
|
||||
Performance: 0.052 tau/day, 5975.895 timesteps/s
|
||||
100.0% CPU use with 1 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.081773 | 0.081773 | 0.081773 | 0.0 | 0.98
|
||||
Output | 3.0396e-05 | 3.0396e-05 | 3.0396e-05 | 0.0 | 0.00
|
||||
Modify | 7.8039 | 7.8039 | 7.8039 | 0.0 | 93.27
|
||||
Other | | 0.4812 | | | 5.75
|
||||
|
||||
Nlocal: 4096.00 ave 4096 max 4096 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 817.000 ave 817 max 817 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 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
|
||||
|
||||
|
||||
#initialisation for the main run
|
||||
|
||||
# MSD
|
||||
compute msd all msd
|
||||
|
||||
|
||||
thermo_style custom step temp epair c_msd[*] c_press
|
||||
|
||||
|
||||
timestep 0.00001
|
||||
thermo 1000
|
||||
|
||||
# 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) = 2.694 | 2.694 | 2.694 Mbytes
|
||||
Step Temp E_pair c_msd[1] c_msd[2] c_msd[3] c_msd[4] c_press
|
||||
0 2.0204192e+10 0 0 0 0 0 3.053092
|
||||
1000 199589.87 0 0.020515986 0.019445659 0.020084014 0.060045658 -10.855191
|
||||
2000 202255.37 0 0.039423387 0.039839097 0.041243768 0.12050625 -7.8484884
|
||||
3000 197932.37 0 0.058057959 0.060232381 0.061205106 0.17949545 14.441358
|
||||
4000 201354.42 0 0.08020915 0.080140903 0.081351793 0.24170185 -2.7088264
|
||||
5000 201599.34 0 0.10065125 0.099423954 0.10004367 0.30011888 10.000367
|
||||
6000 198355.3 0 0.12039302 0.12228165 0.12136204 0.36403671 -7.7331104
|
||||
7000 201842.08 0 0.13901495 0.14401324 0.14070032 0.42372851 1.049759
|
||||
8000 200224.96 0 0.16063188 0.16389028 0.16409878 0.48862093 14.748306
|
||||
9000 198468.62 0 0.18050666 0.18555949 0.18441359 0.55047975 -7.2751841
|
||||
10000 197958.21 0 0.20229316 0.20438608 0.20471694 0.61139618 -8.8169492
|
||||
11000 200852.8 0 0.22244961 0.22388152 0.22417791 0.67050904 -0.56323783
|
||||
12000 200135.84 0 0.24070342 0.23964066 0.24829459 0.72863867 -20.246433
|
||||
13000 199160.33 0 0.26021951 0.2581362 0.27210493 0.79046065 6.2400743
|
||||
14000 200016.58 0 0.28005738 0.27429833 0.29107673 0.84543244 9.8390268
|
||||
15000 200805.25 0 0.30110658 0.29756813 0.31002936 0.90870407 -9.7474978
|
||||
16000 201221.83 0 0.31982855 0.31758368 0.33140494 0.96881716 -11.937757
|
||||
17000 199690.18 0 0.33896659 0.33908504 0.35143307 1.0294847 -11.627541
|
||||
18000 197165.97 0 0.35589324 0.3605585 0.37062404 1.0870758 17.045279
|
||||
19000 202921.66 0 0.36968243 0.3825335 0.38579438 1.1380103 -3.2942099
|
||||
20000 200415.26 0 0.39189833 0.39811842 0.40650662 1.1965234 -0.55679809
|
||||
21000 198882.79 0 0.41394816 0.42495341 0.42660986 1.2655114 0.26572202
|
||||
22000 199352.13 0 0.43200564 0.44555233 0.44842807 1.325986 5.9839474
|
||||
23000 200565.11 0 0.45248704 0.46736174 0.4633773 1.3832261 0.51141727
|
||||
24000 197632.73 0 0.47648878 0.49070755 0.48561349 1.4528098 11.301781
|
||||
25000 203284.74 0 0.50103025 0.5117464 0.50162451 1.5144012 4.0966379
|
||||
26000 200458.42 0 0.52111051 0.53245594 0.52193205 1.5754985 -9.2090189
|
||||
27000 197553.02 0 0.54097922 0.553215 0.5440757 1.6382699 -6.0708365
|
||||
28000 200874.32 0 0.55707209 0.57015852 0.56334199 1.6905726 -17.349073
|
||||
29000 199494.43 0 0.57734934 0.59276209 0.58014344 1.7502549 -7.0407471
|
||||
30000 199898.09 0 0.60161157 0.6156313 0.60516078 1.8224036 -1.3770813
|
||||
31000 200867.7 0 0.62383361 0.64008503 0.6272796 1.8911982 -8.3200079
|
||||
32000 198811.59 0 0.64361889 0.66089102 0.64987533 1.9543852 -0.21133799
|
||||
33000 198224.12 0 0.66695348 0.67669734 0.67749939 2.0211502 0.87075722
|
||||
34000 202553.2 0 0.69470773 0.69930784 0.69600623 2.0900218 1.7428524
|
||||
35000 199227.86 0 0.71973625 0.72255845 0.71991572 2.1622104 2.4381317
|
||||
36000 200719.05 0 0.73772312 0.74086704 0.74438632 2.2229765 -8.8184183
|
||||
37000 198054.45 0 0.75901548 0.76107029 0.76677002 2.2868558 9.0699905
|
||||
38000 198496.92 0 0.77679013 0.78067416 0.78055585 2.3380201 -8.8374756
|
||||
39000 199224.19 0 0.79689891 0.8074183 0.79624042 2.4005576 -2.3099574
|
||||
40000 199514.11 0 0.81493522 0.8271518 0.81729253 2.4593796 8.5411105
|
||||
41000 199926.55 0 0.82918834 0.84950171 0.83531381 2.5140039 -6.8276624
|
||||
42000 201659.21 0 0.84991161 0.87416096 0.86612093 2.5901935 -11.006396
|
||||
43000 201502.51 0 0.87967597 0.89570393 0.88920119 2.6645811 9.4305203
|
||||
44000 194956.27 0 0.9022655 0.91604833 0.90539218 2.723706 1.7636771
|
||||
45000 201853.02 0 0.92302957 0.93661741 0.93225051 2.7918975 9.6483674
|
||||
46000 200572.92 0 0.94062662 0.94933155 0.94953376 2.8394919 1.0882497
|
||||
47000 202008.49 0 0.95397524 0.97655157 0.96534559 2.8958724 -7.9450141
|
||||
48000 199748.89 0 0.97441115 0.99687233 0.98429144 2.9555749 5.0525854
|
||||
49000 203008.7 0 0.99757948 1.0095536 0.99836015 3.0054932 28.410878
|
||||
50000 198810.18 0 1.0182591 1.0287254 1.0110652 3.0580497 22.596989
|
||||
51000 201716 0 1.0414747 1.0448379 1.0361279 3.1224404 -0.3397634
|
||||
52000 200682.63 0 1.0640391 1.0597767 1.0596255 3.1834413 8.9163814
|
||||
53000 201068.23 0 1.0804112 1.0740541 1.077586 3.2320513 8.9631815
|
||||
54000 203379.33 0 1.0965663 1.0832317 1.0981473 3.2779453 -7.8020084
|
||||
55000 197117.81 0 1.1145489 1.1008769 1.1259188 3.3413446 13.723633
|
||||
56000 201100.37 0 1.1420502 1.1311309 1.1425839 3.415765 10.39045
|
||||
57000 199911.61 0 1.1641357 1.1461183 1.1598876 3.4701416 -1.5384143
|
||||
58000 202180.4 0 1.1787793 1.1703422 1.1794284 3.5285498 4.9552552
|
||||
59000 202359.35 0 1.2017068 1.1846051 1.1912556 3.5775675 14.774737
|
||||
60000 196182.35 0 1.2297664 1.2087508 1.2116928 3.65021 12.484104
|
||||
61000 201858.45 0 1.2505094 1.2255583 1.2330327 3.7091004 3.6962199
|
||||
62000 198313.02 0 1.2756755 1.2437434 1.2514738 3.7708927 -6.5600779
|
||||
63000 198423.1 0 1.2998594 1.262369 1.2656961 3.8279245 3.1919497
|
||||
64000 202601.56 0 1.3177969 1.2812961 1.2841642 3.8832572 7.1735962
|
||||
65000 201270.73 0 1.3353533 1.3006911 1.306668 3.9427124 3.7612957
|
||||
66000 199480.32 0 1.3502663 1.3165416 1.3300249 3.9968327 -8.0484056
|
||||
67000 202829.47 0 1.3628924 1.3369328 1.339246 4.0390712 -2.962791
|
||||
68000 200269.53 0 1.3901352 1.3600404 1.3528283 4.1030039 29.395886
|
||||
69000 201514.31 0 1.4135333 1.3834796 1.3719116 4.1689245 9.5358653
|
||||
70000 198898.14 0 1.4323413 1.4056025 1.4015088 4.2394526 13.713608
|
||||
71000 198446.83 0 1.4424061 1.4229225 1.4231698 4.2884984 9.5266069
|
||||
72000 199550.94 0 1.4730822 1.4438093 1.4436864 4.360578 5.4060926
|
||||
73000 201536.57 0 1.4860349 1.4557531 1.4673417 4.4091297 2.7527606
|
||||
74000 201688.68 0 1.5078921 1.4770318 1.4889958 4.4739197 -2.5507167
|
||||
75000 203168.01 0 1.5264946 1.4942757 1.5095341 4.5303045 13.955087
|
||||
76000 198782.86 0 1.5363851 1.5145529 1.5306446 4.5815826 1.6748995
|
||||
77000 199306.04 0 1.55336 1.5308045 1.5501462 4.6343107 1.4463705
|
||||
78000 199420.52 0 1.5679804 1.5445387 1.5765555 4.6890745 -15.926362
|
||||
79000 198711.33 0 1.5818761 1.5651424 1.5966885 4.743707 -23.662716
|
||||
80000 202734.35 0 1.6045976 1.5878467 1.6155179 4.8079623 1.4177771
|
||||
81000 199953.12 0 1.6266918 1.6125491 1.6370893 4.8763302 14.40188
|
||||
82000 201254.27 0 1.6418366 1.6362867 1.6497081 4.9278313 -4.3503
|
||||
83000 200442.35 0 1.6671772 1.6544269 1.6638838 4.9854878 -5.5569751
|
||||
84000 201442.09 0 1.6772029 1.6819435 1.6824211 5.0415676 -17.634517
|
||||
85000 202012.75 0 1.7026782 1.7059915 1.7079961 5.1166659 8.8311485
|
||||
86000 198613.4 0 1.725679 1.7251262 1.7235885 5.1743936 10.070509
|
||||
87000 198650.79 0 1.748013 1.7447674 1.7397963 5.2325767 7.3326989
|
||||
88000 199753.54 0 1.7793864 1.7677848 1.7548586 5.3020298 -3.8460869
|
||||
89000 200851.56 0 1.7992856 1.7843281 1.7671516 5.3507653 10.274664
|
||||
90000 202029.24 0 1.8178051 1.8010651 1.7886511 5.4075213 1.1705818
|
||||
91000 200963.98 0 1.8425785 1.8123748 1.8079693 5.4629226 -8.0390883
|
||||
92000 200443.4 0 1.8743616 1.8355736 1.8340515 5.5439867 -12.186363
|
||||
93000 197974.1 0 1.8911764 1.8440063 1.8442582 5.579441 10.189145
|
||||
94000 201285.41 0 1.9040394 1.8567044 1.8663386 5.6270824 7.8537148
|
||||
95000 198472.81 0 1.9268236 1.8638624 1.8855767 5.6762627 11.556616
|
||||
96000 198171.93 0 1.9378011 1.8811168 1.9024245 5.7213424 -7.3493903
|
||||
97000 200055.5 0 1.9539002 1.9031647 1.9221125 5.7791774 -7.2823252
|
||||
98000 200350.77 0 1.973424 1.9255707 1.9393139 5.8383087 2.3526328
|
||||
99000 198923.17 0 1.9946733 1.9416292 1.9616759 5.8979785 -8.6362233
|
||||
100000 196205.31 0 2.015688 1.967164 1.9801696 5.9630216 1.5261152
|
||||
101000 198842.45 0 2.0402634 1.9858628 1.9939889 6.0201151 6.8070808
|
||||
102000 199060.56 0 2.0583018 2.0040652 2.0225396 6.0849067 5.4626963
|
||||
103000 204892.64 0 2.0788003 2.0245826 2.0405395 6.1439224 -19.988675
|
||||
104000 198709.07 0 2.0990457 2.0519007 2.0571079 6.2080544 -21.365135
|
||||
105000 198916.99 0 2.1089408 2.0758832 2.0899796 6.2748037 4.3696238
|
||||
106000 202191.68 0 2.1172909 2.0923523 2.1208274 6.3304706 8.2072292
|
||||
107000 202428.89 0 2.1387532 2.114944 2.1418816 6.3955788 4.1099611
|
||||
108000 197690.67 0 2.1620575 2.136726 2.1703555 6.469139 5.7183695
|
||||
109000 200098.73 0 2.1814376 2.1464455 2.1828177 6.5107008 7.4366333
|
||||
110000 197901.18 0 2.1955124 2.1764141 2.1994286 6.5713551 -6.4288954
|
||||
111000 199478.54 0 2.2167884 2.1900638 2.2140739 6.6209261 22.379137
|
||||
112000 198391.65 0 2.249996 2.2100316 2.2309406 6.6909682 -20.040892
|
||||
113000 200542.42 0 2.2634106 2.2313768 2.2610988 6.7558862 0.2953844
|
||||
114000 202117.15 0 2.28441 2.2517036 2.2787302 6.8148438 16.75177
|
||||
115000 200004.06 0 2.2957226 2.2730837 2.2901883 6.8589947 -4.3125612
|
||||
116000 200648.11 0 2.3184059 2.2934521 2.3257075 6.9375656 5.7210624
|
||||
117000 198600.57 0 2.3413891 2.3102468 2.3511234 7.0027593 -2.9987639
|
||||
118000 199817.34 0 2.3732051 2.3347741 2.3601752 7.0681544 -3.3658539
|
||||
119000 200556.96 0 2.3873448 2.3595646 2.3774937 7.1244031 20.860601
|
||||
120000 200997.81 0 2.4097258 2.3736031 2.3871081 7.170437 -5.3623487
|
||||
Loop time of 20.5037 on 1 procs for 120000 steps with 4096 atoms
|
||||
|
||||
Performance: 5056.640 tau/day, 5852.593 timesteps/s
|
||||
99.9% CPU use with 1 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.0012425 | 0.0012425 | 0.0012425 | 0.0 | 0.01
|
||||
Comm | 0.015074 | 0.015074 | 0.015074 | 0.0 | 0.07
|
||||
Output | 0.006156 | 0.006156 | 0.006156 | 0.0 | 0.03
|
||||
Modify | 19.243 | 19.243 | 19.243 | 0.0 | 93.85
|
||||
Other | | 1.238 | | | 6.04
|
||||
|
||||
Nlocal: 4096.00 ave 4096 max 4096 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Nghost: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
Neighs: 0.00000 ave 0 max 0 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 0
|
||||
Ave neighs/atom = 0.0000000
|
||||
Neighbor list builds = 205
|
||||
Dangerous builds = 0
|
||||
|
||||
|
||||
Total wall time: 0:00:28
|
||||
@ -38,69 +38,24 @@ using namespace FixConst;
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBrownian::FixBrownian(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
FixBrownianBase(lmp, narg, arg)
|
||||
{
|
||||
virial_flag = 1;
|
||||
|
||||
time_integrate = 1;
|
||||
|
||||
if (narg != 6 && narg != 8)
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
|
||||
gamma_t = utils::numeric(FLERR,arg[3],false,lmp);
|
||||
if (gamma_t <= 0.0)
|
||||
error->all(FLERR,"Fix brownian viscous drag "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
diff_t = utils::numeric(FLERR,arg[4],false,lmp);
|
||||
if (diff_t <= 0.0)
|
||||
error->all(FLERR,"Fix brownian diffusion "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
seed = utils::inumeric(FLERR,arg[5],false,lmp);
|
||||
if (seed <= 0) error->all(FLERR,"Fix brownian seed must be > 0.");
|
||||
|
||||
noise_flag = 1;
|
||||
gaussian_noise_flag = 0;
|
||||
|
||||
if (narg == 8) {
|
||||
|
||||
if (strcmp(arg[6],"rng") == 0) {
|
||||
if (strcmp(arg[7],"uniform") == 0) {
|
||||
noise_flag = 1;
|
||||
} else if (strcmp(arg[7],"gaussian") == 0) {
|
||||
noise_flag = 1;
|
||||
gaussian_noise_flag = 1;
|
||||
} else if (strcmp(arg[7],"none") == 0) {
|
||||
noise_flag = 0;
|
||||
} else {
|
||||
if (dipole_flag || gamma_t_eigen_flag || gamma_r_eigen_flag || gamma_r_flag) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
} else {
|
||||
|
||||
if (!gamma_t_flag) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
}
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
random = new RanMars(lmp,seed + comm->me);
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBrownian::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= POST_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBrownian::~FixBrownian()
|
||||
{
|
||||
delete random;
|
||||
}
|
||||
|
||||
|
||||
@ -110,108 +65,98 @@ FixBrownian::~FixBrownian()
|
||||
void FixBrownian::init()
|
||||
{
|
||||
|
||||
g1 = force->ftm2v/gamma_t;
|
||||
if (noise_flag == 0) {
|
||||
g2 = 0;
|
||||
rng_func = &RanMars::zero_rng;
|
||||
} else if (gaussian_noise_flag == 1) {
|
||||
g2 = gamma_t*sqrt(2 * diff_t)/force->ftm2v;
|
||||
rng_func = &RanMars::gaussian;
|
||||
} else {
|
||||
g2 = gamma_t*sqrt( 24 * diff_t)/force->ftm2v;
|
||||
rng_func = &RanMars::uniform_middle;
|
||||
}
|
||||
FixBrownianBase::init();
|
||||
|
||||
dt = update->dt;
|
||||
sqrtdt = sqrt(dt);
|
||||
}
|
||||
g1 /= gamma_t;
|
||||
g2 *= sqrt(gamma_t);
|
||||
|
||||
void FixBrownian::setup(int vflag)
|
||||
{
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBrownian::initial_integrate(int /* vflag */)
|
||||
void FixBrownian::initial_integrate(int /*vflag */)
|
||||
{
|
||||
|
||||
if (domain->dimension == 2) {
|
||||
if (!noise_flag) {
|
||||
initial_integrate_templated<0,0,1>();
|
||||
} else if (gaussian_noise_flag) {
|
||||
initial_integrate_templated<0,1,1>();
|
||||
} else {
|
||||
initial_integrate_templated<1,0,1>();
|
||||
}
|
||||
} else {
|
||||
if (!noise_flag) {
|
||||
initial_integrate_templated<0,0,0>();
|
||||
} else if (gaussian_noise_flag) {
|
||||
initial_integrate_templated<0,1,0>();
|
||||
} else {
|
||||
initial_integrate_templated<1,0,0>();
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D >
|
||||
void FixBrownian::initial_integrate_templated()
|
||||
{
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
double dx,dy,dz;
|
||||
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
double dx,dy,dz;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
if (Tp_2D) {
|
||||
dz = 0;
|
||||
if (Tp_UNIFORM) {
|
||||
dx = dt * (g1 * f[i][0] + g2 * (random->uniform()-0.5));
|
||||
dy = dt * (g1 * f[i][1] + g2 * (random->uniform()-0.5));
|
||||
|
||||
} else if (Tp_GAUSS) {
|
||||
dx = dt * (g1 * f[i][0] + g2 * random->gaussian());
|
||||
dy = dt * (g1 * f[i][1] + g2 * random->gaussian());
|
||||
} else {
|
||||
dx = dt * g1 * f[i][0];
|
||||
dy = dt * g1 * f[i][1];
|
||||
}
|
||||
} else {
|
||||
if (Tp_UNIFORM) {
|
||||
dx = dt * (g1 * f[i][0] + g2 * (random->uniform()-0.5));
|
||||
dy = dt * (g1 * f[i][1] + g2 * (random->uniform()-0.5));
|
||||
dz = dt * (g1 * f[i][2] + g2 * (random->uniform()-0.5));
|
||||
} else if (Tp_GAUSS) {
|
||||
dx = dt * (g1 * f[i][0] + g2 * random->gaussian());
|
||||
dy = dt * (g1 * f[i][1] + g2 * random->gaussian());
|
||||
dz = dt * (g1 * f[i][2] + g2 * random->gaussian());
|
||||
} else {
|
||||
dx = dt * g1 * f[i][0];
|
||||
dy = dt * g1 * f[i][1];
|
||||
dz = dt * g1 * f[i][2];
|
||||
}
|
||||
}
|
||||
|
||||
x[i][0] += dx;
|
||||
v[i][0] = dx/dt;
|
||||
|
||||
dy = dt * g1 * f[i][1];
|
||||
|
||||
x[i][1] += dy;
|
||||
v[i][1] = dy/dt;
|
||||
|
||||
dz = dt * g1 * f[i][2];
|
||||
|
||||
x[i][2] += dz;
|
||||
v[i][2] = dz/dt;
|
||||
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
apply random force, stolen from MISC/fix_efield.cpp
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixBrownian::post_force(int vflag)
|
||||
{
|
||||
double **f = atom->f;
|
||||
double **x = atom->x;
|
||||
int *mask = atom->mask;
|
||||
imageint *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
// virial setup
|
||||
|
||||
if (vflag) v_setup(vflag);
|
||||
else evflag = 0;
|
||||
|
||||
double fx,fy,fz;
|
||||
double v[6];
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
fx = g2 * (random->*rng_func)()/sqrtdt;
|
||||
fy = g2 * (random->*rng_func)()/sqrtdt;
|
||||
fz = g2 * (random->*rng_func)()/sqrtdt;
|
||||
f[i][0] += fx;
|
||||
f[i][1] += fy;
|
||||
f[i][2] += fz;
|
||||
|
||||
if (evflag) {
|
||||
v[0] = fx*x[i][0];
|
||||
v[1] = fy*x[i][1];
|
||||
v[2] = fz*x[i][2];
|
||||
v[3] = fx*x[i][1];
|
||||
v[4] = fx*x[i][2];
|
||||
v[5] = fy*x[i][2];
|
||||
v_tally(i, v);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void FixBrownian::reset_dt()
|
||||
{
|
||||
|
||||
dt = update->dt;
|
||||
sqrtdt = sqrt(dt);
|
||||
}
|
||||
|
||||
@ -20,38 +20,22 @@ FixStyle(brownian,FixBrownian)
|
||||
#ifndef LMP_FIX_BROWNIAN_H
|
||||
#define LMP_FIX_BROWNIAN_H
|
||||
|
||||
#include "fix.h"
|
||||
#include "fix_brownian_base.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixBrownian : public Fix {
|
||||
class FixBrownian : public FixBrownianBase {
|
||||
public:
|
||||
FixBrownian(class LAMMPS *, int, char **);
|
||||
virtual ~FixBrownian();
|
||||
void init();
|
||||
void initial_integrate(int);
|
||||
void setup(int);
|
||||
void post_force(int);
|
||||
int setmask();
|
||||
void reset_dt();
|
||||
|
||||
private:
|
||||
int seed; // RNG seed
|
||||
double dt, sqrtdt; // time step interval and its sqrt
|
||||
template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D >
|
||||
void initial_integrate_templated();
|
||||
|
||||
|
||||
double gamma_t; // translational damping param
|
||||
double diff_t; // translational diffusion coeff
|
||||
|
||||
double g1,g2; // prefactors in time stepping
|
||||
int noise_flag; // 0/1 for noise off/on
|
||||
int gaussian_noise_flag; // 0/1 for uniform/gaussian noise
|
||||
|
||||
protected:
|
||||
class RanMars *random;
|
||||
typedef double (RanMars::*rng_member)();
|
||||
rng_member rng_func; // placeholder for RNG function
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
@ -36,104 +36,34 @@
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
#define SMALL 1e-14
|
||||
|
||||
enum{NODIPOLE,DIPOLE};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBrownianAsphere::FixBrownianAsphere(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
FixBrownianBase(lmp, narg, arg)
|
||||
{
|
||||
virial_flag = 1;
|
||||
|
||||
time_integrate = 1;
|
||||
|
||||
dipole_flag = NODIPOLE;
|
||||
|
||||
|
||||
if (narg > 11 || narg < 8 )
|
||||
error->all(FLERR,"Illegal fix brownian/asphere command.");
|
||||
|
||||
if (!atom->sphere_flag)
|
||||
error->all(FLERR,"Fix brownian/asphere requires atom style sphere");
|
||||
|
||||
gamma_t = utils::numeric(FLERR,arg[3],false,lmp);
|
||||
if (gamma_t <= 0.0)
|
||||
error->all(FLERR,"Fix brownian/asphere translational viscous drag "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
gamma_r = utils::numeric(FLERR,arg[4],false,lmp);
|
||||
if (gamma_t <= 0.0)
|
||||
error->all(FLERR,"Fix brownian/asphere rotational viscous drag "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
|
||||
diff_t = utils::numeric(FLERR,arg[5],false,lmp);
|
||||
if (diff_t <= 0.0)
|
||||
error->all(FLERR,"Fix brownian/asphere translational diffusion "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
diff_r = utils::numeric(FLERR,arg[6],false,lmp);
|
||||
if (diff_r <= 0.0)
|
||||
error->all(FLERR,"Fix brownian/asphere rotational diffusion "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
seed = utils::inumeric(FLERR,arg[7],false,lmp);
|
||||
if (seed <= 0) error->all(FLERR,"Fix brownian/asphere seed must be > 0.");
|
||||
|
||||
noise_flag = 1;
|
||||
gaussian_noise_flag = 0;
|
||||
|
||||
int iarg = 8;
|
||||
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"rng") == 0) {
|
||||
if (narg == iarg + 1) {
|
||||
error->all(FLERR,"Illegal fix/brownian/asphere command.");
|
||||
}
|
||||
if (strcmp(arg[iarg + 1],"uniform") == 0) {
|
||||
noise_flag = 1;
|
||||
} else if (strcmp(arg[iarg + 1],"gaussian") == 0) {
|
||||
noise_flag = 1;
|
||||
gaussian_noise_flag = 1;
|
||||
} else if (strcmp(arg[iarg + 1],"none") == 0) {
|
||||
noise_flag = 0;
|
||||
} else {
|
||||
error->all(FLERR,"Illegal fix/brownian/asphere command.");
|
||||
}
|
||||
iarg = iarg + 2;
|
||||
} else if (strcmp(arg[iarg],"dipole") == 0) {
|
||||
dipole_flag = DIPOLE;
|
||||
iarg = iarg + 1;
|
||||
} else {
|
||||
error->all(FLERR,"Illegal fix/brownian/asphere command.");
|
||||
}
|
||||
if (!gamma_t_eigen_flag || !gamma_r_eigen_flag) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
|
||||
if (dipole_flag == DIPOLE && !atom->mu_flag)
|
||||
if (gamma_t_flag || gamma_r_flag) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
|
||||
|
||||
|
||||
if (dipole_flag && !atom->mu_flag)
|
||||
error->all(FLERR,"Fix brownian/asphere dipole requires atom attribute mu");
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
random = new RanMars(lmp,seed + comm->me);
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBrownianAsphere::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= POST_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBrownianAsphere::~FixBrownianAsphere()
|
||||
{
|
||||
delete random;
|
||||
}
|
||||
|
||||
|
||||
@ -162,9 +92,8 @@ void FixBrownianAsphere::init()
|
||||
error->one(FLERR,"Fix brownian/asphere requires extended particles");
|
||||
|
||||
|
||||
if (dipole_flag == DIPOLE) {
|
||||
if (dipole_flag) {
|
||||
|
||||
double f_act[3] = { 1.0, 0.0, 0.0 };
|
||||
double f_rot[3];
|
||||
double *quat;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
@ -178,7 +107,7 @@ void FixBrownianAsphere::init()
|
||||
if (mask[i] & groupbit) {
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
MathExtra::quat_to_mat( quat, Q );
|
||||
MathExtra::matvec( Q, f_act, f_rot );
|
||||
MathExtra::matvec( Q, dipole_body, f_rot );
|
||||
|
||||
mu[i][0] = f_rot[0];
|
||||
mu[i][1] = f_rot[1];
|
||||
@ -188,205 +117,243 @@ void FixBrownianAsphere::init()
|
||||
}
|
||||
}
|
||||
|
||||
g1 = force->ftm2v/gamma_t;
|
||||
g3 = force->ftm2v/gamma_r;
|
||||
if (noise_flag == 0) {
|
||||
g2 = 0;
|
||||
g4 = 0;
|
||||
rng_func = &RanMars::zero_rng;
|
||||
} else if (gaussian_noise_flag == 1) {
|
||||
g2 = gamma_t*sqrt(2 * diff_t)/force->ftm2v;
|
||||
g4 = gamma_r*sqrt(2 * diff_r)/force->ftm2v;
|
||||
rng_func = &RanMars::gaussian;
|
||||
} else {
|
||||
g2 = gamma_t*sqrt( 24 * diff_t)/force->ftm2v;
|
||||
g4 = gamma_r*sqrt( 24 * diff_r )/force->ftm2v;
|
||||
rng_func = &RanMars::uniform_middle;
|
||||
}
|
||||
FixBrownianBase::init();
|
||||
|
||||
dt = update->dt;
|
||||
sqrtdt = sqrt(dt);
|
||||
}
|
||||
|
||||
void FixBrownianAsphere::setup(int vflag)
|
||||
{
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBrownianAsphere::initial_integrate(int /*vflag */)
|
||||
{
|
||||
|
||||
void FixBrownianAsphere::initial_integrate(int /* vflag */)
|
||||
if (domain->dimension == 2) {
|
||||
if (dipole_flag) {
|
||||
if (!noise_flag) {
|
||||
initial_integrate_templated<0,0,1,1>();
|
||||
} else if (gaussian_noise_flag) {
|
||||
initial_integrate_templated<0,1,1,1>();
|
||||
} else {
|
||||
initial_integrate_templated<1,0,1,1>();
|
||||
}
|
||||
} else {
|
||||
if (!noise_flag) {
|
||||
initial_integrate_templated<0,0,0,1>();
|
||||
} else if (gaussian_noise_flag) {
|
||||
initial_integrate_templated<0,1,0,1>();
|
||||
} else {
|
||||
initial_integrate_templated<1,0,0,1>();
|
||||
}
|
||||
}
|
||||
} else {
|
||||
if (dipole_flag) {
|
||||
if (!noise_flag) {
|
||||
initial_integrate_templated<0,0,1,0>();
|
||||
} else if (gaussian_noise_flag) {
|
||||
initial_integrate_templated<0,1,1,0>();
|
||||
} else {
|
||||
initial_integrate_templated<1,0,1,0>();
|
||||
}
|
||||
} else {
|
||||
if (!noise_flag) {
|
||||
initial_integrate_templated<0,0,0,0>();
|
||||
} else if (gaussian_noise_flag) {
|
||||
initial_integrate_templated<0,1,0,0>();
|
||||
} else {
|
||||
initial_integrate_templated<1,0,0,0>();
|
||||
}
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D >
|
||||
void FixBrownianAsphere::initial_integrate_templated()
|
||||
{
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double **omega = atom->omega;
|
||||
double **torque = atom->torque;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
|
||||
|
||||
int d3rot; // whether to compute angular momentum in xy plane
|
||||
|
||||
if (domain->dimension==2) {
|
||||
d3rot = 0;
|
||||
} else {
|
||||
d3rot = 1;
|
||||
}
|
||||
double **mu = atom->mu;
|
||||
double **torque = atom->torque;
|
||||
|
||||
if (dipole_flag == DIPOLE) {
|
||||
|
||||
// if dipole is being tracked, then update it along with
|
||||
// quaternions accordingly along with angular velocity
|
||||
|
||||
double wq[4];
|
||||
double qw[4];
|
||||
|
||||
double *quat;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
|
||||
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
|
||||
|
||||
// project dipole along x axis of quat
|
||||
double f_act[3] = { 1.0, 0.0, 0.0 };
|
||||
double f_rot[3];
|
||||
|
||||
double Q[3][3];
|
||||
|
||||
double **mu = atom->mu;
|
||||
double rotationmatrix_transpose[3][3];
|
||||
|
||||
double tmp[3];
|
||||
double dv[3];
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
update_x_and_omega(x[i],v[i],omega[i],f[i],torque[i],d3rot);
|
||||
// update orientation first
|
||||
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
|
||||
MathExtra::vecquat(omega[i],quat,wq);
|
||||
MathExtra::quat_to_mat_trans( quat, rotationmatrix_transpose );
|
||||
|
||||
quat[0] = quat[0] + 0.5*dt*wq[0];
|
||||
quat[1] = quat[1] + 0.5*dt*wq[1];
|
||||
quat[2] = quat[2] + 0.5*dt*wq[2];
|
||||
quat[3] = quat[3] + 0.5*dt*wq[3];
|
||||
// tmp holds angular velocity in body frame
|
||||
MathExtra::matvec(rotationmatrix_transpose,torque[i],tmp);
|
||||
|
||||
if (Tp_2D) {
|
||||
tmp[0] = tmp[1] = 0.0;
|
||||
if (Tp_UNIFORM) {
|
||||
tmp[2] = g1*tmp[2]*gamma_r_inv[2] + gamma_r_invsqrt[2]*(random->uniform()-0.5)*g2;
|
||||
|
||||
} else if (Tp_GAUSS) {
|
||||
tmp[2] = g1*tmp[2]*gamma_r_inv[2] + gamma_r_invsqrt[2]*random->gaussian()*g2;
|
||||
|
||||
} else {
|
||||
tmp[2] = g1*tmp[2]*gamma_r_inv[2];
|
||||
|
||||
}
|
||||
} else {
|
||||
if (Tp_UNIFORM) {
|
||||
|
||||
tmp[0] = g1*tmp[0]*gamma_r_inv[0] + gamma_r_invsqrt[0]*(random->uniform()-0.5)*g2;
|
||||
tmp[1] = g1*tmp[1]*gamma_r_inv[1] + gamma_r_invsqrt[1]*(random->uniform()-0.5)*g2;
|
||||
tmp[2] = g1*tmp[2]*gamma_r_inv[2] + gamma_r_invsqrt[2]*(random->uniform()-0.5)*g2;
|
||||
|
||||
} else if (Tp_GAUSS) {
|
||||
|
||||
tmp[0] = g1*tmp[0]*gamma_r_inv[0] + gamma_r_invsqrt[0]*random->gaussian()*g2;
|
||||
tmp[1] = g1*tmp[1]*gamma_r_inv[1] + gamma_r_invsqrt[1]*random->gaussian()*g2;
|
||||
tmp[2] = g1*tmp[2]*gamma_r_inv[2] + gamma_r_invsqrt[2]*random->gaussian()*g2;
|
||||
|
||||
} else {
|
||||
tmp[0] = g1*tmp[0]*gamma_r_inv[0];
|
||||
tmp[1] = g1*tmp[1]*gamma_r_inv[1];
|
||||
tmp[2] = g1*tmp[2]*gamma_r_inv[2];
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
// convert body frame angular velocity to quaternion
|
||||
MathExtra::quatvec(quat,tmp,qw);
|
||||
quat[0] = quat[0] + 0.5*dt*qw[0];
|
||||
quat[1] = quat[1] + 0.5*dt*qw[1];
|
||||
quat[2] = quat[2] + 0.5*dt*qw[2];
|
||||
quat[3] = quat[3] + 0.5*dt*qw[3];
|
||||
|
||||
// normalisation introduces the stochastic drift term
|
||||
// to recover the Boltzmann distribution for the case of conservative torques
|
||||
MathExtra::qnormalize(quat);
|
||||
|
||||
MathExtra::quat_to_mat( quat, Q );
|
||||
MathExtra::matvec( Q, f_act, f_rot );
|
||||
// next, update centre of mass positions and velocities
|
||||
|
||||
// tmp now holds force in body frame
|
||||
MathExtra::matvec(rotationmatrix_transpose,f[i],tmp);
|
||||
// and then converts to gamma_t^{-1} * F (velocity) in body frame
|
||||
|
||||
if (Tp_2D) {
|
||||
tmp[2] = 0.0;
|
||||
if (Tp_UNIFORM) {
|
||||
tmp[0] = g1*tmp[0]*gamma_t_inv[0]
|
||||
+ gamma_t_invsqrt[0]*(random->uniform()-0.5)*g2;
|
||||
tmp[1] = g1*tmp[1]*gamma_t_inv[1]
|
||||
+ gamma_t_invsqrt[1]*(random->uniform()-0.5)*g2;
|
||||
|
||||
} else if (Tp_GAUSS) {
|
||||
tmp[0] = g1*tmp[0]*gamma_t_inv[0]
|
||||
+ gamma_t_invsqrt[0]*random->gaussian()*g2;
|
||||
tmp[1] = g1*tmp[1]*gamma_t_inv[1]
|
||||
+ gamma_t_invsqrt[1]*random->gaussian()*g2;
|
||||
|
||||
} else {
|
||||
tmp[0] = g1*tmp[0]*gamma_t_inv[0];
|
||||
tmp[1] = g1*tmp[1]*gamma_t_inv[1];
|
||||
|
||||
}
|
||||
} else {
|
||||
if (Tp_UNIFORM) {
|
||||
|
||||
tmp[0] = g1*tmp[0]*gamma_t_inv[0]
|
||||
+ gamma_t_invsqrt[0]*(random->uniform()-0.5)*g2;
|
||||
tmp[1] = g1*tmp[1]*gamma_t_inv[1]
|
||||
+ gamma_t_invsqrt[1]*(random->uniform()-0.5)*g2;
|
||||
tmp[2] = g1*tmp[2]*gamma_t_inv[2]
|
||||
+ gamma_t_invsqrt[2]*(random->uniform()-0.5)*g2;
|
||||
|
||||
} else if (Tp_GAUSS) {
|
||||
|
||||
tmp[0] = g1*tmp[0]*gamma_t_inv[0]
|
||||
+ gamma_t_invsqrt[0]*random->gaussian()*g2;
|
||||
tmp[1] = g1*tmp[1]*gamma_t_inv[1]
|
||||
+ gamma_t_invsqrt[1]*random->gaussian()*g2;
|
||||
tmp[2] = g1*tmp[2]*gamma_t_inv[2]
|
||||
+ gamma_t_invsqrt[2]*random->gaussian()*g2;
|
||||
|
||||
|
||||
} else {
|
||||
|
||||
tmp[0] = g1*tmp[0]*gamma_t_inv[0];
|
||||
tmp[1] = g1*tmp[1]*gamma_t_inv[1];
|
||||
tmp[2] = g1*tmp[2]*gamma_t_inv[2];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
// finally, convert this back to lab-frame velocity and store in dv
|
||||
MathExtra::transpose_matvec(rotationmatrix_transpose, tmp, dv );
|
||||
|
||||
|
||||
/*
|
||||
if (Tp_UNIFORM) {
|
||||
dv[0] = dt * (g1 * f[i][0] + g2 * (random->uniform()-0.5));
|
||||
dv[1] = dt * (g1 * f[i][1] + g2 * (random->uniform()-0.5));
|
||||
dv[2] = dt * (g1 * f[i][2] + g2 * (random->uniform()-0.5));
|
||||
|
||||
} else if (Tp_GAUSS) {
|
||||
dv[0] = dt * (g1 * f[i][0] + g2 * random->gaussian());
|
||||
dv[1] = dt * (g1 * f[i][1] + g2 * random->gaussian());
|
||||
dv[2] = dt * (g1 * f[i][2] + g2 * random->gaussian());
|
||||
|
||||
} else {
|
||||
dv[0] = dt * g1 * f[i][0];
|
||||
dv[1] = dt * g1 * f[i][1];
|
||||
dv[2] = dt * g1 * f[i][2];
|
||||
|
||||
}
|
||||
*/
|
||||
v[i][0] = dv[0];
|
||||
v[i][1] = dv[1];
|
||||
v[i][2] = dv[2];
|
||||
|
||||
x[i][0] += dv[0]*dt;
|
||||
x[i][1] += dv[1]*dt;
|
||||
x[i][2] += dv[2]*dt;
|
||||
|
||||
if (Tp_DIPOLE) {
|
||||
|
||||
MathExtra::quat_to_mat_trans( quat, rotationmatrix_transpose );
|
||||
MathExtra::transpose_matvec(rotationmatrix_transpose, dipole_body, f_rot );
|
||||
|
||||
mu[i][0] = f_rot[0];
|
||||
mu[i][1] = f_rot[1];
|
||||
mu[i][2] = f_rot[2];
|
||||
|
||||
}
|
||||
}
|
||||
} else {
|
||||
|
||||
// if no dipole, just update quaternions and
|
||||
// angular velocity
|
||||
|
||||
|
||||
double wq[4];
|
||||
|
||||
double *quat;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
update_x_and_omega(x[i],v[i],omega[i],f[i],torque[i],d3rot);
|
||||
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
|
||||
MathExtra::vecquat(omega[i],quat,wq);
|
||||
|
||||
quat[0] = quat[0] + 0.5*dt*wq[0];
|
||||
quat[1] = quat[1] + 0.5*dt*wq[1];
|
||||
quat[2] = quat[2] + 0.5*dt*wq[2];
|
||||
quat[3] = quat[3] + 0.5*dt*wq[3];
|
||||
MathExtra::qnormalize(quat);
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
void FixBrownianAsphere::update_x_and_omega(double *x, double *v, double *omega,
|
||||
double *f, double *torque, int d3rot)
|
||||
{
|
||||
double dx, dy, dz;
|
||||
|
||||
dx = dt * g1 * f[0];
|
||||
x[0] += dx;
|
||||
v[0] = dx/dt;
|
||||
|
||||
dy = dt * g1 * f[1];
|
||||
x[1] += dy;
|
||||
v[1] = dy/dt;
|
||||
|
||||
dz = dt * g1 * f[2];
|
||||
x[2] += dz;
|
||||
v[2] = dz/dt;
|
||||
|
||||
omega[0] = d3rot * g3* torque[0];
|
||||
omega[1] = d3rot * g3* torque[1];
|
||||
omega[2] = g3* torque[2];
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
apply random force, stolen from MISC/fix_efield.cpp
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixBrownianAsphere::post_force(int vflag)
|
||||
{
|
||||
double **f = atom->f;
|
||||
double **x = atom->x;
|
||||
double **torque = atom->torque;
|
||||
int *mask = atom->mask;
|
||||
imageint *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
// virial setup
|
||||
|
||||
if (vflag) v_setup(vflag);
|
||||
else evflag = 0;
|
||||
|
||||
double fx,fy,fz;
|
||||
double v[6];
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
fx = g2 * (random->*rng_func)()/sqrtdt;
|
||||
fy = g2 * (random->*rng_func)()/sqrtdt;
|
||||
fz = g2 * (random->*rng_func)()/sqrtdt;
|
||||
f[i][0] += fx;
|
||||
f[i][1] += fy;
|
||||
f[i][2] += fz;
|
||||
|
||||
torque[i][0] = g4*(random->*rng_func)()/sqrtdt;
|
||||
torque[i][1] = g4*(random->*rng_func)()/sqrtdt;
|
||||
torque[i][2] = g4*(random->*rng_func)()/sqrtdt;
|
||||
|
||||
if (evflag) {
|
||||
v[0] = fx*x[i][0];
|
||||
v[1] = fy*x[i][1];
|
||||
v[2] = fz*x[i][2];
|
||||
v[3] = fx*x[i][1];
|
||||
v[4] = fx*x[i][2];
|
||||
v[5] = fy*x[i][2];
|
||||
v_tally(i, v);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void FixBrownianAsphere::reset_dt()
|
||||
{
|
||||
|
||||
dt = update->dt;
|
||||
sqrtdt = sqrt(dt);
|
||||
}
|
||||
|
||||
@ -20,43 +20,28 @@ FixStyle(brownian/asphere,FixBrownianAsphere)
|
||||
#ifndef LMP_FIX_BROWNIAN_ASPHERE_H
|
||||
#define LMP_FIX_BROWNIAN_ASPHERE_H
|
||||
|
||||
#include "fix.h"
|
||||
#include "fix_brownian_base.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixBrownianAsphere : public Fix {
|
||||
class FixBrownianAsphere : public FixBrownianBase {
|
||||
public:
|
||||
FixBrownianAsphere(class LAMMPS *, int, char **);
|
||||
virtual ~FixBrownianAsphere();
|
||||
void init();
|
||||
void initial_integrate(int);
|
||||
void setup(int);
|
||||
void post_force(int);
|
||||
int setmask();
|
||||
void reset_dt();
|
||||
void update_x_and_omega(double *, double *, double *,
|
||||
double *, double *, int );
|
||||
|
||||
void init();
|
||||
|
||||
private:
|
||||
int seed; // RNG seed
|
||||
int dipole_flag; // set if dipole is used
|
||||
double dt, sqrtdt; // time step interval and its sqrt
|
||||
|
||||
|
||||
double gamma_t,gamma_r; // translational and rotational damping params
|
||||
double diff_t,diff_r; // translational and rotational diffusion coeffs
|
||||
|
||||
double g1,g2, g3, g4; // prefactors in time stepping
|
||||
int noise_flag; // 0/1 for noise off/on
|
||||
int gaussian_noise_flag; // 0/1 for uniform/gaussian noise
|
||||
|
||||
protected:
|
||||
class AtomVecEllipsoid *avec;
|
||||
|
||||
protected:
|
||||
class RanMars *random;
|
||||
typedef double (RanMars::*rng_member)();
|
||||
rng_member rng_func; // placeholder for RNG function
|
||||
private:
|
||||
template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D >
|
||||
void initial_integrate_templated();
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
};
|
||||
|
||||
|
||||
250
src/USER-BROWNIAN/fix_brownian_base.cpp
Normal file
250
src/USER-BROWNIAN/fix_brownian_base.cpp
Normal file
@ -0,0 +1,250 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Originally modified from USER-CGDNA/fix_nve_dotc_langevin.cpp.
|
||||
|
||||
Contributing author: Sam Cameron (University of Bristol)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "fix_brownian.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include "math_extra.h"
|
||||
#include "atom.h"
|
||||
#include "force.h"
|
||||
#include "update.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "random_mars.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
|
||||
time_integrate = 1;
|
||||
|
||||
noise_flag = 1;
|
||||
gaussian_noise_flag = 0;
|
||||
gamma_t_flag = gamma_r_flag = 0;
|
||||
gamma_t_eigen_flag = gamma_r_eigen_flag = 0;
|
||||
dipole_flag = 0;
|
||||
|
||||
if (narg < 5)
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
|
||||
temp = utils::numeric(FLERR,arg[3],false,lmp);
|
||||
if (temp <= 0) error->all(FLERR,"Fix brownian temp must be > 0.");
|
||||
|
||||
seed = utils::inumeric(FLERR,arg[4],false,lmp);
|
||||
if (seed <= 0) error->all(FLERR,"Fix brownian seed must be > 0.");
|
||||
|
||||
|
||||
int iarg = 5;
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"rng") == 0) {
|
||||
if (narg == iarg + 1) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
if (strcmp(arg[iarg + 1],"uniform") == 0) {
|
||||
noise_flag = 1;
|
||||
} else if (strcmp(arg[iarg + 1],"gaussian") == 0) {
|
||||
noise_flag = 1;
|
||||
gaussian_noise_flag = 1;
|
||||
} else if (strcmp(arg[iarg + 1],"none") == 0) {
|
||||
noise_flag = 0;
|
||||
} else {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
iarg = iarg + 2;
|
||||
} else if (strcmp(arg[iarg],"dipole") == 0) {
|
||||
if (narg == iarg + 3) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
|
||||
dipole_flag = 1;
|
||||
dipole_body = new double[3];
|
||||
|
||||
dipole_body[0] = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
dipole_body[1] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
dipole_body[2] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
iarg = iarg + 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"gamma_t_eigen") == 0) {
|
||||
if (narg == iarg + 3) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
|
||||
gamma_t_eigen_flag = 1;
|
||||
gamma_t_inv = new double[3];
|
||||
gamma_t_invsqrt = new double[3];
|
||||
gamma_t_inv[0] = 1./utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
gamma_t_inv[1] = 1./utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
|
||||
if (domain->dimension == 2) {
|
||||
if (strcmp(arg[iarg+3],"inf") != 0) {
|
||||
error->all(FLERR,"Fix brownian gamma_t_eigen third value must be inf for 2D system.");
|
||||
}
|
||||
gamma_t_inv[2] = 0;
|
||||
} else {
|
||||
gamma_t_inv[2] = 1./utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
}
|
||||
|
||||
if (gamma_t_inv[0] < 0 || gamma_t_inv[1] < 0 || gamma_t_inv[2] < 0)
|
||||
error->all(FLERR,"Fix brownian gamma_t_eigen values must be > 0.");
|
||||
|
||||
gamma_t_invsqrt[0] = sqrt(gamma_t_inv[0]);
|
||||
gamma_t_invsqrt[1] = sqrt(gamma_t_inv[1]);
|
||||
gamma_t_invsqrt[2] = sqrt(gamma_t_inv[2]);
|
||||
iarg = iarg + 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"gamma_r_eigen") == 0) {
|
||||
if (narg == iarg + 3) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
|
||||
gamma_r_eigen_flag = 1;
|
||||
gamma_r_inv = new double[3];
|
||||
gamma_r_invsqrt = new double[3];
|
||||
|
||||
|
||||
if (domain->dimension == 2) {
|
||||
if (strcmp(arg[iarg+1],"inf") != 0) {
|
||||
error->all(FLERR,"Fix brownian gamma_r_eigen first value must be inf for 2D system.");
|
||||
}
|
||||
gamma_r_inv[0] = 0;
|
||||
|
||||
if (strcmp(arg[iarg+2],"inf") != 0) {
|
||||
error->all(FLERR,"Fix brownian gamma_r_eigen second value must be inf for 2D system.");
|
||||
}
|
||||
gamma_r_inv[1] = 0;
|
||||
} else {
|
||||
|
||||
gamma_r_inv[0] = 1./utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
gamma_r_inv[1] = 1./utils::numeric(FLERR,arg[iarg+2],false,lmp);
|
||||
|
||||
}
|
||||
|
||||
gamma_r_inv[2] = 1./utils::numeric(FLERR,arg[iarg+3],false,lmp);
|
||||
|
||||
if (gamma_r_inv[0] < 0 || gamma_r_inv[1] < 0 || gamma_r_inv[2] < 0)
|
||||
error->all(FLERR,"Fix brownian gamma_r_eigen values must be > 0.");
|
||||
|
||||
gamma_r_invsqrt[0] = sqrt(gamma_r_inv[0]);
|
||||
gamma_r_invsqrt[1] = sqrt(gamma_r_inv[1]);
|
||||
gamma_r_invsqrt[2] = sqrt(gamma_r_inv[2]);
|
||||
iarg = iarg + 4;
|
||||
|
||||
} else if (strcmp(arg[iarg],"gamma_t") == 0) {
|
||||
if (narg == iarg + 1) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
|
||||
gamma_t_flag = 1;
|
||||
gamma_t = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
if (gamma_t <= 0)
|
||||
error->all(FLERR,"Fix brownian gamma_t must be > 0.");
|
||||
iarg = iarg + 2;
|
||||
|
||||
} else if (strcmp(arg[iarg],"gamma_r") == 0) {
|
||||
if (narg == iarg + 1) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
|
||||
gamma_r_flag = 1;
|
||||
gamma_r = utils::numeric(FLERR,arg[iarg+1],false,lmp);
|
||||
if (gamma_r <= 0)
|
||||
error->all(FLERR,"Fix brownian gamma_r must be > 0.");
|
||||
iarg = iarg + 2;
|
||||
|
||||
} else {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
}
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
random = new RanMars(lmp,seed + comm->me);
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBrownianBase::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBrownianBase::~FixBrownianBase()
|
||||
{
|
||||
|
||||
if (gamma_t_eigen_flag) {
|
||||
delete [] gamma_t_inv;
|
||||
delete [] gamma_t_invsqrt;
|
||||
}
|
||||
if (gamma_r_eigen_flag) {
|
||||
delete [] gamma_r_inv;
|
||||
delete [] gamma_r_invsqrt;
|
||||
}
|
||||
|
||||
if (dipole_flag) {
|
||||
delete [] dipole_body;
|
||||
}
|
||||
|
||||
|
||||
delete random;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBrownianBase::init()
|
||||
{
|
||||
dt = update->dt;
|
||||
sqrtdt = sqrt(dt);
|
||||
|
||||
g1 = force->ftm2v;
|
||||
if (noise_flag == 0) {
|
||||
g2 = 0;
|
||||
} else if (gaussian_noise_flag == 1) {
|
||||
g2 = sqrt(2 * force->boltz*temp/dt/force->mvv2e);
|
||||
} else {
|
||||
g2 = sqrt( 24 * force->boltz*temp/dt/force->mvv2e);
|
||||
}
|
||||
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
void FixBrownianBase::reset_dt()
|
||||
{
|
||||
double sqrtdt_old = sqrtdt;
|
||||
dt = update->dt;
|
||||
sqrtdt = sqrt(dt);
|
||||
g2 *= sqrtdt_old/sqrtdt;
|
||||
|
||||
}
|
||||
100
src/USER-BROWNIAN/fix_brownian_base.h
Normal file
100
src/USER-BROWNIAN/fix_brownian_base.h
Normal file
@ -0,0 +1,100 @@
|
||||
/* -*- 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.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
|
||||
#ifndef LMP_FIX_BROWNIAN_BASE_H
|
||||
#define LMP_FIX_BROWNIAN_BASE_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
|
||||
class FixBrownianBase : public Fix {
|
||||
public:
|
||||
FixBrownianBase(class LAMMPS *, int, char **);
|
||||
virtual ~FixBrownianBase();
|
||||
void init();
|
||||
int setmask();
|
||||
void reset_dt();
|
||||
|
||||
protected:
|
||||
int seed; // RNG seed
|
||||
double dt, sqrtdt; // time step interval and its sqrt
|
||||
|
||||
int gamma_t_flag; // 0/1 if isotropic translational damping is unset/set
|
||||
int gamma_r_flag; // 0/1 if isotropic rotational damping is unset/set
|
||||
double gamma_t,gamma_r; // translational and rotational (isotropic) damping params
|
||||
|
||||
|
||||
int gamma_t_eigen_flag; // 0/1 if anisotropic translational damping is unset/set
|
||||
int gamma_r_eigen_flag; // 0/1 if anisotropic rotational damping is unset/set
|
||||
|
||||
double *gamma_t_inv; // anisotropic damping parameter eigenvalues
|
||||
double *gamma_r_inv;
|
||||
double *gamma_t_invsqrt;
|
||||
double *gamma_r_invsqrt;
|
||||
|
||||
|
||||
int dipole_flag; // set if dipole is used for asphere
|
||||
double *dipole_body; // direction dipole is slaved to in body frame
|
||||
|
||||
int noise_flag; // 0/1 for noise off/on
|
||||
int gaussian_noise_flag; // 0/1 for uniform/gaussian noise
|
||||
|
||||
double temp; // temperature
|
||||
|
||||
|
||||
|
||||
double g1,g2; // prefactors in time stepping
|
||||
|
||||
|
||||
class RanMars *random;
|
||||
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal fix brownian command.
|
||||
|
||||
Wrong number/type of input arguments.
|
||||
|
||||
E: Fix brownian gamma_t_eigen values must be > 0.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix brownian gamma_r_eigen values must be > 0.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix brownian seed must be > 0.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix brownian temp must be > 0.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix brownian gamma_t must be > 0.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix brownian gamma_r must be > 0.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
*/
|
||||
@ -35,103 +35,34 @@
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
#define SMALL 1e-14
|
||||
|
||||
enum{NONE,DIPOLE};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBrownianSphere::FixBrownianSphere(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
FixBrownianBase(lmp, narg, arg)
|
||||
{
|
||||
virial_flag = 1;
|
||||
|
||||
time_integrate = 1;
|
||||
|
||||
extra = NONE;
|
||||
|
||||
if (narg > 11 || narg < 8 )
|
||||
error->all(FLERR,"Illegal fix brownian/sphere command.");
|
||||
|
||||
if (!atom->sphere_flag)
|
||||
error->all(FLERR,"Fix brownian/sphere requires atom style sphere");
|
||||
|
||||
gamma_t = utils::numeric(FLERR,arg[3],false,lmp);
|
||||
if (gamma_t <= 0.0)
|
||||
error->all(FLERR,"Fix brownian/sphere translational viscous drag "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
gamma_r = utils::numeric(FLERR,arg[4],false,lmp);
|
||||
if (gamma_t <= 0.0)
|
||||
error->all(FLERR,"Fix brownian/sphere rotational viscous drag "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
|
||||
diff_t = utils::numeric(FLERR,arg[5],false,lmp);
|
||||
if (diff_t <= 0.0)
|
||||
error->all(FLERR,"Fix brownian/sphere translational diffusion "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
diff_r = utils::numeric(FLERR,arg[6],false,lmp);
|
||||
if (diff_r <= 0.0)
|
||||
error->all(FLERR,"Fix brownian/sphere rotational diffusion "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
seed = utils::inumeric(FLERR,arg[7],false,lmp);
|
||||
if (seed <= 0) error->all(FLERR,"Fix brownian/sphere seed must be > 0.");
|
||||
|
||||
noise_flag = 1;
|
||||
gaussian_noise_flag = 0;
|
||||
|
||||
int iarg = 8;
|
||||
|
||||
while (iarg < narg) {
|
||||
if (strcmp(arg[iarg],"rng") == 0) {
|
||||
if (narg == iarg + 1) {
|
||||
error->all(FLERR,"Illegal fix/brownian/sphere command.");
|
||||
}
|
||||
if (strcmp(arg[iarg + 1],"uniform") == 0) {
|
||||
noise_flag = 1;
|
||||
} else if (strcmp(arg[iarg + 1],"gaussian") == 0) {
|
||||
noise_flag = 1;
|
||||
gaussian_noise_flag = 1;
|
||||
} else if (strcmp(arg[iarg + 1],"none") == 0) {
|
||||
noise_flag = 0;
|
||||
} else {
|
||||
error->all(FLERR,"Illegal fix/brownian/sphere command.");
|
||||
}
|
||||
iarg = iarg + 2;
|
||||
} else if (strcmp(arg[iarg],"dipole") == 0) {
|
||||
extra = DIPOLE;
|
||||
iarg = iarg + 1;
|
||||
} else {
|
||||
error->all(FLERR,"Illegal fix/brownian/sphere command.");
|
||||
}
|
||||
if (gamma_t_eigen_flag || gamma_r_eigen_flag) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
|
||||
if (extra == DIPOLE && !atom->mu_flag)
|
||||
error->all(FLERR,"Fix brownian/sphere update dipole requires atom attribute mu");
|
||||
if (!gamma_t_flag || !gamma_r_flag) {
|
||||
error->all(FLERR,"Illegal fix brownian command.");
|
||||
}
|
||||
|
||||
|
||||
if (!atom->mu_flag)
|
||||
error->all(FLERR,"Fix brownian/sphere requires atom attribute mu");
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
random = new RanMars(lmp,seed + comm->me);
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBrownianSphere::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= POST_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBrownianSphere::~FixBrownianSphere()
|
||||
{
|
||||
delete random;
|
||||
}
|
||||
|
||||
|
||||
@ -141,181 +72,161 @@ FixBrownianSphere::~FixBrownianSphere()
|
||||
void FixBrownianSphere::init()
|
||||
{
|
||||
|
||||
g1 = force->ftm2v/gamma_t;
|
||||
g3 = force->ftm2v/gamma_r;
|
||||
if (noise_flag == 0) {
|
||||
g2 = 0;
|
||||
g4 = 0;
|
||||
rng_func = &RanMars::zero_rng;
|
||||
} else if (gaussian_noise_flag == 1) {
|
||||
g2 = gamma_t*sqrt(2 * diff_t)/force->ftm2v;
|
||||
g4 = gamma_r*sqrt(2 * diff_r)/force->ftm2v;
|
||||
rng_func = &RanMars::gaussian;
|
||||
} else {
|
||||
g2 = gamma_t*sqrt( 24 * diff_t)/force->ftm2v;
|
||||
g4 = gamma_r*sqrt( 24 * diff_r )/force->ftm2v;
|
||||
rng_func = &RanMars::uniform_middle;
|
||||
}
|
||||
FixBrownianBase::init();
|
||||
|
||||
g3 = g1/gamma_r;
|
||||
g4 = g2/sqrt(gamma_r);
|
||||
|
||||
g1 /= gamma_t;
|
||||
g2 /= sqrt(gamma_t);
|
||||
|
||||
dt = update->dt;
|
||||
sqrtdt = sqrt(dt);
|
||||
}
|
||||
|
||||
void FixBrownianSphere::setup(int vflag)
|
||||
{
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBrownianSphere::initial_integrate(int /* vflag */)
|
||||
void FixBrownianSphere::initial_integrate(int /*vflag */)
|
||||
{
|
||||
double **x = atom->x;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
double **omega = atom->omega;
|
||||
double **torque = atom->torque;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
double dx,dy,dz;
|
||||
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
if (domain->dimension == 2) {
|
||||
|
||||
int d3rot; // whether to compute angular momentum in xy plane
|
||||
|
||||
if (domain->dimension==2) {
|
||||
d3rot = 0;
|
||||
if (!noise_flag) {
|
||||
initial_integrate_templated<0,0,1>();
|
||||
} else if (gaussian_noise_flag) {
|
||||
initial_integrate_templated<0,1,1>();
|
||||
} else {
|
||||
d3rot = 1;
|
||||
initial_integrate_templated<1,0,1>();
|
||||
}
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
dx = dt * g1 * f[i][0];
|
||||
x[i][0] += dx;
|
||||
v[i][0] = dx/dt;
|
||||
|
||||
dy = dt * g1 * f[i][1];
|
||||
x[i][1] += dy;
|
||||
v[i][1] = dy/dt;
|
||||
|
||||
dz = dt * g1 * f[i][2];
|
||||
x[i][2] += dz;
|
||||
v[i][2] = dz/dt;
|
||||
|
||||
omega[i][0] = d3rot * g3* torque[i][0];
|
||||
omega[i][1] = d3rot * g3* torque[i][1];
|
||||
omega[i][2] = g3* torque[i][2];
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
if (extra == DIPOLE) {
|
||||
|
||||
double **mu = atom->mu;
|
||||
double dtheta;
|
||||
double mux,muy,muz,mu_tmp,wx,wy,wz;
|
||||
double prefac_1, prefac_2;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
dtheta = sqrt((omega[i][0]*dt)*(omega[i][0]*dt)
|
||||
+(omega[i][1]*dt)*(omega[i][1]*dt)
|
||||
+(omega[i][2]*dt)*(omega[i][2]*dt));
|
||||
|
||||
|
||||
if (fabs(dtheta) < SMALL) {
|
||||
prefac_1 = dt;
|
||||
prefac_2 = 0.5*dt*dt;
|
||||
} else {
|
||||
prefac_1 = dt*sin(dtheta)/dtheta;
|
||||
prefac_2 = dt*dt*(1-cos(dtheta))/(dtheta*dtheta);
|
||||
}
|
||||
|
||||
mux = mu[i][0];
|
||||
muy = mu[i][1];
|
||||
muz = mu[i][2];
|
||||
|
||||
wx = omega[i][0];
|
||||
wy = omega[i][1];
|
||||
wz = omega[i][2];
|
||||
|
||||
mu[i][0] = (mux + prefac_1 * ( -wz*muy + wy*muz )
|
||||
+ prefac_2 * ( -1*( wz*wz + wy*wy ) * mux
|
||||
+ ( wz*muz + wy*muy ) * wx));
|
||||
|
||||
mu[i][1] = (muy + prefac_1 * ( wz*mux - wx*muz )
|
||||
+ prefac_2 * ( -1*(wz*wz + wx*wx) * muy
|
||||
+ ( wz*muz + wx*mux ) * wy));
|
||||
|
||||
mu[i][2] = (muz + prefac_1 * ( -wy*mux + wx*muy )
|
||||
+ prefac_2 * ( -1*( wx*wx + wy*wy ) * muz
|
||||
+ ( wy*muy + wx*mux ) * wz));
|
||||
|
||||
mu_tmp = sqrt(mu[i][0]*mu[i][0]+mu[i][1]*mu[i][1]+mu[i][2]*mu[i][2]);
|
||||
|
||||
mu[i][0] = mu[i][0]/mu_tmp;
|
||||
mu[i][1] = mu[i][1]/mu_tmp;
|
||||
mu[i][2] = mu[i][2]/mu_tmp;
|
||||
}
|
||||
if (!noise_flag) {
|
||||
initial_integrate_templated<0,0,0>();
|
||||
} else if (gaussian_noise_flag) {
|
||||
initial_integrate_templated<0,1,0>();
|
||||
} else {
|
||||
initial_integrate_templated<1,0,0>();
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
apply random force, stolen from MISC/fix_efield.cpp
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void FixBrownianSphere::post_force(int vflag)
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D >
|
||||
void FixBrownianSphere::initial_integrate_templated()
|
||||
{
|
||||
double **f = atom->f;
|
||||
double **x = atom->x;
|
||||
double **torque = atom->torque;
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
int *mask = atom->mask;
|
||||
imageint *image = atom->image;
|
||||
int nlocal = atom->nlocal;
|
||||
double wx,wy,wz;
|
||||
double **torque = atom->torque;
|
||||
double **mu = atom->mu;
|
||||
double mux,muy,muz,mulen;
|
||||
|
||||
// virial setup
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
if (vflag) v_setup(vflag);
|
||||
else evflag = 0;
|
||||
|
||||
double fx,fy,fz;
|
||||
double v[6];
|
||||
double dx,dy,dz;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
|
||||
fx = g2 * (random->*rng_func)()/sqrtdt;
|
||||
fy = g2 * (random->*rng_func)()/sqrtdt;
|
||||
fz = g2 * (random->*rng_func)()/sqrtdt;
|
||||
f[i][0] += fx;
|
||||
f[i][1] += fy;
|
||||
f[i][2] += fz;
|
||||
if (Tp_2D) {
|
||||
dz = 0;
|
||||
wx = wy = 0;
|
||||
if (Tp_UNIFORM) {
|
||||
dx = dt * (g1 * f[i][0] + g2 * (random->uniform()-0.5));
|
||||
dy = dt * (g1 * f[i][1] + g2 * (random->uniform()-0.5));
|
||||
wz = (random->uniform()-0.5)*g4;
|
||||
|
||||
torque[i][0] = g4*(random->*rng_func)()/sqrtdt;
|
||||
torque[i][1] = g4*(random->*rng_func)()/sqrtdt;
|
||||
torque[i][2] = g4*(random->*rng_func)()/sqrtdt;
|
||||
} else if (Tp_GAUSS) {
|
||||
dx = dt * (g1 * f[i][0] + g2 * random->gaussian());
|
||||
dy = dt * (g1 * f[i][1] + g2 * random->gaussian());
|
||||
wz = random->gaussian()*g4;
|
||||
|
||||
|
||||
} else {
|
||||
dx = dt * g1 * f[i][0];
|
||||
dy = dt * g1 * f[i][1];
|
||||
wz = 0;
|
||||
|
||||
}
|
||||
} else {
|
||||
if (Tp_UNIFORM) {
|
||||
dx = dt * (g1 * f[i][0] + g2 * (random->uniform()-0.5));
|
||||
dy = dt * (g1 * f[i][1] + g2 * (random->uniform()-0.5));
|
||||
dz = dt * (g1 * f[i][2] + g2 * (random->uniform()-0.5));
|
||||
wx = (random->uniform()-0.5)*g4;
|
||||
wy = (random->uniform()-0.5)*g4;
|
||||
wz = (random->uniform()-0.5)*g4;
|
||||
|
||||
} else if (Tp_GAUSS) {
|
||||
dx = dt * (g1 * f[i][0] + g2 * random->gaussian());
|
||||
dy = dt * (g1 * f[i][1] + g2 * random->gaussian());
|
||||
dz = dt * (g1 * f[i][2] + g2 * random->gaussian());
|
||||
wx = random->gaussian()*g4;
|
||||
wy = random->gaussian()*g4;
|
||||
wz = random->gaussian()*g4;
|
||||
|
||||
|
||||
} else {
|
||||
dx = dt * g1 * f[i][0];
|
||||
dy = dt * g1 * f[i][1];
|
||||
dz = dt * g1 * f[i][2];
|
||||
wx = wy = wz = 0;
|
||||
|
||||
if (evflag) {
|
||||
v[0] = fx*x[i][0];
|
||||
v[1] = fy*x[i][1];
|
||||
v[2] = fz*x[i][2];
|
||||
v[3] = fx*x[i][1];
|
||||
v[4] = fx*x[i][2];
|
||||
v[5] = fy*x[i][2];
|
||||
v_tally(i, v);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void FixBrownianSphere::reset_dt()
|
||||
{
|
||||
x[i][0] += dx;
|
||||
v[i][0] = dx/dt;
|
||||
|
||||
dt = update->dt;
|
||||
sqrtdt = sqrt(dt);
|
||||
|
||||
x[i][1] += dy;
|
||||
v[i][1] = dy/dt;
|
||||
|
||||
|
||||
x[i][2] += dz;
|
||||
v[i][2] = dz/dt;
|
||||
|
||||
|
||||
|
||||
wx += g3*torque[i][0];
|
||||
wy += g3*torque[i][1];
|
||||
wz += g3*torque[i][2];
|
||||
|
||||
|
||||
// store length of dipole as we need to convert it to a unit vector and
|
||||
// then back again
|
||||
|
||||
mulen = sqrt(mu[i][0]*mu[i][0] + mu[i][1]*mu[i][1] + mu[i][2]*mu[i][2]);
|
||||
|
||||
// unit vector at time t
|
||||
mux = mu[i][0]/mulen;
|
||||
muy = mu[i][1]/mulen;
|
||||
muz = mu[i][2]/mulen;
|
||||
|
||||
|
||||
|
||||
// un-normalised unit vector at time t + dt
|
||||
mu[i][0] = mux + (wy*muz - wz*muy)*dt;
|
||||
mu[i][1] = muy + (wz*mux - wx*muz)*dt;
|
||||
mu[i][2] = muz + (wx*muy - wy*mux)*dt;
|
||||
|
||||
// normalisation introduces the stochastic drift term due to changing from
|
||||
// Stratonovich to Ito interpretation
|
||||
MathExtra::norm3(mu[i]);
|
||||
|
||||
// multiply by original magnitude to obtain dipole of same length
|
||||
mu[i][0] = mu[i][0]*mulen;
|
||||
mu[i][1] = mu[i][1]*mulen;
|
||||
mu[i][2] = mu[i][2]*mulen;
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
@ -20,39 +20,22 @@ FixStyle(brownian/sphere,FixBrownianSphere)
|
||||
#ifndef LMP_FIX_BROWNIAN_SPHERE_H
|
||||
#define LMP_FIX_BROWNIAN_SPHERE_H
|
||||
|
||||
#include "fix.h"
|
||||
#include "fix_brownian_base.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixBrownianSphere : public Fix {
|
||||
class FixBrownianSphere : public FixBrownianBase {
|
||||
public:
|
||||
FixBrownianSphere(class LAMMPS *, int, char **);
|
||||
virtual ~FixBrownianSphere();
|
||||
void init();
|
||||
void initial_integrate(int);
|
||||
void setup(int);
|
||||
void post_force(int);
|
||||
int setmask();
|
||||
void reset_dt();
|
||||
void initial_integrate(int /*vflag */);
|
||||
|
||||
|
||||
private:
|
||||
int seed; // RNG seed
|
||||
int extra; // set if dipole is used
|
||||
double dt, sqrtdt; // time step interval and its sqrt
|
||||
|
||||
|
||||
double gamma_t,gamma_r; // translational and rotational damping params
|
||||
double diff_t,diff_r; // translational and rotational diffusion coeffs
|
||||
|
||||
double g1,g2, g3, g4; // prefactors in time stepping
|
||||
int noise_flag; // 0/1 for noise off/on
|
||||
int gaussian_noise_flag; // 0/1 for uniform/gaussian noise
|
||||
|
||||
protected:
|
||||
class RanMars *random;
|
||||
typedef double (RanMars::*rng_member)();
|
||||
rng_member rng_func; // placeholder for RNG function
|
||||
|
||||
template < int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D >
|
||||
void initial_integrate_templated();
|
||||
double g3,g4;
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
@ -94,21 +94,6 @@ double RanMars::uniform()
|
||||
return uni;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
uniform RN shifted to be symmetric about zero (for fix bd/sphere).
|
||||
------------------------------------------------------------------------- */
|
||||
double RanMars::uniform_middle()
|
||||
{
|
||||
return uniform()-0.5;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Return 0 (for fix/bd/sphere).
|
||||
------------------------------------------------------------------------- */
|
||||
double RanMars::zero_rng()
|
||||
{
|
||||
return 0.0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
gaussian RN
|
||||
|
||||
@ -23,8 +23,6 @@ class RanMars : protected Pointers {
|
||||
RanMars(class LAMMPS *, int);
|
||||
~RanMars();
|
||||
double uniform();
|
||||
double uniform_middle();
|
||||
double zero_rng();
|
||||
double gaussian();
|
||||
double gaussian(double mu, double sigma);
|
||||
double rayleigh(double sigma);
|
||||
|
||||
Reference in New Issue
Block a user