diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 1be5348bfe..8be0908a1a 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -40,6 +40,7 @@ OPT. * :doc:`aveforce ` * :doc:`balance ` * :doc:`bd/sphere ` + * :doc:`bd/asphere ` * :doc:`bocs ` * :doc:`bond/break ` * :doc:`bond/create ` diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 0993a885a5..c395489a58 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -181,8 +181,9 @@ accelerated styles exist. * :doc:`ave/histo/weight ` - weighted version of fix ave/histo * :doc:`ave/time ` - compute/output global time-averaged quantities * :doc:`aveforce ` - add an averaged force to each atom -* :doc:`bd/sphere ` - overdamped translational and rotational brownian dynamics * :doc:`balance ` - perform dynamic load-balancing +* :doc:`bd/asphere ` - integrate positions and orientations in overdamped motion +* :doc:`bd/sphere ` - overdamped translational and rotational brownian dynamics * :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_bd_asphere.rst b/doc/src/fix_bd_asphere.rst new file mode 100644 index 0000000000..8c0337cb5e --- /dev/null +++ b/doc/src/fix_bd_asphere.rst @@ -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 ` 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) `), :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 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) ` 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-MISC 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 bd/sphere `, :doc:`fix langevin `, +:doc:`fix nve/asphere `, :doc:`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). + + diff --git a/examples/USER/misc/bd_asphere/README.txt b/examples/USER/misc/bd_asphere/README.txt new file mode 100644 index 0000000000..7cb5d2e8e6 --- /dev/null +++ b/examples/USER/misc/bd_asphere/README.txt @@ -0,0 +1,2 @@ +The input file in2d.bd demonstrates how to run a 2d simulation +of ellipsoidal particles undergoing overdamped brownian motion. diff --git a/examples/USER/misc/bd_asphere/in3d.bd b/examples/USER/misc/bd_asphere/in3d.bd new file mode 100644 index 0000000000..6163eb0c58 --- /dev/null +++ b/examples/USER/misc/bd_asphere/in3d.bd @@ -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 diff --git a/examples/USER/misc/bd_asphere/log_gaussian_4_1_7_13_3d.lammps.log b/examples/USER/misc/bd_asphere/log_gaussian_4_1_7_13_3d.lammps.log new file mode 100644 index 0000000000..c6926dbccf --- /dev/null +++ b/examples/USER/misc/bd_asphere/log_gaussian_4_1_7_13_3d.lammps.log @@ -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 diff --git a/src/USER-MISC/README b/src/USER-MISC/README index 2dd16d23ac..ac41d60dfd 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -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 diff --git a/src/USER-MISC/fix_bd_asphere.cpp b/src/USER-MISC/fix_bd_asphere.cpp new file mode 100644 index 0000000000..88772b9237 --- /dev/null +++ b/src/USER-MISC/fix_bd_asphere.cpp @@ -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 +#include +#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); +} diff --git a/src/USER-MISC/fix_bd_asphere.h b/src/USER-MISC/fix_bd_asphere.h new file mode 100644 index 0000000000..696b659729 --- /dev/null +++ b/src/USER-MISC/fix_bd_asphere.h @@ -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. + +*/