Added bd integrator for ellipsoidal particles as well.
This commit is contained in:
Sam Cameron
@ -40,6 +40,7 @@ OPT.
|
||||
* :doc:`aveforce <fix_aveforce>`
|
||||
* :doc:`balance <fix_balance>`
|
||||
* :doc:`bd/sphere <fix_bd_sphere>`
|
||||
* :doc:`bd/asphere <fix_bd_asphere>`
|
||||
* :doc:`bocs <fix_bocs>`
|
||||
* :doc:`bond/break <fix_bond_break>`
|
||||
* :doc:`bond/create <fix_bond_create>`
|
||||
|
||||
@ -181,8 +181,9 @@ accelerated styles exist.
|
||||
* :doc:`ave/histo/weight <fix_ave_histo>` - weighted version of fix ave/histo
|
||||
* :doc:`ave/time <fix_ave_time>` - compute/output global time-averaged quantities
|
||||
* :doc:`aveforce <fix_aveforce>` - add an averaged force to each atom
|
||||
* :doc:`bd/sphere <fix_bd_sphere>` - overdamped translational and rotational brownian dynamics
|
||||
* :doc:`balance <fix_balance>` - perform dynamic load-balancing
|
||||
* :doc:`bd/asphere <fix_bd_asphere>` - integrate positions and orientations in overdamped motion
|
||||
* :doc:`bd/sphere <fix_bd_sphere>` - overdamped translational and rotational brownian dynamics
|
||||
* :doc:`bocs <fix_bocs>` - NPT style time integration with pressure correction
|
||||
* :doc:`bond/break <fix_bond_break>` - break bonds on the fly
|
||||
* :doc:`bond/create <fix_bond_create>` - create bonds on the fly
|
||||
|
||||
149
doc/src/fix_bd_asphere.rst
Normal file
149
doc/src/fix_bd_asphere.rst
Normal file
@ -0,0 +1,149 @@
|
||||
.. index:: fix bd/asphere
|
||||
|
||||
fix bd/asphere command
|
||||
======================
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. parsed-literal::
|
||||
|
||||
fix ID group-ID bd/asphere gamma_t gamma_r diff_t diff_r seed keyword args
|
||||
|
||||
* ID, group-ID are documented in :doc:`fix <fix>` command
|
||||
* bd/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 bd/asphere 1.0 1.0 1.0 1.0 1294019
|
||||
fix 1 all bd/asphere 1.0 1.0 1.0 1.0 19581092 rng none dipole
|
||||
fix 1 all bd/asphere 1.0 1.0 1.0 1.0 19581092 rng uniform
|
||||
fix 1 all bd/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) <GoldsteinCM1>`), :math:`dW_t` and
|
||||
:math:`dW_r` are Wiener processes (see e.g. :ref:`(Gardiner) <GardinerC1>`).
|
||||
The quaternions :math:`q` of the ellipsoid are updated each timestep from
|
||||
the angular velocity vector.
|
||||
|
||||
See :doc:`fix bd/sphere <fix_bd_sphere>` for discussion on the
|
||||
values of :math:`\gamma_t`, :math:`\gamma_r`, :math:`D_t`,
|
||||
and :math:`D_r` 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) <Dunweg7>` for why this works). This is the same method
|
||||
of noise generation as used in :doc:`fix_langevin <fix_langevin>`.
|
||||
|
||||
If the *rng* keyword is used with the *gaussian* value, then the noise
|
||||
is generated from a gaussian distribution. Typically this added
|
||||
complexity is unnecessary, and one should be fine using the *uniform*
|
||||
value for reasons argued in :ref:`(Dunweg) <Dunweg7>`.
|
||||
|
||||
If the *rng* keyword is used with the *none* value, then the noise
|
||||
terms are set to zero.
|
||||
|
||||
If the *dipole* keyword is used, then the dipole moments of the particles
|
||||
are updated by setting them along the x axis of the ellipsoidal frames of
|
||||
reference.
|
||||
|
||||
----------
|
||||
|
||||
.. include:: accel_styles.rst
|
||||
|
||||
----------
|
||||
|
||||
Restart, fix_modify, output, run start/stop, minimize info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
No information about this fix is written to :doc:`binary restart files <restart>`.
|
||||
No global or per-atom quantities are stored
|
||||
by this fix for access by various :doc:`output commands <Howto_output>`.
|
||||
|
||||
The :doc:`fix_modify <fix_modify>` *virial* option is supported by this
|
||||
fix to add the contribution due to the added forces on atoms to the
|
||||
system's virial as part of :doc:`thermodynamic output <thermo_style>`.
|
||||
The default is *virial no*.
|
||||
|
||||
No parameter of this fix can be used with the *start/stop* keywords of
|
||||
the :doc:`run <run>` command. This fix is not invoked during
|
||||
:doc:`energy minimization <minimize>`.
|
||||
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This fix requires that atoms store torque and angular velocity (omega)
|
||||
as defined by the :doc:`atom_style sphere <atom_style>` command, as well
|
||||
as atoms which have a definite orientation as defined by the
|
||||
:doc:`atom_style ellipsoid <atom_style>` command.
|
||||
Optionally, they can also store a dipole moment as defined by the
|
||||
:doc:`atom_style dipole <atom_style>` command.
|
||||
|
||||
This fix is part of the USER-MISC package. It is only enabled if
|
||||
LAMMPS was built with that package. See the :doc:`Build package <Build_package>`
|
||||
doc page for more info.
|
||||
|
||||
All particles in the group must be finite-size ellipsoids. They cannot
|
||||
be point particles.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
:doc:`fix bd/sphere <fix_bd_sphere>`, :doc:`fix langevin <fix_langevin>`,
|
||||
:doc:`fix nve/asphere <fix_nve_asphere>`, :doc:`atom style <atom_style>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
The default for *rng* is *uniform*.
|
||||
|
||||
----------
|
||||
|
||||
.. _GoldsteinCM1:
|
||||
|
||||
**(Goldstein)** Goldstein, Poole, and Safko, Classical Mechanics, 3rd Ed. (2001).
|
||||
|
||||
.. _GardinerC1:
|
||||
|
||||
**(Gardiner)** Gardiner, A Handbook for the Natural and Social Sciences 4th Ed. (2009).
|
||||
|
||||
.. _Dunweg7:
|
||||
|
||||
**(Dunweg)** Dunweg and Paul, Int J of Modern Physics C, 2, 817-27 (1991).
|
||||
|
||||
|
||||
2
examples/USER/misc/bd_asphere/README.txt
Normal file
2
examples/USER/misc/bd_asphere/README.txt
Normal file
@ -0,0 +1,2 @@
|
||||
The input file in2d.bd demonstrates how to run a 2d simulation
|
||||
of ellipsoidal particles undergoing overdamped brownian motion.
|
||||
71
examples/USER/misc/bd_asphere/in3d.bd
Normal file
71
examples/USER/misc/bd_asphere/in3d.bd
Normal file
@ -0,0 +1,71 @@
|
||||
# 3d overdamped brownian dynamics
|
||||
|
||||
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}_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
|
||||
|
||||
|
||||
neighbor 1.0 bin
|
||||
neigh_modify every 1 delay 1 check yes
|
||||
|
||||
|
||||
|
||||
pair_style none
|
||||
|
||||
|
||||
fix 1 all bd/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
|
||||
|
||||
|
||||
#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
|
||||
@ -0,0 +1,154 @@
|
||||
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 bd/asphere ${gamma_t} ${gamma_r} ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all bd/asphere 4 ${gamma_r} ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all bd/asphere 4 1 ${D_t} ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all bd/asphere 4 1 7 ${D_r} ${seed} rng ${rng} dipole
|
||||
fix 1 all bd/asphere 4 1 7 13 ${seed} rng ${rng} dipole
|
||||
fix 1 all bd/asphere 4 1 7 13 1974019 rng ${rng} dipole
|
||||
fix 1 all bd/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.68636 on 1 procs for 50000 steps with 512 atoms
|
||||
|
||||
Performance: 0.076 tau/day, 8792.977 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.11589 | 0.11589 | 0.11589 | 0.0 | 2.04
|
||||
Output | 2.1155e-05 | 2.1155e-05 | 2.1155e-05 | 0.0 | 0.00
|
||||
Modify | 5.4911 | 5.4911 | 5.4911 | 0.0 | 96.57
|
||||
Other | | 0.07936 | | | 1.40
|
||||
|
||||
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}_2d.lammpstrj id type x y xu yu mux muy muz fx fy fz
|
||||
dump 1 all custom 1000 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 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) = 7.103 | 7.103 | 7.103 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.595 on 1 procs for 120000 steps with 512 atoms
|
||||
|
||||
Performance: 7626.329 tau/day, 8826.770 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.00079148 | 0.00079148 | 0.00079148 | 0.0 | 0.01
|
||||
Comm | 0.029132 | 0.029132 | 0.029132 | 0.0 | 0.21
|
||||
Output | 0.15178 | 0.15178 | 0.15178 | 0.0 | 1.12
|
||||
Modify | 13.222 | 13.222 | 13.222 | 0.0 | 97.25
|
||||
Other | | 0.1916 | | | 1.41
|
||||
|
||||
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:19
|
||||
@ -53,7 +53,8 @@ dihedral_style table/cut, Mike Salerno, ksalerno@pha.jhu.edu, 11 May 18
|
||||
fix accelerate/cos, Zheng Gong (ENS de Lyon), z.gong@outlook.com, 24 Apr 20
|
||||
fix addtorque, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11
|
||||
fix ave/correlate/long, Jorge Ramirez (UPM Madrid), jorge.ramirez at upm.es, 21 Oct 2015
|
||||
fix bd/sphere, Sam Cameron (U of Bristol), samuel.j.m.cameron at gmail.com, 10 Dec 2020
|
||||
fix bd/sphere, Sam Cameron (U of Bristol), samuel.j.m.cameron at gmail.com, 13 Dec 2020
|
||||
fix bd/asphere, Sam Cameron (U of Bristol), samuel.j.m.cameron at gmail.com, 13 Dec 2020
|
||||
fix electron/stopping/fit, James Stewart (SNL), jstewa .at. sandia.gov, 23 Sep 2020
|
||||
fix electron/stopping, Konstantin Avchaciov, k.avchachov at gmail.com, 26 Feb 2019
|
||||
fix ffl, David Wilkins (EPFL Lausanne), david.wilkins @ epfl.ch, 28 Sep 2018
|
||||
|
||||
392
src/USER-MISC/fix_bd_asphere.cpp
Normal file
392
src/USER-MISC/fix_bd_asphere.cpp
Normal file
@ -0,0 +1,392 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
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_bd_asphere.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
#include "math_extra.h"
|
||||
#include "atom.h"
|
||||
#include "atom_vec_ellipsoid.h"
|
||||
#include "force.h"
|
||||
#include "update.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "random_mars.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
#define SMALL 1e-14
|
||||
|
||||
enum{NODIPOLE,DIPOLE};
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBdAsphere::FixBdAsphere(LAMMPS *lmp, int narg, char **arg) :
|
||||
Fix(lmp, narg, arg)
|
||||
{
|
||||
virial_flag = 1;
|
||||
|
||||
time_integrate = 1;
|
||||
|
||||
dipole_flag = NODIPOLE;
|
||||
|
||||
|
||||
if (narg > 11 || narg < 8 )
|
||||
error->all(FLERR,"Illegal fix bd/asphere command.");
|
||||
|
||||
if (!atom->sphere_flag)
|
||||
error->all(FLERR,"Fix bd/asphere requires atom style sphere");
|
||||
|
||||
gamma_t = utils::numeric(FLERR,arg[3],false,lmp);
|
||||
if (gamma_t <= 0.0)
|
||||
error->all(FLERR,"Fix bd/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 bd/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 bd/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 bd/asphere rotational diffusion "
|
||||
"coefficient must be > 0.");
|
||||
|
||||
seed = utils::inumeric(FLERR,arg[7],false,lmp);
|
||||
if (seed <= 0) error->all(FLERR,"Fix bd/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/bd/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/bd/asphere command.");
|
||||
}
|
||||
iarg = iarg + 2;
|
||||
} else if (strcmp(arg[iarg],"dipole") == 0) {
|
||||
dipole_flag = DIPOLE;
|
||||
iarg = iarg + 1;
|
||||
} else {
|
||||
error->all(FLERR,"Illegal fix/bd/asphere command.");
|
||||
}
|
||||
}
|
||||
|
||||
if (dipole_flag == DIPOLE && !atom->mu_flag)
|
||||
error->all(FLERR,"Fix bd/asphere dipole requires atom attribute mu");
|
||||
|
||||
// initialize Marsaglia RNG with processor-unique seed
|
||||
random = new RanMars(lmp,seed + comm->me);
|
||||
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
int FixBdAsphere::setmask()
|
||||
{
|
||||
int mask = 0;
|
||||
mask |= INITIAL_INTEGRATE;
|
||||
mask |= POST_FORCE;
|
||||
return mask;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
FixBdAsphere::~FixBdAsphere()
|
||||
{
|
||||
delete random;
|
||||
}
|
||||
|
||||
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixBdAsphere::init()
|
||||
{
|
||||
|
||||
|
||||
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
|
||||
if (!avec)
|
||||
error->all(FLERR,"Compute bd/asphere requires "
|
||||
"atom style ellipsoid");
|
||||
|
||||
// check that all particles are finite-size ellipsoids
|
||||
// no point particles allowed, spherical is OK
|
||||
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
for (int i = 0; i < nlocal; i++)
|
||||
if (mask[i] & groupbit)
|
||||
if (ellipsoid[i] < 0)
|
||||
error->one(FLERR,"Fix bd/asphere requires extended particles");
|
||||
|
||||
|
||||
if (dipole_flag == DIPOLE) {
|
||||
|
||||
double f_act[3] = { 1.0, 0.0, 0.0 };
|
||||
double f_rot[3];
|
||||
double *quat;
|
||||
int *ellipsoid = atom->ellipsoid;
|
||||
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
|
||||
|
||||
double Q[3][3];
|
||||
|
||||
double **mu = atom->mu;
|
||||
|
||||
for (int i = 0; i < nlocal; i++) {
|
||||
if (mask[i] & groupbit) {
|
||||
quat = bonus[ellipsoid[i]].quat;
|
||||
MathExtra::quat_to_mat( quat, Q );
|
||||
MathExtra::matvec( Q, f_act, f_rot );
|
||||
|
||||
mu[i][0] = f_rot[0];
|
||||
mu[i][1] = f_rot[1];
|
||||
mu[i][2] = f_rot[2];
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
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;
|
||||
}
|
||||
|
||||
dt = update->dt;
|
||||
sqrtdt = sqrt(dt);
|
||||
}
|
||||
|
||||
void FixBdAsphere::setup(int vflag)
|
||||
{
|
||||
post_force(vflag);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
|
||||
void FixBdAsphere::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;
|
||||
|
||||
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
|
||||
|
||||
int d3rot; // whether to compute angular momentum in xy plane
|
||||
|
||||
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 *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 Q[3][3];
|
||||
|
||||
double **mu = atom->mu;
|
||||
|
||||
|
||||
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);
|
||||
|
||||
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 FixBdAsphere::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 FixBdAsphere::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 FixBdAsphere::reset_dt()
|
||||
{
|
||||
|
||||
dt = update->dt;
|
||||
sqrtdt = sqrt(dt);
|
||||
}
|
||||
103
src/USER-MISC/fix_bd_asphere.h
Normal file
103
src/USER-MISC/fix_bd_asphere.h
Normal file
@ -0,0 +1,103 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef FIX_CLASS
|
||||
|
||||
FixStyle(bd/asphere,FixBdAsphere)
|
||||
|
||||
#else
|
||||
|
||||
#ifndef LMP_FIX_BD_ASPHERE_H
|
||||
#define LMP_FIX_BD_ASPHERE_H
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class FixBdAsphere : public Fix {
|
||||
public:
|
||||
FixBdAsphere(class LAMMPS *, int, char **);
|
||||
virtual ~FixBdAsphere();
|
||||
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 );
|
||||
|
||||
|
||||
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
|
||||
|
||||
class AtomVecEllipsoid *avec;
|
||||
|
||||
protected:
|
||||
class RanMars *random;
|
||||
typedef double (RanMars::*rng_member)();
|
||||
rng_member rng_func; // placeholder for RNG function
|
||||
|
||||
};
|
||||
|
||||
}
|
||||
#endif
|
||||
#endif
|
||||
|
||||
/* ERROR/WARNING messages:
|
||||
|
||||
E: Illegal fix bd/asphere command.
|
||||
|
||||
Wrong number/type of input arguments.
|
||||
|
||||
E: Compute bd/asphere requires atom style sphere
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Compute bd/asphere requires atom style ellipsoid
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Compute bd/asphere dipole requires atom attribute mu
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix bd/asphere translational viscous drag coefficient must be > 0.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix bd/asphere rotational viscous drag coefficient must be > 0.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix bd/asphere translational diffusion coefficient must be > 0.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix bd/asphere rotational diffusion coefficient must be > 0.
|
||||
|
||||
Self-explanatory.
|
||||
|
||||
E: Fix bd/asphere seed must be > 0.
|
||||
|
||||
*/
|
||||
Reference in New Issue
Block a user