diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 9ee699faea..4a0b44fd90 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -40,8 +40,8 @@ OPT. * :doc:`aveforce ` * :doc:`balance ` * :doc:`brownian ` - * :doc:`brownian/asphere ` - * :doc:`brownian/sphere ` + * :doc:`brownian/asphere ` + * :doc:`brownian/sphere ` * :doc:`bocs ` * :doc:`bond/break ` * :doc:`bond/create ` diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 373e5ac29f..69c381bfa2 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -183,8 +183,8 @@ accelerated styles exist. * :doc:`aveforce ` - add an averaged force to each atom * :doc:`balance ` - perform dynamic load-balancing * :doc:`brownian ` - overdamped translational brownian motion -* :doc:`brownian/asphere ` - overdamped translational and rotational brownian motion for ellipsoids -* :doc:`brownian/sphere ` - overdamped translational and rotational brownian motion for spheres +* :doc:`brownian/asphere ` - overdamped translational and rotational brownian motion for ellipsoids +* :doc:`brownian/sphere ` - overdamped translational and rotational brownian motion for spheres * :doc:`bocs ` - NPT style time integration with pressure correction * :doc:`bond/break ` - break bonds on the fly * :doc:`bond/create ` - create bonds on the fly diff --git a/doc/src/fix_brownian.rst b/doc/src/fix_brownian.rst index e3f3a2679f..7a82c2ffea 100644 --- a/doc/src/fix_brownian.rst +++ b/doc/src/fix_brownian.rst @@ -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 ` 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) `). +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) ` or :ref:`(Delong) `). + +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 ` and +:doc:`atom_style sphere `. + +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) Temperature computation using the :doc:`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) ` 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 ` 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) ` for why this works). This is the same method +:ref:`(Dunweg) ` for why this works). This is the same method of noise generation as used in :doc:`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) `. +value for reasons argued in :ref:`(Dunweg) `. 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 No global or per-atom quantities are stored by this fix for access by various :doc:`output commands `. -The :doc:`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 `. -The default is *virial no*. No parameter of this fix can be used with the *start/stop* keywords of the :doc:`run ` command. This fix is not invoked during @@ -117,6 +176,13 @@ the :doc:`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 ` command. +The style *brownian/asphere* fix requires that atoms store torque and quaternions +as defined by the :doc:`atom_style ellipsoid ` command. +If the *dipole* keyword is used, they must also store a dipole moment +as defined by the :doc:`atom_style dipole ` 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 ` doc page for more info. @@ -125,26 +191,28 @@ doc page for more info. Related commands """""""""""""""" +:doc:`fix propel/self `, :doc:`fix langevin `, :doc:`fix nve/sphere `, -:doc:`fix brownian/sphere `, -:doc:`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). diff --git a/doc/src/fix_brownian_asphere.rst b/doc/src/fix_brownian_asphere.rst deleted file mode 100644 index c59774ca65..0000000000 --- a/doc/src/fix_brownian_asphere.rst +++ /dev/null @@ -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 ` 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) `), :math:`dW_t` and -:math:`dW_r` are Wiener processes (see e.g. :ref:`(Gardiner) `). -The quaternions :math:`q` of the ellipsoid are updated each timestep from -the angular velocity vector. - -See :doc:`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) ` for why this works). This is the same method -of noise generation as used in :doc:`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) `. - -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 `. -No global or per-atom quantities are stored -by this fix for access by various :doc:`output commands `. - -The :doc:`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 `. -The default is *virial no*. - -No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is not invoked during -:doc:`energy minimization `. - - -Restrictions -"""""""""""" - -This fix requires that atoms store torque and angular velocity (omega) -as defined by the :doc:`atom_style sphere ` command, as well -as atoms which have a definite orientation as defined by the -:doc:`atom_style ellipsoid ` command. -Optionally, they can also store a dipole moment as defined by the -:doc:`atom_style dipole ` 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 ` -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 `, :doc:`fix brownian/sphere `, -:doc:`fix propel/self `, :doc:`fix langevin `, -:doc:`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). - - diff --git a/doc/src/fix_brownian_sphere.rst b/doc/src/fix_brownian_sphere.rst deleted file mode 100644 index 7e73226e08..0000000000 --- a/doc/src/fix_brownian_sphere.rst +++ /dev/null @@ -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 ` 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) `), :math:`dW_t` and -:math:`dW_r` are Wiener processes (see e.g. :ref:`(Gardiner) `). -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) `, -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 ` -with this fix in :doc:`fix brownian `. - ---------- - -.. note:: - The diffusion coefficient :math:`D_t` and the translational - drag coefficient :math:`\gamma_t` are discussed in - :doc:`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 ` 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) ` for why this works). This is the same method -of noise generation as used in :doc:`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) `. - -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 `. -No global or per-atom quantities are stored -by this fix for access by various :doc:`output commands `. - -The :doc:`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 `. -The default is *virial no*. - -No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. This fix is not invoked during -:doc:`energy minimization `. - -Restrictions -"""""""""""" - -This fix requires that atoms store torque and angular velocity (omega) -as defined by the :doc:`atom_style sphere ` command. -If the *dipole* keyword is used, they must also store a dipole moment -as defined by the :doc:`atom_style dipole ` 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 ` -doc page for more info. - - -Related commands -"""""""""""""""" - -:doc:`fix brownian `, :doc:`fix brownian/asphere `, -:doc:`fix propel/self `, -:doc:`fix langevin `, :doc:`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). - - diff --git a/examples/USER/brownian/2d_velocity/in.2d_velocity b/examples/USER/brownian/2d_velocity/in2d.velocity similarity index 80% rename from examples/USER/brownian/2d_velocity/in.2d_velocity rename to examples/USER/brownian/2d_velocity/in2d.velocity index cbd7d543cb..454396b959 100644 --- a/examples/USER/brownian/2d_velocity/in.2d_velocity +++ b/examples/USER/brownian/2d_velocity/in2d.velocity @@ -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 diff --git a/examples/USER/brownian/2d_velocity/log_1_1_4_2d.lammps.log b/examples/USER/brownian/2d_velocity/log_1_1_4_2d.lammps.log index f6b234e9ee..655b04f665 100644 --- a/examples/USER/brownian/2d_velocity/log_1_1_4_2d.lammps.log +++ b/examples/USER/brownian/2d_velocity/log_1_1_4_2d.lammps.log @@ -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 diff --git a/examples/USER/brownian/asphere/README.txt b/examples/USER/brownian/asphere/README.txt deleted file mode 100644 index 47e79353c6..0000000000 --- a/examples/USER/brownian/asphere/README.txt +++ /dev/null @@ -1,2 +0,0 @@ -The input file in3d.brownian demonstrates how to run a 3d simulation -of ellipsoidal particles undergoing overdamped brownian motion. diff --git a/examples/USER/brownian/asphere/in2d.ellipsoid b/examples/USER/brownian/asphere/in2d.ellipsoid new file mode 100644 index 0000000000..f36ab3b19b --- /dev/null +++ b/examples/USER/brownian/asphere/in2d.ellipsoid @@ -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 diff --git a/examples/USER/brownian/asphere/in3d.brownian b/examples/USER/brownian/asphere/in3d.brownian deleted file mode 100644 index 19487978ce..0000000000 --- a/examples/USER/brownian/asphere/in3d.brownian +++ /dev/null @@ -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 diff --git a/examples/USER/brownian/asphere/in3d.ellipsoid b/examples/USER/brownian/asphere/in3d.ellipsoid new file mode 100644 index 0000000000..b8aa230a61 --- /dev/null +++ b/examples/USER/brownian/asphere/in3d.ellipsoid @@ -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 diff --git a/examples/USER/brownian/asphere/log_gaussian_1_0.33_1_3_3d.lammps.log b/examples/USER/brownian/asphere/log_gaussian_1_0.33_1_3_3d.lammps.log deleted file mode 100644 index 2ed81ec8ba..0000000000 --- a/examples/USER/brownian/asphere/log_gaussian_1_0.33_1_3_3d.lammps.log +++ /dev/null @@ -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 diff --git a/examples/USER/brownian/asphere/log_gaussian_4_1_7_13_3d.lammps.log b/examples/USER/brownian/asphere/log_gaussian_4_1_7_13_3d.lammps.log deleted file mode 100644 index 21bd04f7da..0000000000 --- a/examples/USER/brownian/asphere/log_gaussian_4_1_7_13_3d.lammps.log +++ /dev/null @@ -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 diff --git a/examples/USER/brownian/point/in2d.point b/examples/USER/brownian/point/in2d.point new file mode 100644 index 0000000000..e502bc0b98 --- /dev/null +++ b/examples/USER/brownian/point/in2d.point @@ -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 + diff --git a/examples/USER/brownian/point/in3d.point b/examples/USER/brownian/point/in3d.point new file mode 100644 index 0000000000..b30d01e4d1 --- /dev/null +++ b/examples/USER/brownian/point/in3d.point @@ -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 + diff --git a/examples/USER/brownian/sphere/README.txt b/examples/USER/brownian/sphere/README.txt deleted file mode 100644 index 8987a23601..0000000000 --- a/examples/USER/brownian/sphere/README.txt +++ /dev/null @@ -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 /(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 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. diff --git a/examples/USER/brownian/sphere/in2d.sphere b/examples/USER/brownian/sphere/in2d.sphere new file mode 100644 index 0000000000..f0d5c26262 --- /dev/null +++ b/examples/USER/brownian/sphere/in2d.sphere @@ -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 diff --git a/examples/USER/brownian/sphere/in2ddipole.brownian b/examples/USER/brownian/sphere/in2ddipole.brownian deleted file mode 100644 index e7a8374984..0000000000 --- a/examples/USER/brownian/sphere/in2ddipole.brownian +++ /dev/null @@ -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 - - diff --git a/examples/USER/brownian/sphere/in3d.sphere b/examples/USER/brownian/sphere/in3d.sphere new file mode 100644 index 0000000000..06ca649c4c --- /dev/null +++ b/examples/USER/brownian/sphere/in3d.sphere @@ -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 diff --git a/examples/USER/brownian/sphere/in3d_virial_on.brownian b/examples/USER/brownian/sphere/in3d_virial_on.brownian deleted file mode 100644 index 2fb08cbb17..0000000000 --- a/examples/USER/brownian/sphere/in3d_virial_on.brownian +++ /dev/null @@ -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 - - diff --git a/examples/USER/brownian/sphere/log_gaussian_4_1_7_13_2d.lammps.log b/examples/USER/brownian/sphere/log_gaussian_4_1_7_13_2d.lammps.log deleted file mode 100644 index b10f77f0cd..0000000000 --- a/examples/USER/brownian/sphere/log_gaussian_4_1_7_13_2d.lammps.log +++ /dev/null @@ -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 diff --git a/examples/USER/brownian/sphere/log_uniform_10_1_5_15_3d.lammps.log b/examples/USER/brownian/sphere/log_uniform_10_1_5_15_3d.lammps.log deleted file mode 100644 index 57f90193ac..0000000000 --- a/examples/USER/brownian/sphere/log_uniform_10_1_5_15_3d.lammps.log +++ /dev/null @@ -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 diff --git a/examples/USER/brownian/spherical_ABP/in.2d_abp b/examples/USER/brownian/spherical_ABP/in2d.abp similarity index 81% rename from examples/USER/brownian/spherical_ABP/in.2d_abp rename to examples/USER/brownian/spherical_ABP/in2d.abp index 9db5a46cf6..7c53717038 100644 --- a/examples/USER/brownian/spherical_ABP/in.2d_abp +++ b/examples/USER/brownian/spherical_ABP/in2d.abp @@ -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 diff --git a/examples/USER/brownian/spherical_ABP/in.3d_ideal_abp b/examples/USER/brownian/spherical_ABP/in3d.ideal_abp similarity index 78% rename from examples/USER/brownian/spherical_ABP/in.3d_ideal_abp rename to examples/USER/brownian/spherical_ABP/in3d.ideal_abp index 4c3d403f20..dfe782525c 100644 --- a/examples/USER/brownian/spherical_ABP/in.3d_ideal_abp +++ b/examples/USER/brownian/spherical_ABP/in3d.ideal_abp @@ -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 diff --git a/examples/USER/brownian/translational/in.brownian b/examples/USER/brownian/translational/in.brownian deleted file mode 100644 index fbf43e68c6..0000000000 --- a/examples/USER/brownian/translational/in.brownian +++ /dev/null @@ -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 - - diff --git a/examples/USER/brownian/translational/log_1_1.lammps.log b/examples/USER/brownian/translational/log_1_1.lammps.log deleted file mode 100644 index 0fdba989bd..0000000000 --- a/examples/USER/brownian/translational/log_1_1.lammps.log +++ /dev/null @@ -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 diff --git a/src/USER-BROWNIAN/fix_brownian.cpp b/src/USER-BROWNIAN/fix_brownian.cpp index f25f6d6b32..128bd92042 100644 --- a/src/USER-BROWNIAN/fix_brownian.cpp +++ b/src/USER-BROWNIAN/fix_brownian.cpp @@ -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) + if (dipole_flag || gamma_t_eigen_flag || gamma_r_eigen_flag || gamma_r_flag) { 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 { - error->all(FLERR,"Illegal fix brownian command."); - } - } else { - error->all(FLERR,"Illegal fix brownian command."); - } + 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; } @@ -109,109 +64,99 @@ FixBrownian::~FixBrownian() void FixBrownian::init() { + + FixBrownianBase::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; - } + g1 /= gamma_t; + g2 *= sqrt(gamma_t); - dt = update->dt; - sqrtdt = sqrt(dt); -} - -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]; + } + } - 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; - + + } } - - 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); -} diff --git a/src/USER-BROWNIAN/fix_brownian.h b/src/USER-BROWNIAN/fix_brownian.h index fb56edb39e..7260eb15e3 100644 --- a/src/USER-BROWNIAN/fix_brownian.h +++ b/src/USER-BROWNIAN/fix_brownian.h @@ -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 - }; } diff --git a/src/USER-BROWNIAN/fix_brownian_asphere.cpp b/src/USER-BROWNIAN/fix_brownian_asphere.cpp index 123a7a5f5f..98c201230a 100644 --- a/src/USER-BROWNIAN/fix_brownian_asphere.cpp +++ b/src/USER-BROWNIAN/fix_brownian_asphere.cpp @@ -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 (gamma_t_flag || gamma_r_flag) { + error->all(FLERR,"Illegal fix brownian command."); + } + - if (dipole_flag == DIPOLE && !atom->mu_flag) + + 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; + + AtomVecEllipsoid::Bonus *bonus = avec->bonus; + + + double **mu = atom->mu; + double **torque = atom->torque; + + double qw[4]; + + double *quat; + int *ellipsoid = atom->ellipsoid; if (igroup == atom->firstgroup) nlocal = atom->nfirst; - int d3rot; // whether to compute angular momentum in xy plane + + + // project dipole along x axis of quat + double f_rot[3]; - if (domain->dimension==2) { - d3rot = 0; - } else { - d3rot = 1; - } - - if (dipole_flag == DIPOLE) { - - // if dipole is being tracked, then update it along with - // quaternions accordingly along with angular velocity - - double wq[4]; + double rotationmatrix_transpose[3][3]; - double *quat; - int *ellipsoid = atom->ellipsoid; - AtomVecEllipsoid::Bonus *bonus = avec->bonus; - - // project dipole along x axis of quat - double f_act[3] = { 1.0, 0.0, 0.0 }; - double f_rot[3]; + double tmp[3]; + double dv[3]; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + + // update orientation first + + quat = bonus[ellipsoid[i]].quat; + + MathExtra::quat_to_mat_trans( quat, rotationmatrix_transpose ); + + // tmp holds angular velocity in body frame + MathExtra::matvec(rotationmatrix_transpose,torque[i],tmp); - double Q[3][3]; - - double **mu = atom->mu; - - - for (int i = 0; i < nlocal; i++) { - if (mask[i] & groupbit) { + 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; - update_x_and_omega(x[i],v[i],omega[i],f[i],torque[i],d3rot); + } 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); + + // 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 - quat = bonus[ellipsoid[i]].quat; + 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)); - MathExtra::vecquat(omega[i],quat,wq); + } 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 ); - 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); - - MathExtra::quat_to_mat( quat, Q ); - MathExtra::matvec( Q, f_act, 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); -} diff --git a/src/USER-BROWNIAN/fix_brownian_asphere.h b/src/USER-BROWNIAN/fix_brownian_asphere.h index 31ccf4fc8b..c229ff1c8b 100644 --- a/src/USER-BROWNIAN/fix_brownian_asphere.h +++ b/src/USER-BROWNIAN/fix_brownian_asphere.h @@ -20,44 +20,29 @@ 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(); + + + + }; } diff --git a/src/USER-BROWNIAN/fix_brownian_base.cpp b/src/USER-BROWNIAN/fix_brownian_base.cpp new file mode 100644 index 0000000000..1a571c7a74 --- /dev/null +++ b/src/USER-BROWNIAN/fix_brownian_base.cpp @@ -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 +#include +#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; + +} diff --git a/src/USER-BROWNIAN/fix_brownian_base.h b/src/USER-BROWNIAN/fix_brownian_base.h new file mode 100644 index 0000000000..356742f314 --- /dev/null +++ b/src/USER-BROWNIAN/fix_brownian_base.h @@ -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. + +*/ diff --git a/src/USER-BROWNIAN/fix_brownian_sphere.cpp b/src/USER-BROWNIAN/fix_brownian_sphere.cpp index c3e0cd2457..a01eaac719 100644 --- a/src/USER-BROWNIAN/fix_brownian_sphere.cpp +++ b/src/USER-BROWNIAN/fix_brownian_sphere.cpp @@ -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 (!gamma_t_flag || !gamma_r_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"); - - // initialize Marsaglia RNG with processor-unique seed - random = new RanMars(lmp,seed + comm->me); + if (!atom->mu_flag) + error->all(FLERR,"Fix brownian/sphere requires atom attribute mu"); + } -/* ---------------------------------------------------------------------- */ - -int FixBrownianSphere::setmask() -{ - int mask = 0; - mask |= INITIAL_INTEGRATE; - mask |= POST_FORCE; - return mask; -} /* ---------------------------------------------------------------------- */ FixBrownianSphere::~FixBrownianSphere() { - delete random; } @@ -140,182 +71,162 @@ FixBrownianSphere::~FixBrownianSphere() void FixBrownianSphere::init() { + + FixBrownianBase::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; - } + 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; - int d3rot; // whether to compute angular momentum in xy plane + if (domain->dimension == 2) { - if (domain->dimension==2) { - d3rot = 0; - } else { - d3rot = 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 (!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>(); } - } - - 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; - } + } 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; } -/* ---------------------------------------------------------------------- - 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 (vflag) v_setup(vflag); - else evflag = 0; - - double fx,fy,fz; - double v[6]; - - for (int i = 0; i < nlocal; i++) + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + + double dx,dy,dz; + + 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); + 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; + + } 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; + + } + } + + x[i][0] += dx; + v[i][0] = dx/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; + } -} + } -void FixBrownianSphere::reset_dt() -{ - - dt = update->dt; - sqrtdt = sqrt(dt); + return; } diff --git a/src/USER-BROWNIAN/fix_brownian_sphere.h b/src/USER-BROWNIAN/fix_brownian_sphere.h index 8d30001fdd..e518d28db3 100644 --- a/src/USER-BROWNIAN/fix_brownian_sphere.h +++ b/src/USER-BROWNIAN/fix_brownian_sphere.h @@ -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; }; } diff --git a/src/random_mars.cpp b/src/random_mars.cpp index f31a92624a..2844bde417 100644 --- a/src/random_mars.cpp +++ b/src/random_mars.cpp @@ -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 diff --git a/src/random_mars.h b/src/random_mars.h index 7028ef54d2..1bcd16b051 100644 --- a/src/random_mars.h +++ b/src/random_mars.h @@ -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);