diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index a778c89d42..8dd5da17db 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -145,6 +145,7 @@ set(STANDARD_PACKAGES AWPMD BOCS BODY + BPM BROWNIAN CG-DNA CG-SDK diff --git a/doc/src/Commands_bond.rst b/doc/src/Commands_bond.rst index 1a0876e88f..b79e8ee174 100644 --- a/doc/src/Commands_bond.rst +++ b/doc/src/Commands_bond.rst @@ -32,6 +32,7 @@ OPT. * * * + * :doc:`bpm/rotational ` * :doc:`class2 (ko) ` * :doc:`fene (iko) ` * :doc:`fene/expand (o) ` diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index 9dfb28fa8b..4254c4f65c 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -89,6 +89,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`msd ` * :doc:`msd/chunk ` * :doc:`msd/nongauss ` + * :doc:`nbond/atom ` * :doc:`omega/chunk ` * :doc:`orientorder/atom (k) ` * :doc:`pair ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 45a75ff394..ef9c666751 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -140,6 +140,7 @@ OPT. * :doc:`nve/manifold/rattle ` * :doc:`nve/noforce ` * :doc:`nve/sphere (ko) ` + * :doc:`nve/sphere/bpm ` * :doc:`nve/spin ` * :doc:`nve/tri ` * :doc:`nvk ` @@ -221,6 +222,7 @@ OPT. * :doc:`spring/rg ` * :doc:`spring/self ` * :doc:`srd ` + * :doc:`store/local ` * :doc:`store/force ` * :doc:`store/state ` * :doc:`tdpd/source ` @@ -238,6 +240,7 @@ OPT. * :doc:`ttm ` * :doc:`ttm/mod ` * :doc:`tune/kspace ` + * :doc:`update/special/bonds ` * :doc:`vector ` * :doc:`viscosity ` * :doc:`viscous ` diff --git a/doc/src/Examples.rst b/doc/src/Examples.rst index 649be52ab7..159de8aea5 100644 --- a/doc/src/Examples.rst +++ b/doc/src/Examples.rst @@ -54,6 +54,8 @@ Lowercase directories +-------------+------------------------------------------------------------------+ | body | body particles, 2d system | +-------------+------------------------------------------------------------------+ +| bpm | bonded particle models of pouring, crushing, and fracture | ++-------------+------------------------------------------------------------------+ | cmap | CMAP 5-body contributions to CHARMM force field | +-------------+------------------------------------------------------------------+ | colloid | big colloid particles in a small particle solvent, 2d system | diff --git a/doc/src/Howto_bpm.rst b/doc/src/Howto_bpm.rst new file mode 100755 index 0000000000..24dc77d44b --- /dev/null +++ b/doc/src/Howto_bpm.rst @@ -0,0 +1,77 @@ +Bonded particle models +=============== + +Bonded particle models are used to simulate mesoscale solids. +Solids are constructed as a collection of particles which each +represent a coarse-grained region of space much larger than the +atomistic scale. Particles within a solid region are then connected +by a network of bonds to provide solid elasticity. + +Unlike traditional bonds in molecular dynamics, the equilibrium +bond length can vary between bonds. Bonds store the reference state. +This includes setting the equilibrium length equal to the initial +distance between the two particles but can also include data on the +bond orientation for rotational models. This produces a stress free +initial state. Furthermore, bonds are allowed to break under large +strains producing fracture. + +Bonds can be created using a :doc:`read data ` +or :doc:`create bond ` command. Alternatively, a +:doc:`molecule ` template with bonds can be used with +:doc:`fix deposit ` or :doc:`fix pour ` to +create solid grains. +In this implementation, bonds store their reference state when they +are first computed in the setup of a simulation run. Data is then +preserved across run commands and is written to :doc:`binary restart files ` +such that restarting the system will not reset the reference state of a bond. + +As bonds can be broken between neighbor list builds, :doc:`special_bonds ` +work differently for BPM bond styles. There are two possible special +bond settings which determine how pair interactions work between bonded +particles. First, one can simply overlay pair interactions such that all +bonded particles also feel pair interactions. This can be accomplished by +simply turning off all special bonds by setting + +.. code-block:: LAMMPS + + special_bonds lj/coul 1 1 1 + +Alternatively, one can censor all pair interactions between bonded particles. +Unlike :doc:`bond quartic `, this is not done by subtracting +pair forces during the bond computation but rather by dynamically updating +the special bond list. To do this, one must both define an instance of +:doc:`fix update/special/bonds ` and have the special bond +settings + +.. code-block:: LAMMPS + + special_bonds lj 0 1 1 coul 1 1 1 + +This fix ensures the 1-2 special bond list remains updated as bonds break. The fix +also requires :doc:`newton ` bond off such that whena bond breaks between +atoms across multiple processors, all processors are aware of the event. +The special bond settings then accomplish two tasks. First, they turns off 1-3 and +1-4 special bond lists, which are not currently supported for BPMs. As BPMs often +have dense bond networks, generating 1-3 and 1-4 special bond lists is expensive. +By setting the lj weight for 1-2 bonds to zero, this censors pairwise interactions. +However, setting a nonzero coul weight for 1-2 bonds ensures all bonded +neighbors are included in the neighbor list. All bonded neighbors must be included +in neighbor lists as they could become unbonded at any timestep. + +Currently there are two types of bonds included in this package. The first +bond style, :doc:`bond bpm/simple `, only applies pairwise, +central body forces. Point particles must have :doc:`bond atom style ` +and may be thought of as nodes in a spring network. Alternatively, +the second bond style, :doc:`bond bpm/rotational `, +resolves tangential forces and torques arising with the shearing, bending, +and twisting of the bond due to rotation or displacement of particles. +Particles are similar to those used in the :doc:`granular package `, +:doc:`atom style sphere `. However, they must also track the +current orientation of particles and therefore use a derived :doc:`atom style sphere/bpm `. +This also requires a unique integrator :doc:`fix nve/sphere/bpm ` +which numerically integrates orientation similar to :doc:`fix nve/asphere `. + +To monitor the fracture of bonds in the system, all BPM bond styles +can be associated with an instance of :doc:`fix store/local ` +to record all instances of bond breakage for output. Additionally, one can use +:doc:`compute nbond/atom ` to tally the current number of bonds per atom. diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst index b4e50fad75..e788da552d 100644 --- a/doc/src/Packages_details.rst +++ b/doc/src/Packages_details.rst @@ -32,6 +32,7 @@ page gives those details. * :ref:`AWPMD ` * :ref:`BOCS ` * :ref:`BODY ` + * :ref:`BPM ` * :ref:`BROWNIAN ` * :ref:`CG-DNA ` * :ref:`CG-SDK ` @@ -284,6 +285,32 @@ overview. ---------- +.. _PKG-BPM: + +BPM package +------------ + +**Contents:** + +Pair styles, bond styles, fixes, and computes for bonded particle +models for mesoscale simulations of solids and fracture. See the +:doc:`Howto bpm ` page for an overview. + +**Authors:** Joel T. Clemmer (Sandia National Labs) + +**Supporting info:** + +* src/BPM filenames -> commands +* :doc:`Howto_bpm ` +* :doc:`atom_style sphere/bpm ` +* :doc:`bond_style bpm/rotational ` +* :doc:`compute nbonds/atom ` +* :doc:`fix nve/sphere/bpm ` +* :doc:`fix update/special/bonds ` +* examples/bpm + +---------- + .. _PKG-BROWNIAN: BROWNIAN package @@ -1564,7 +1591,6 @@ listing, "ls src/MISC", to see the list of commands. * :doc:`pair_style list ` * :doc:`pair_style srp ` * :doc:`pair_style tracker ` -* :doc:`fix pair/tracker ` ---------- diff --git a/doc/src/Packages_list.rst b/doc/src/Packages_list.rst index 22dfa2c69e..4ab3f179e7 100644 --- a/doc/src/Packages_list.rst +++ b/doc/src/Packages_list.rst @@ -58,6 +58,11 @@ whether an extra library is needed to build and use the package: - :doc:`Howto body ` - body - no + * - :ref:`BPM ` + - bonded particle models + - :doc:`Howto bpm ` + - bpm + - no * - :ref:`BROWNIAN ` - Brownian dynamics, self-propelled particles - :doc:`fix brownian `, :doc:`fix propel/self ` diff --git a/doc/src/atom_style.rst b/doc/src/atom_style.rst index bade8c2f79..09eca474ae 100644 --- a/doc/src/atom_style.rst +++ b/doc/src/atom_style.rst @@ -10,7 +10,7 @@ Syntax atom_style style args -* style = *angle* or *atomic* or *body* or *bond* or *charge* or *dipole* or *dpd* or *edpd* or *electron* or *ellipsoid* or *full* or *line* or *mdpd* or *molecular* or *oxdna* or *peri* or *smd* or *sph* or *sphere* or *spin* or *tdpd* or *tri* or *template* or *hybrid* +* style = *angle* or *atomic* or *body* or *bond* or *charge* or *dipole* or *dpd* or *edpd* or *electron* or *ellipsoid* or *full* or *line* or *mdpd* or *molecular* or *oxdna* or *peri* or *smd* or *sph* or *sphere* or *sphere/bpm* or *spin* or *tdpd* or *tri* or *template* or *hybrid* .. parsed-literal:: @@ -21,6 +21,7 @@ Syntax see the :doc:`Howto body ` doc page for details *sphere* arg = 0/1 (optional) for static/dynamic particle radii + *sphere/bpm* arg = 0/1 (optional) for static/dynamic particle radii *tdpd* arg = Nspecies Nspecies = # of chemical species *template* arg = template-ID @@ -120,6 +121,8 @@ quantities. +--------------+-----------------------------------------------------+--------------------------------------+ | *sphere* | diameter, mass, angular velocity | granular models | +--------------+-----------------------------------------------------+--------------------------------------+ +| *sphere/bpm* | diameter, mass, angular velocity, quaternion | granular bonded particle models (bpm)| ++--------------+-----------------------------------------------------+--------------------------------------+ | *spin* | magnetic moment | system with magnetic particles | +--------------+-----------------------------------------------------+--------------------------------------+ | *tdpd* | chemical concentration | tDPD particles | @@ -141,8 +144,9 @@ quantities. output the custom values. All of the above styles define point particles, except the *sphere*, -*ellipsoid*, *electron*, *peri*, *wavepacket*, *line*, *tri*, and -*body* styles, which define finite-size particles. See the :doc:`Howto spherical ` page for an overview of using +*sphere/bpm*, *ellipsoid*, *electron*, *peri*, *wavepacket*, *line*, +*tri*, and *body* styles, which define finite-size particles. See the +:doc:`Howto spherical ` page for an overview of using finite-size particle models with LAMMPS. All of the point-particle styles assign mass to particles on a @@ -150,15 +154,15 @@ per-type basis, using the :doc:`mass ` command, The finite-size particle styles assign mass to individual particles on a per-particle basis. -For the *sphere* style, the particles are spheres and each stores a +For the *sphere* and *sphere/bpm* styles, the particles are spheres and each stores a per-particle diameter and mass. If the diameter > 0.0, the particle is a finite-size sphere. If the diameter = 0.0, it is a point particle. Note that by use of the *disc* keyword with the :doc:`fix nve/sphere `, :doc:`fix nvt/sphere `, :doc:`fix nph/sphere `, :doc:`fix npt/sphere -` commands, spheres can be effectively treated as 2d +` commands for the *sphere* style, spheres can be effectively treated as 2d discs for a 2d simulation if desired. See also the :doc:`set -density/disc ` command. The *sphere* style takes an optional 0 +density/disc ` command. The *sphere* and *sphere/bpm* styles take an optional 0 or 1 argument. A value of 0 means the radius of each sphere is constant for the duration of the simulation. A value of 1 means the radii may vary dynamically during the simulation, e.g. due to use of @@ -195,6 +199,8 @@ position, which is represented by the eradius = electron size. For the *peri* style, the particles are spherical and each stores a per-particle mass and volume. +The *sphere/bpm* style is part of the BPM package. + The *oxdna* style is for coarse-grained nucleotides and stores the 3'-to-5' polarity of the nucleotide strand, which is set through the bond topology in the data file. The first (second) atom in a diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst new file mode 100644 index 0000000000..632e340da6 --- /dev/null +++ b/doc/src/bond_bpm_rotational.rst @@ -0,0 +1,197 @@ +.. index:: bond_style bpm/rotational + +bond_style bpm/rotational command +========================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + bond_style bpm/rotational keyword value attribute1 attribute2 ... + +* optional keyword = *store/local* + + .. parsed-literal:: + + *store/local* values = ID of associated fix store/local followed by one or more attributes + + *id1, id2* = IDs of 2 atoms in the bond + *time* = the time the bond broke + *x, y, z* = the center of mass position of the 2 atoms when the bond broke + *x/ref, y/ref, z/ref* = the inintial center of mass position of the 2 atoms + + +Examples +"""""""" + +.. code-block:: LAMMPS + + bond_style bpm/rotational + bond_coeff 1 + + bond_style bpm/rotational myfix time id1 id2 + fix myfix all store/local 1000 3 + dump 1 all local 1000 dump.broken f_myfix[1] f_myfix[2] f_myfix[3] + dump_modify 1 write_header no + +Description +""""""""""" + +The *bpm/rotational* bond style computes forces and torques based +on deviations from the initial reference state of the two atoms. +The reference state is stored by each bond when it is first computed +in the setup of a run. Data is then preserved across run commands and +is written to :doc:`binary restart files ` such that restarting +the system will not reset the reference state of a bond. + +Forces include a normal and tangential component. The base normal force +has a magnitude of + +.. math:: + + F_r = K_r (r - r_0) + +where :math:`K_r` is a stiffness and :math:`r` is the current distance and +:math:`r_0` is the initial distance between the two particles. + +A tangential force is applied perpendicular to the normal direction +which is proportional to the tangential shear displacement with a stiffness +of :math:`K_s`. This tangential force also induces a torque. +In addition, bending and twisting torques are also applied to particles +which are proportional to angular bending and twisting displacements with +stiffnesses of :math`K_b` and :math:`K_t', respectively. +Details on the calculations of shear displacements and angular displacements +can be found in :ref:`(Wang) ` and :ref:`(Wang and Mora) `. + +Bonds will break under sufficient stress. A breaking criteria is calculated + +.. math:: + + B = \alpha(r, r_0) \frac{F_r}{F_{r,c}} + \frac{F_s}{F_{s,c}} + + \frac{\tau_b}{\tau_{b,c}} + \frac{\tau_t}{\tau_{t,c}} + +where :math:`F_s` is the magnitude of the shear force and +:math:`\tau_b` and :math:`\tau_t` are the magnitudes of the bending and +twisting forces, respectively. The corresponding variables :math:`F_{r,c}` +:math:`F_{s,c}`, :math:`\tau_{b,c}`, and :math:`\tau_{t,c}` are critical +limits to each force or torque. The term :math:`\alpha` is simply one in +extension and zero in compression such that the normal force component +does not contribute to the breaking criteria in compression. +If :math:`B` is ever equal to or exceeds one, the bond will break. +This is done by setting by setting its type to 0 such that forces and +torques are no longer computed. + +After computing the base magnitudes of the forces and torques, they are +all multiplied by an extra factor :math:`w` to smoothly interpolate +forces and torques to zero as the bond breaks. This term is calculated +as :math:`w = (1.0 - B)^2`. + +Finally, additional damping forces and torques are applied to the two +particles. A force is applied proportional to the difference in the +normal velocity of particles using a similar construction as +dissipative particle dynamics (:ref:`(Groot) `): + +.. math:: + + F_D = - \gamma_n w (\hat{r} \bullet \vec{v}) + +where :math:`\gamma_n` is the damping strength, :math:`\hat{r}` is the +radial normal vector, and :math:`\vec{v}` is the velocity difference +between the two particles. Similarly, tangential forces are applied to +each atom proportional to the relative differences in sliding velocities +with a constant prefactor :math:`\gamma_s` (:ref:`(Wang et al.) ) +along with their associated torques. The rolling and twisting components of +the relative angular velocities of the two atoms are also damped by applying +torques with prefactors of :math:`\gamma_r` and :math:`\gamma_t`, respectively. + +The following coefficients must be defined for each bond type via the +:doc:`bond_coeff ` command as in the example above, or in +the data file or restart files read by the :doc:`read_data ` +or :doc:`read_restart ` commands: + +* :math:`K_r` (force/distance units) +* :math:`K_s` (force/distance units) +* :math:`K_t` (force units) +* :math:`K_b` (force units) +* :math:`F_{r,c}` (force units) +* :math:`F_{s,c}` (force units) +* :math:`\tau_{b,c}` (force*distance units) +* :math:`\tau_{t,c}` (force*distance units) +* :math:`\gamma_n` (force/velocity units) +* :math:`\gamma_s` (force/velocity units) +* :math:`\gamma_r` (distance*force/seconds/radians units) +* :math:`\gamma_t` (distance*force/seconds/radians units) + +As bonds can be broken between neighbor list builds, particular +:doc:`special_bonds ` are required. See the `:doc: how to ` +page on BPMs or `:doc: fix update/special/bonds ` +for details. + +This bond style tracks broken bonds and can record them using an instance of +:doc:`fix store/local ` if the *store/local* keyword is +used followed by the ID of the fix and then a series of bond attributes. + +Note that when bonds are dumped to a file via the :doc:`dump local ` +command, bonds with type 0 (broken bonds) are not included. The +:doc:`delete_bonds ` command can also be used to query the +status of broken bonds or permanently delete them, e.g.: + +.. code-block:: LAMMPS + + delete_bonds all stats + delete_bonds all bond 0 remove + + +---------- + +Restart +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +This bond style writes the reference state of each bond to +:doc:`binary restart files `. Loading a restart +file will properly resume bonds. + +Restrictions +"""""""""""" + +This bond style can only be used if LAMMPS was built with the BPM +package. See the :doc:`Build package ` doc page for more +info. + +The *bpm/rotational* style requires 1-3 and 1-4 :doc:`special_bonds ` +be turned off using the :doc:`special_bonds ` command. + +The *bpm/rotational* style requires :doc:`atom style sphere/bpm `. + +Related commands +"""""""""""""""" + +:doc:`bond_coeff `, :doc:`fix store/local `, +:doc:`fix nve/sphere/bpm ` + +Default +""""""" + +none + + +.. _Wang2009: + +**(Wang)** Wang, Acta Geotechnica, 4, +p 117-127 (2009). + +.. _Wang2009b: + +**(Wang and Mora)** Wang, Mora, Advances in Geocomputing, +119, p 183-228 (2009). + +.. _Groot1: + +**(Groot)** Groot and Warren, J Chem Phys, 107, 4423-35 (1997). + +.. _Wang2015: + +**(Wang et al, 2015)** Wang, Y., Alonso-Marroquin, F., & Guo, +W. W. (2015). Rolling and sliding in 3-D discrete element +models. Particuology, 23, 49-55. \ No newline at end of file diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 71d3bada76..0b8249cc7d 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -234,6 +234,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`msd ` - mean-squared displacement of group of atoms * :doc:`msd/chunk ` - mean-squared displacement for each chunk * :doc:`msd/nongauss ` - MSD and non-Gaussian parameter of group of atoms +* :doc:`nbond/atom ` - calculates number of bonds per atom * :doc:`omega/chunk ` - angular velocity for each chunk * :doc:`orientorder/atom ` - Steinhardt bond orientational order parameters Ql * :doc:`pair ` - values computed by a pair style diff --git a/doc/src/compute_nbonds_atom.rst b/doc/src/compute_nbonds_atom.rst new file mode 100644 index 0000000000..fd7ff4f8a9 --- /dev/null +++ b/doc/src/compute_nbonds_atom.rst @@ -0,0 +1,56 @@ +.. index:: compute nbonds/atom + +compute nbonds/atom command +======================= + +Syntax +"""""" + +.. parsed-literal:: + + compute ID group-ID nbonds/atom + +* ID, group-ID are documented in :doc:`compute ` command +* nbonds/atom = style name of this compute command + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all nbonds/atom + +Description +""""""""""" + +Define a computation that computes the number of bonds per-atom. +Bonds which are broken are not counted in the tally. +See :doc:`bond_style quartic ` or the +:doc:`Howto bpm ` page. The number of bonds will be zero +for atoms not in the specified compute group. + +Output info +""""""""""" + +This compute calculates a per-atom vector, which can be accessed by +any command that uses per-atom values from a compute as input. See +the :doc:`Howto output ` doc page for an overview of +LAMMPS output options. + +Restrictions +"""""""""""" + +This fix can only be used if LAMMPS was built with the BPM +package. See the :doc:`Build package ` doc page for more +info. + +Related commands +"""""""""""""""" + +Default +""""""" + +none + +---------- + diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 81c0c87320..9c749b2f38 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -283,6 +283,7 @@ accelerated styles exist. * :doc:`nve/manifold/rattle ` - * :doc:`nve/noforce ` - NVE without forces (v only) * :doc:`nve/sphere ` - NVE for spherical particles +* :doc:`nve/sphere/bpm ` - NVE for spherical particles used in bonded particle models * :doc:`nve/spin ` - NVE for a spin or spin-lattice system * :doc:`nve/tri ` - NVE for triangles * :doc:`nvk ` - constant kinetic energy time integration @@ -300,7 +301,6 @@ accelerated styles exist. * :doc:`orient/fcc ` - add grain boundary migration force for FCC * :doc:`orient/eco ` - add generalized grain boundary migration force * :doc:`pafi ` - constrained force averages on hyper-planes to compute free energies (PAFI) -* :doc:`pair/tracker ` - track properties of pairwise interactions * :doc:`phonon ` - calculate dynamical matrix from MD simulations * :doc:`pimd ` - Feynman path integral molecular dynamics * :doc:`planeforce ` - constrain atoms to move in a plane @@ -364,6 +364,7 @@ accelerated styles exist. * :doc:`spring/rg ` - spring on radius of gyration of group of atoms * :doc:`spring/self ` - spring from each atom to its origin * :doc:`srd ` - stochastic rotation dynamics (SRD) +* :doc:`store/local ` - store local data for output * :doc:`store/force ` - store force on each atom * :doc:`store/state ` - store attributes for each atom * :doc:`tdpd/source ` - @@ -381,6 +382,7 @@ accelerated styles exist. * :doc:`ttm ` - two-temperature model for electronic/atomic coupling * :doc:`ttm/mod ` - enhanced two-temperature model with additional options * :doc:`tune/kspace ` - auto-tune KSpace parameters +* :doc:`update/special/bonds ` - update special bond lists for BPM bond styles that allow for bond breakage * :doc:`vector ` - accumulate a global vector every N timesteps * :doc:`viscosity ` - Muller-Plathe momentum exchange for viscosity calculation * :doc:`viscous ` - viscous damping for granular simulations diff --git a/doc/src/fix_nve_sphere_bpm.rst b/doc/src/fix_nve_sphere_bpm.rst new file mode 100644 index 0000000000..1afa53a4de --- /dev/null +++ b/doc/src/fix_nve_sphere_bpm.rst @@ -0,0 +1,82 @@ +.. index:: fix nve/sphere/bpm + +fix nve/sphere/bpm command +====================== + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group-ID nve/sphere/bpm + +* ID, group-ID are documented in :doc:`fix ` command +* nve/sphere/bpm = style name of this fix command +* zero or more keyword/value pairs may be appended +* keyword = *disc* + + .. parsed-literal:: + + *disc* value = none = treat particles as 2d discs, not spheres + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all nve/sphere/bpm + fix 1 all nve/sphere/bpm disc + +Description +""""""""""" + +Perform constant NVE integration to update position, velocity, angular velocity, +and quaternion orientation for finite-size spherical particles in the group each +timestep. V is volume; E is energy. This creates a system trajectory +consistent with the microcanonical ensemble. + +This fix differs from the :doc:`fix nve ` command, which +assumes point particles and only updates their position and velocity. +It also differs from the :doc:`fix nve/sphere ` command, which +does not evaluate a particles orientation or quaternion. + +If the *disc* keyword is used, then each particle is treated as a 2d +disc (circle) instead of as a sphere. This is only possible for 2d +simulations, as defined by the :doc:`dimension ` keyword. +The only difference between discs and spheres in this context is their +moment of inertia, as used in the time integration. + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files `. None of the :doc:`fix_modify ` options +are relevant to this fix. No global or per-atom quantities are stored +by this fix for access by various :doc:`output commands `. +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, angular velocity (omega), +a radius, and a quaternion as defined by the :doc:`atom_style sphere/bpm ` +command. + +All particles in the group must be finite-size spheres with quaternions. They cannot +be point particles. + +Use of the *disc* keyword is only allowed for 2d simulations, as +defined by the :doc:`dimension ` keyword. + +Related commands +"""""""""""""""" + +:doc:`fix nve `, :doc:`fix nve/sphere ` + +Default +""""""" + +none + diff --git a/doc/src/fix_pair_tracker.rst b/doc/src/fix_pair_tracker.rst deleted file mode 100644 index 5c2ecf5774..0000000000 --- a/doc/src/fix_pair_tracker.rst +++ /dev/null @@ -1,124 +0,0 @@ -.. index:: fix pair/tracker - -fix pair/tracker command -======================== - -Syntax -"""""" - -.. parsed-literal:: - - fix ID group-ID pair/tracker N attribute1 attribute2 ... keyword values ... - -* ID, group-ID are documented in :doc:`fix ` command -* pair/tracker = style name of this fix command -* N = prepare data for output every this many timesteps -* one or more attributes may be appended - - .. parsed-literal:: - - possible attributes = id1 id2 time/created time/broken time/total - rmin rave x y z - - .. parsed-literal:: - - id1, id2 = IDs of the 2 atoms in each pair interaction - time/created = the time that the 2 atoms began interacting - time/broken = the time that the 2 atoms stopped interacting - time/total = the total time the 2 atoms interacted - r/min = the minimum radial distance between the 2 atoms during the interaction - r/ave = the average radial distance between the 2 atoms during the interaction - x, y, z = the center of mass position of the 2 atoms when they stopped interacting - -* zero or more keyword/value pairs may be appended -* keyword = *time/min* or *type/include* - - .. parsed-literal:: - - *time/min* value = T - T = minimum interaction time - *type/include* value = arg1 arg2 - arg = separate lists of types (see below) - -Examples -"""""""" - -.. code-block:: LAMMPS - - fix 1 all pair/tracker 1000 id1 id2 time/min 100 - fix 1 all pair/tracker 1000 time/created time/broken type/include 1 * type/include 2 3,4 - -Description -""""""""""" - -Tracks properties of pairwise interactions between two atoms and records data -whenever the atoms move beyond the interaction cutoff. -Must be used in conjunction with :doc:`pair tracker `. -Data is accumulated over a span of *N* timesteps before being deleted. -The number of datums generated, aggregated across all processors, equals -the number of broken interactions. Interactions are only included if both -atoms are included in the specified fix group. Additional filters can be -applied using the *time/min* or *type/include* keywords described below. - -.. note:: - - For extremely long-lived interactions, the calculation of *r/ave* may not be - correct due to double overflow. - -The *time/min* keyword sets a minimum amount of time that an interaction must -persist to be included. This setting can be used to censor short-lived interactions. -The *type/include* keyword filters interactions based on the types of the two atoms. -Data is only saved for interactions between atoms with types in the two lists. -Each list consists of a series of type -ranges separated by commas. The range can be specified as a -single numeric value, or a wildcard asterisk can be used to specify a range -of values. This takes the form "\*" or "\*n" or "n\*" or "m\*n". For -example, if M = the number of atom types, then an asterisk with no numeric -values means all types from 1 to M. A leading asterisk means all types -from 1 to n (inclusive). A trailing asterisk means all types from n to M -(inclusive). A middle asterisk means all types from m to n (inclusive). -Multiple *type/include* keywords may be added. - ----------- - -Restart, fix_modify, run start/stop, minimize info -""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" - -No information about this fix is written to :doc:`binary restart files `. -None of the :doc:`fix_modify ` options are -relevant to this fix. -No parameter of this fix can be used with the *start/stop* keywords of -the :doc:`run ` command. - -Output info -""""""""""" - -This compute calculates a local vector or local array depending on the -number of input values. The length of the vector or number of rows in -the array is the number of recorded, lost interactions. If a single input is -specified, a local vector is produced. If two or more inputs are -specified, a local array is produced where the number of columns = the -number of inputs. The vector or array can be accessed by any command -that uses local values from a compute as input. See the :doc:`Howto output ` page for an overview of LAMMPS output -options. - -The vector or array values will be doubles that correspond to the -specified attribute. - -Restrictions -"""""""""""" - -Must be used in conjunction with :doc:`pair style tracker `. - -This fix is part of the MISC package. It is only enabled if LAMMPS -was built with that package. See the :doc:`Build package ` page for more info. - -Related commands -"""""""""""""""" - -:doc:`pair tracker ` - -Default -""""""" - -none diff --git a/doc/src/fix_store_local.rst b/doc/src/fix_store_local.rst new file mode 100644 index 0000000000..6048222754 --- /dev/null +++ b/doc/src/fix_store_local.rst @@ -0,0 +1,77 @@ +.. index:: fix store/local + +fix store/local command +======================== + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group-ID store/local N nvalues + +* ID, group-ID are documented in :doc:`fix ` command +* store/local = style name of this fix command +* N = prepare data for output every this many timesteps +* nvalues = number of values stored by this fix + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all store/local 1000 2 + dump 1 all local 1000 dump.local f_1[1] f_1[2] + +Description +""""""""""" + +This fix provides the ability to store local data produced by +some LAMMPS commands including some pair and bond styles so it can be output. +Data is accumulated over a span of *N* timesteps before being deleted. +The number of datums generated, aggregated across all processors, depends on +the associated commands. Data is only included if it is generated from atoms +within the fix group-ID. + +---------- + +Restart, fix_modify, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +Accumulated local data is written to :doc:`binary restart files `. +None of the :doc:`fix_modify ` options are +relevant to this fix. +No parameter of this fix can be used with the *start/stop* keywords of +the :doc:`run ` command. + +Output info +""""""""""" + +This compute calculates a local vector or local array depending on the +number of input values. The length of the vector or number of rows in +the array is the number of recorded, lost interactions. If a single input is +specified, a local vector is produced. If two or more inputs are +specified, a local array is produced where the number of columns = the +number of inputs. The vector or array can be accessed by any command +that uses local values from a compute as input. See the :doc:`Howto output ` page for an overview of LAMMPS output +options. + +The vector or array values will be doubles that correspond to the +specified attribute. + +Restrictions +"""""""""""" + +Must be used in conjunction with another LAMMPS class which outputs local data. + +Related commands +"""""""""""""""" + +:doc:`pair tracker ` +:doc:`bond bpm/rotational ` +:doc:`dump local ` + +Default +""""""" + +none diff --git a/doc/src/fix_update_special_bonds.rst b/doc/src/fix_update_special_bonds.rst new file mode 100755 index 0000000000..f129b537d6 --- /dev/null +++ b/doc/src/fix_update_special_bonds.rst @@ -0,0 +1,55 @@ +.. index:: fix update/special/bonds + +fix update/special/bonds command +====================== + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group-ID update/special/bonds + +* ID, group-ID are documented in :doc:`fix ` command +* update/special/bonds = style name of this fix command + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all update/special/bonds + +Description +""""""""""" + +This fix is used to update the 1-2 special bond list for BPM bond styles. +This feature is used to censor pair forces between bonded particles. +See the :doc:`BPM how to ` for more information. + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +No information about this fix is written to :doc:`binary restart files `. None of the :doc:`fix_modify ` options +are relevant to this fix. No global or per-atom quantities are stored +by this fix for access by various :doc:`output commands `. +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 :doc:`newton ` bond off. + +Related commands +"""""""""""""""" + +:doc:`bond bpm/rotational ` + +Default +""""""" + +none + diff --git a/doc/src/pair_tracker.rst b/doc/src/pair_tracker.rst index d2cee1d879..c176c6d132 100644 --- a/doc/src/pair_tracker.rst +++ b/doc/src/pair_tracker.rst @@ -8,29 +8,51 @@ Syntax .. code-block:: LAMMPS - pair_style tracker keyword + pair_style tracker fix_ID keyword values attribute1 attribute2 ... -* zero or more keyword/arg pairs may be appended -* keyword = *finite* +* fix_ID = ID of associated fix store/local +* zero or more keywords may be appended +* keyword = *finite* or *time/min* or *type/include* .. parsed-literal:: *finite* value = none pair style uses atomic diameters to identify contacts + *time/min* value = T + T = minimum interaction time + *type/include* value = arg1 arg2 + arg = separate lists of types (see below) + +* one or more attributes may be appended + + .. parsed-literal:: + + possible attributes = id1 id2 time/created time/broken time/total + r/min r/ave x y z + + .. parsed-literal:: + + id1, id2 = IDs of the 2 atoms in each pair interaction + time/created = the time that the 2 atoms began interacting + time/broken = the time that the 2 atoms stopped interacting + time/total = the total time the 2 atoms interacted + r/min = the minimum radial distance between the 2 atoms during the interaction + r/ave = the average radial distance between the 2 atoms during the interaction + x, y, z = the center of mass position of the 2 atoms when they stopped interacting Examples """""""" .. code-block:: LAMMPS - pair_style hybrid/overlay tracker ... + pair_style hybrid/overlay tracker myfix id1 id2 type/include 1 * type/include 2 3,4 pair_coeff 1 1 tracker 2.0 - pair_style hybrid/overlay tracker finite ... + pair_style hybrid/overlay tracker myfix finite x y z time/min 100 pair_coeff * * tracker - fix 1 all pair/tracker 1000 time/created time/broken - dump 1 all local 1000 dump.local f_1[1] f_1[2] + fix myfix all store/local 1000 3 + dump 1 all local 1000 dump.local f_myfix[1] f_myfix[2] f_myfix[3] dump_modify 1 write_header no Description @@ -40,8 +62,11 @@ Style *tracker* monitors information about pairwise interactions. It does not calculate any forces on atoms. :doc:`Pair hybrid/overlay ` can be used to combine this pair style with another pair style. Style *tracker* must be used in conjunction -with about :doc:`fix pair_tracker ` which contains -information on what data can be output. +with :doc:`fix store/local ` which contains +information on outputting data. Whenever two neighboring atoms move beyond +the interaction cutoffPairwise data is processed and transferred to the associated +instance of :doc:`fix store/local `. Additional filters can +be applied using the *time/min* or *type/include* keywords described below. If the *finite* keyword is not defined, the following coefficients must be defined for each pair of atom types via the :doc:`pair_coeff ` @@ -55,6 +80,24 @@ If the *finite* keyword is defined, no coefficients may be defined. Interaction cutoffs are alternatively calculated based on the diameter of finite particles. +.. note:: + + For extremely long-lived interactions, the calculation of *r/ave* may not be + correct due to double overflow. + +The *time/min* keyword sets a minimum amount of time that an interaction must +persist to be included. This setting can be used to censor short-lived interactions. +The *type/include* keyword filters interactions based on the types of the two atoms. +Data is only saved for interactions between atoms with types in the two lists. +Each list consists of a series of type +ranges separated by commas. The range can be specified as a +single numeric value, or a wildcard asterisk can be used to specify a range +of values. This takes the form "\*" or "\*n" or "n\*" or "m\*n". For +example, if M = the number of atom types, then an asterisk with no numeric +values means all types from 1 to M. A leading asterisk means all types +from 1 to n (inclusive). A trailing asterisk means all types from n to M +(inclusive). A middle asterisk means all types from m to n (inclusive). +Multiple *type/include* keywords may be added. Mixing, shift, table, tail correction, restart, rRESPA info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/src/.gitignore b/src/.gitignore index 6c0a838c1b..1ed30868ab 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -251,6 +251,23 @@ /pair_mesont_tpm.cpp /pair_mesont_tpm.h +/atom_vec_sphere_bpm.cpp +/atom_vec_sphere_bpm.h +/bond_bpm.cpp +/bond_bpm.h +/bond_bpm_rotational.cpp +/bond_bpm_rotational.h +/bond_bpm_simple.cpp +/bond_bpm_simple.h +/compute_nbond_atom.cpp +/compute_nbond_atom.h +/fix_bond_history.cpp +/fix_bond_history.h +/fix_nve_sphere_bpm.cpp +/fix_nve_sphere_bpm.h +/fix_update_special_bonds.cpp +/fix_update_special_bonds.h + /compute_adf.cpp /compute_adf.h /compute_contact_atom.cpp @@ -779,8 +796,6 @@ /fix_orient_eco.h /fix_orient_fcc.cpp /fix_orient_fcc.h -/fix_pair_tracker.cpp -/fix_pair_tracker.h /fix_peri_neigh.cpp /fix_peri_neigh.h /fix_phonon.cpp diff --git a/src/BPM/atom_vec_sphere_bpm.cpp b/src/BPM/atom_vec_sphere_bpm.cpp new file mode 100644 index 0000000000..3fe372e340 --- /dev/null +++ b/src/BPM/atom_vec_sphere_bpm.cpp @@ -0,0 +1,254 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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. +------------------------------------------------------------------------- */ + +#include "atom_vec_sphere_bpm.h" + +#include "atom.h" +#include "error.h" +#include "fix.h" +#include "fix_adapt.h" +#include "math_const.h" +#include "modify.h" +#include "utils.h" + +#include + +using namespace LAMMPS_NS; +using namespace MathConst; + +/* ---------------------------------------------------------------------- */ + +AtomVecSphereBPM::AtomVecSphereBPM(LAMMPS *lmp) : AtomVec(lmp) +{ + mass_type = PER_ATOM; + molecular = Atom::MOLECULAR; + bonds_allow = 1; + + atom->molecule_flag = 1; + atom->sphere_flag = 1; + atom->radius_flag = atom->rmass_flag = atom->omega_flag = + atom->torque_flag = atom->quat_flag = 1; + + // strings with peratom variables to include in each AtomVec method + // strings cannot contain fields in corresponding AtomVec default strings + // order of fields in a string does not matter + // except: fields_data_atom & fields_data_vel must match data file + + fields_grow = (char *) + "molecule num_bond bond_type bond_atom nspecial special radius rmass omega torque quat"; + fields_copy = (char *) + "molecule num_bond bond_type bond_atom nspecial special radius rmass omega quat"; + fields_comm = (char *) ""; + fields_comm_vel = (char *) "omega quat"; + fields_reverse = (char *) "torque"; + fields_border = (char *) "molecule radius rmass"; + fields_border_vel = (char *) "molecule radius rmass omega quat"; + fields_exchange = (char *) + "molecule num_bond bond_type bond_atom nspecial special radius rmass omega quat"; + fields_restart = (char *) + "molecule num_bond bond_type bond_atom radius rmass omega quat"; + fields_create = (char *) + "molecule num_bond nspecial radius rmass omega quat"; + fields_data_atom = (char *) "id molecule type radius rmass x"; + fields_data_vel = (char *) "id v omega"; + + bond_per_atom = 0; + bond_negative = NULL; +} + +/* ---------------------------------------------------------------------- + process sub-style args + optional arg = 0/1 for static/dynamic particle radii +------------------------------------------------------------------------- */ + +void AtomVecSphereBPM::process_args(int narg, char **arg) +{ + if (narg != 0 && narg != 1) + error->all(FLERR,"Illegal atom_style sphere/bpm command"); + + radvary = 0; + if (narg == 1) { + radvary = utils::numeric(FLERR,arg[0],true,lmp); + if (radvary < 0 || radvary > 1) + error->all(FLERR,"Illegal atom_style sphere/bpm command"); + } + + // dynamic particle radius and mass must be communicated every step + + if (radvary) { + fields_comm = (char *) "radius rmass"; + fields_comm_vel = (char *) "radius rmass omega"; + } + + // delay setting up of fields until now + + setup_fields(); +} + +/* ---------------------------------------------------------------------- */ + +void AtomVecSphereBPM::init() +{ + AtomVec::init(); + + // check if optional radvary setting should have been set to 1 + + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"adapt") == 0) { + FixAdapt *fix = (FixAdapt *) modify->fix[i]; + if (fix->diamflag && radvary == 0) + error->all(FLERR,"Fix adapt changes particle radii " + "but atom_style sphere is not dynamic"); + } +} + +/* ---------------------------------------------------------------------- + set local copies of all grow ptrs used by this class, except defaults + needed in replicate when 2 atom classes exist and it calls pack_restart() +------------------------------------------------------------------------- */ + +void AtomVecSphereBPM::grow_pointers() +{ + radius = atom->radius; + rmass = atom->rmass; + omega = atom->omega; + quat = atom->quat; + + num_bond = atom->num_bond; + bond_type = atom->bond_type; + nspecial = atom->nspecial; +} + + +/* ---------------------------------------------------------------------- + initialize non-zero atom quantities +------------------------------------------------------------------------- */ + +void AtomVecSphereBPM::create_atom_post(int ilocal) +{ + radius[ilocal] = 0.5; + rmass[ilocal] = 4.0*MY_PI/3.0 * 0.5*0.5*0.5; + + quat[ilocal][0] = 1.0; + quat[ilocal][1] = 0.0; + quat[ilocal][2] = 0.0; + quat[ilocal][3] = 0.0; + +} + +/* ---------------------------------------------------------------------- + modify values for AtomVec::pack_restart() to pack +------------------------------------------------------------------------- */ + +void AtomVecSphereBPM::pack_restart_pre(int ilocal) +{ + // insure bond_negative vector is needed length + + if (bond_per_atom < atom->bond_per_atom) { + delete [] bond_negative; + bond_per_atom = atom->bond_per_atom; + bond_negative = new int[bond_per_atom]; + } + + // flip any negative types to positive and flag which ones + + any_bond_negative = 0; + for (int m = 0; m < num_bond[ilocal]; m++) { + if (bond_type[ilocal][m] < 0) { + bond_negative[m] = 1; + bond_type[ilocal][m] = -bond_type[ilocal][m]; + any_bond_negative = 1; + } else bond_negative[m] = 0; + } +} + +/* ---------------------------------------------------------------------- + unmodify values packed by AtomVec::pack_restart() +------------------------------------------------------------------------- */ + +void AtomVecSphereBPM::pack_restart_post(int ilocal) +{ + // restore the flagged types to their negative values + + if (any_bond_negative) { + for (int m = 0; m < num_bond[ilocal]; m++) + if (bond_negative[m]) bond_type[ilocal][m] = -bond_type[ilocal][m]; + } +} + +/* ---------------------------------------------------------------------- + initialize other atom quantities after AtomVec::unpack_restart() +------------------------------------------------------------------------- */ + +void AtomVecSphereBPM::unpack_restart_init(int ilocal) +{ + nspecial[ilocal][0] = 0; + nspecial[ilocal][1] = 0; + nspecial[ilocal][2] = 0; +} + +/* ---------------------------------------------------------------------- + modify what AtomVec::data_atom() just unpacked + or initialize other atom quantities +------------------------------------------------------------------------- */ + +void AtomVecSphereBPM::data_atom_post(int ilocal) +{ + radius_one = 0.5 * atom->radius[ilocal]; + radius[ilocal] = radius_one; + if (radius_one > 0.0) + rmass[ilocal] *= 4.0*MY_PI/3.0 * radius_one*radius_one*radius_one; + + if (rmass[ilocal] <= 0.0) + error->one(FLERR,"Invalid density in Atoms section of data file"); + + omega[ilocal][0] = 0.0; + omega[ilocal][1] = 0.0; + omega[ilocal][2] = 0.0; + + quat[ilocal][0] = 1.0; + quat[ilocal][1] = 0.0; + quat[ilocal][2] = 0.0; + quat[ilocal][3] = 0.0; + + num_bond[ilocal] = 0; + nspecial[ilocal][0] = 0; + nspecial[ilocal][1] = 0; + nspecial[ilocal][2] = 0; +} + +/* ---------------------------------------------------------------------- + modify values for AtomVec::pack_data() to pack +------------------------------------------------------------------------- */ + +void AtomVecSphereBPM::pack_data_pre(int ilocal) +{ + radius_one = radius[ilocal]; + rmass_one = rmass[ilocal]; + + radius[ilocal] *= 2.0; + if (radius_one!= 0.0) + rmass[ilocal] = + rmass_one / (4.0*MY_PI/3.0 * radius_one*radius_one*radius_one); +} + +/* ---------------------------------------------------------------------- + unmodify values packed by AtomVec::pack_data() +------------------------------------------------------------------------- */ + +void AtomVecSphereBPM::pack_data_post(int ilocal) +{ + radius[ilocal] = radius_one; + rmass[ilocal] = rmass_one; +} diff --git a/src/BPM/atom_vec_sphere_bpm.h b/src/BPM/atom_vec_sphere_bpm.h new file mode 100644 index 0000000000..3a84e3999e --- /dev/null +++ b/src/BPM/atom_vec_sphere_bpm.h @@ -0,0 +1,83 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 ATOM_CLASS +// clang-format off +AtomStyle(sphere/bpm,AtomVecSphereBPM) +// clang-format on +#else + +#ifndef LMP_ATOM_VEC_SPHERE_BPM_H +#define LMP_ATOM_VEC_SPHERE_BPM_H + +#include "atom_vec.h" + +namespace LAMMPS_NS { + +class AtomVecSphereBPM : public AtomVec { + public: + AtomVecSphereBPM(class LAMMPS *); + void process_args(int, char **); + void init(); + + void grow_pointers(); + void create_atom_post(int); + void pack_restart_pre(int); + void pack_restart_post(int); + void unpack_restart_init(int); + void data_atom_post(int); + void pack_data_pre(int); + void pack_data_post(int); + + + private: + int *num_bond; + int **bond_type; + int **nspecial; + + double *radius,*rmass; + double **omega, **torque, **quat; + + int any_bond_negative; + int bond_per_atom; + int *bond_negative; + + int radvary; + double radius_one,rmass_one; +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Per-processor system is too big + +The number of owned atoms plus ghost atoms on a single +processor must fit in 32-bit integer. + +E: Invalid atom type in Atoms section of data file + +Atom types must range from 1 to specified # of types. + +E: Invalid radius in Atoms section of data file + +Radius must be >= 0.0. + +E: Invalid density in Atoms section of data file + +Density value cannot be <= 0.0. + +*/ diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp new file mode 100644 index 0000000000..7a0d4fedd6 --- /dev/null +++ b/src/BPM/bond_bpm.cpp @@ -0,0 +1,349 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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. +------------------------------------------------------------------------- */ + +#include "bond_bpm.h" + +#include "atom.h" +#include "domain.h" +#include "error.h" +#include "fix_store_local.h" +#include "fix_update_special_bonds.h" +#include "force.h" +#include "memory.h" +#include "modify.h" +#include "update.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp) +{ + id_fix_store_local = nullptr; + id_fix_prop_atom = nullptr; + fix_store_local = nullptr; + fix_update_special_bonds = nullptr; + + prop_atom_flag = 0; + nvalues = 0; + output_data = nullptr; + pack_choice = nullptr; + + r0_max_estimate = 0.0; + max_stretch = 1.0; +} + +/* ---------------------------------------------------------------------- */ + +BondBPM::~BondBPM() +{ + delete [] pack_choice; + delete [] id_fix_store_local; + delete [] id_fix_prop_atom; + memory->destroy(output_data); +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::init_style() +{ + int ifix; + if (id_fix_store_local) { + ifix = modify->find_fix(id_fix_store_local); + if (ifix < 0) error->all(FLERR, "Cannot find fix store/local"); + if (strcmp(modify->fix[ifix]->style, "store/local") != 0) + error->all(FLERR, "Incorrect fix style matched, not store/local"); + fix_store_local = (FixStoreLocal *) modify->fix[ifix]; + fix_store_local->nvalues = nvalues; + } + + ifix = modify->find_fix_by_style("update/special/bonds"); + if (ifix >= 0) + fix_update_special_bonds = (FixUpdateSpecialBonds *) modify->fix[ifix]; + else + fix_update_special_bonds = nullptr; + + + if (force->angle || force->dihedral || force->improper) + error->all(FLERR, + "Bond style bpm cannot be used with 3,4-body interactions"); + if (atom->molecular == 2) + error->all(FLERR, + "Bond style bpm cannot be used with atom style template"); + + // special 1-3 and 1-4 weights must be 1 to prevent building 1-3 and 1-4 special bond lists + if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || + force->special_coul[2] != 1.0 || force->special_coul[3] != 1.0) + error->all(FLERR,"Bond style bpm requires 1-3 and 1-4 special weights of 1.0"); +} + +/* ---------------------------------------------------------------------- + global settings + All args before store/local command are saved for potential args + for specific bond BPM substyles + All args after optional store/local command are variables stored + in the compute store/local +------------------------------------------------------------------------- */ + +void BondBPM::settings(int narg, char **arg) +{ + int iarg = 0; + while (iarg < narg) { + if (id_fix_store_local) { + if (strcmp(arg[iarg], "id1") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_id1; + } else if (strcmp(arg[iarg], "id2") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_id2; + } else if (strcmp(arg[iarg], "time") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_time; + } else if (strcmp(arg[iarg], "x") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_x; + } else if (strcmp(arg[iarg], "y") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_y; + } else if (strcmp(arg[iarg], "z") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_z; + } else if (strcmp(arg[iarg], "x/ref") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_x_ref; + prop_atom_flag = 1; + } else if (strcmp(arg[iarg], "y/ref") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_y_ref; + prop_atom_flag = 1; + } else if (strcmp(arg[iarg], "z/ref") == 0) { + pack_choice[nvalues++] = &BondBPM::pack_z_ref; + prop_atom_flag = 1; + } else { + error->all(FLERR, "Illegal bond_style command"); + } + } else if (strcmp(arg[iarg], "store/local") == 0) { + id_fix_store_local = utils::strdup(arg[iarg+1]); + iarg ++; + nvalues = 0; + pack_choice = new FnPtrPack[narg - iarg - 1]; + } + iarg ++; + } + + if (id_fix_store_local) { + if (nvalues == 0) error->all(FLERR, + "Bond style bpm/rotational must include at least one value to output"); + memory->create(output_data, nvalues, "bond/bpm:output_data"); + + // Use store property to save reference positions as it can transfer to ghost atoms + if (prop_atom_flag == 1) { + + id_fix_prop_atom = utils::strdup("BPM_PROPERTY_ATOM"); + int ifix = modify->find_fix(id_fix_prop_atom); + if (ifix < 0) { + modify->add_fix(fmt::format("{} all property/atom " + "d_BPM_X_REF d_BPM_Y_REF d_BPM_Z_REF ghost yes", id_fix_prop_atom)); + ifix = modify->find_fix(id_fix_prop_atom); + } + + int type_flag; + int col_flag; + index_x_ref = atom->find_custom("BPM_X_REF", type_flag, col_flag); + index_y_ref = atom->find_custom("BPM_Y_REF", type_flag, col_flag); + index_z_ref = atom->find_custom("BPM_Z_REF", type_flag, col_flag); + + if (modify->fix[ifix]->restart_reset) { + modify->fix[ifix]->restart_reset = 0; + } else { + double *x_ref = atom->dvector[index_x_ref]; + double *y_ref = atom->dvector[index_y_ref]; + double *z_ref = atom->dvector[index_z_ref]; + + double **x = atom->x; + for (int i = 0; i < atom->nlocal; i++) { + x_ref[i] = x[i][0]; + y_ref[i] = x[i][1]; + z_ref[i] = x[i][2]; + } + } + } + } +} + +/* ---------------------------------------------------------------------- + used to check bond communiction cutoff - not perfect, estimates based on local-local only +------------------------------------------------------------------------- */ + +double BondBPM::equilibrium_distance(int i) +{ + // Ghost atoms may not yet be communicated, this may only be an estimate + if (r0_max_estimate == 0) { + int type, j; + double delx, dely, delz, r; + double **x = atom->x; + for (int i = 0; i < atom->nlocal; i ++) { + for (int m = 0; m < atom->num_bond[i]; m ++) { + type = atom->bond_type[i][m]; + if (type == 0) continue; + + j = atom->map(atom->bond_atom[i][m]); + if(j == -1) continue; + + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + domain->minimum_image(delx, dely, delz); + + r = sqrt(delx*delx + dely*dely + delz*delz); + if(r > r0_max_estimate) r0_max_estimate = r; + } + } + + double temp; + MPI_Allreduce(&r0_max_estimate,&temp,1,MPI_DOUBLE,MPI_MAX,world); + r0_max_estimate = temp; + + //if (comm->me == 0) + // utils::logmesg(lmp,fmt::format("Estimating longest bond = {}\n",r0_max_estimate)); + } + + // Divide out heuristic prefactor added in comm class + return max_stretch*r0_max_estimate/1.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::process_broken(int i, int j) +{ + if (fix_store_local) { + for (int n = 0; n < nvalues; n++) + (this->*pack_choice[n])(n, i, j); + + fix_store_local->add_data(output_data, i, j); + } + + if (fix_update_special_bonds) + fix_update_special_bonds->add_broken_bond(i, j); + + // Manually search and remove from atom arrays + // need to remove in case special bonds arrays rebuilt + int m, n; + int nlocal = atom->nlocal; + + tagint *tag = atom->tag; + tagint **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + int *num_bond = atom->num_bond; + + if (i < nlocal) { + for (m = 0; m < num_bond[i]; m++) { + // if(tag[i] == 25 or tag[j] == 25) printf("\t 1 ATOM LOOP %d-%d, %d/%d \n", tag[i], bond_atom[i][m], m, num_bond[i]); + if (bond_atom[i][m] == tag[j]) { + // if(tag[i] == 25 or tag[j] == 25) printf("\t DELETED\n"); + bond_type[i][m] = 0; + n = num_bond[i]; + bond_type[i][m] = bond_type[i][n-1]; + bond_atom[i][m] = bond_atom[i][n-1]; + num_bond[i]--; + break; + } + } + } + + if (j < nlocal) { + for (m = 0; m < num_bond[j]; m++) { + // if(tag[i] == 25 or tag[j] == 25) printf("\t 2 ATOM LOOP %d-%d, %d/%d \n", tag[j], bond_atom[j][m], m, num_bond[j]); + if (bond_atom[j][m] == tag[i]) { + // if(tag[j] == 25 or tag[j] == 25) printf("\t DELETED\n"); + bond_type[j][m] = 0; + n = num_bond[j]; + bond_type[j][m] = bond_type[j][n-1]; + bond_atom[j][m] = bond_atom[j][n-1]; + num_bond[j]--; + break; + } + } + } + + // if((tag[i] == 10 and tag[j] == 23)or (tag[i] == 23 and tag[j] == 10)) printf("Broken bond %d-%d (%d %d/%d)\n", tag[i], tag[j], i, j, nlocal); + +} + +/* ---------------------------------------------------------------------- + one method for every keyword bond bpm can output + the atom property is packed into array or vector +------------------------------------------------------------------------- */ + +void BondBPM::pack_id1(int n, int i, int j) +{ + tagint *tag = atom->tag; + output_data[n] = tag[i]; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_id2(int n, int i, int j) +{ + tagint *tag = atom->tag; + output_data[n] = tag[j]; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_time(int n, int i, int j) +{ + bigint time = update->ntimestep; + output_data[n] = time; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_x(int n, int i, int j) +{ + double **x = atom->x; + output_data[n] = (x[i][0] + x[j][0])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_y(int n, int i, int j) +{ + double **x = atom->x; + output_data[n] = (x[i][1] + x[j][1])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_z(int n, int i, int j) +{ + double **x = atom->x; + output_data[n] = (x[i][2] + x[j][2])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_x_ref(int n, int i, int j) +{ + double *x = atom->dvector[index_x_ref]; + output_data[n] = (x[i] + x[j])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_y_ref(int n, int i, int j) +{ + double *y = atom->dvector[index_y_ref]; + output_data[n] = (y[i] + y[j])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void BondBPM::pack_z_ref(int n, int i, int j) +{ + double *z = atom->dvector[index_z_ref]; + output_data[n] = (z[i] + z[j])*0.5; +} diff --git a/src/BPM/bond_bpm.h b/src/BPM/bond_bpm.h new file mode 100644 index 0000000000..d77e44e1d6 --- /dev/null +++ b/src/BPM/bond_bpm.h @@ -0,0 +1,90 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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_BOND_BPM_H +#define LMP_BOND_BPM_H + +#include "bond.h" + +namespace LAMMPS_NS { + +class BondBPM : public Bond { + public: + BondBPM(class LAMMPS *); + virtual ~BondBPM(); + virtual void compute(int, int) = 0; + virtual void coeff(int, char **) = 0; + virtual void init_style(); + void settings(int, char **); + double equilibrium_distance(int); + void write_restart(FILE *){}; + void read_restart(FILE *){}; + void write_data(FILE *) {}; + double single(int, double, int, int, double &) = 0; + + protected: + double r0_max_estimate; + double max_stretch; + + char *id_fix_store_local, *id_fix_prop_atom; + class FixStoreLocal *fix_store_local; + class FixUpdateSpecialBonds *fix_update_special_bonds; + + void process_broken(int, int); + typedef void (BondBPM::*FnPtrPack)(int,int,int); + FnPtrPack *pack_choice; // ptrs to pack functions + double *output_data; + + int prop_atom_flag, nvalues; + int index_x_ref, index_y_ref, index_z_ref; + + void pack_id1(int,int,int); + void pack_id2(int,int,int); + void pack_time(int,int,int); + void pack_x(int,int,int); + void pack_y(int,int,int); + void pack_z(int,int,int); + void pack_x_ref(int,int,int); + void pack_y_ref(int,int,int); + void pack_z_ref(int,int,int); +}; + +} // namespace LAMMPS_NS + +#endif + +/* ERROR/WARNING messages: + +E: Cannot find fix store/local + +Fix id cannot be found. + +E: Illegal bond_style command + +Self-explanatory. + +E: Bond style bpm/rotational must include at least one value to output + +Must include at least one bond property to store in fix store/local + +E: Bond style bpm/rotational cannot be used with 3,4-body interactions + +No angle, dihedral, or improper styles can be defined when using +bond style bpm/rotational. + +E: Bond style bpm/rotational cannot be used with atom style template + +This bond style can change the bond topology which is not +allowed with this atom style. + +*/ diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp new file mode 100644 index 0000000000..bee7a69b6a --- /dev/null +++ b/src/BPM/bond_bpm_rotational.cpp @@ -0,0 +1,727 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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. +------------------------------------------------------------------------- */ + +#include "bond_bpm_rotational.h" + +#include "atom.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "fix_bond_history.h" +#include "force.h" +#include "math_const.h" +#include "math_extra.h" +#include "memory.h" +#include "modify.h" +#include "neighbor.h" +#include "pair.h" + +#include +#include +#include + +#define EPSILON 1e-10 + +using namespace LAMMPS_NS; +using namespace MathExtra; + + + + + +#include "update.h" +/* ---------------------------------------------------------------------- */ + +BondBPMRotational::BondBPMRotational(LAMMPS *lmp) : BondBPM(lmp) +{ + partial_flag = 1; + fix_bond_history = nullptr; +} + +/* ---------------------------------------------------------------------- */ + +BondBPMRotational::~BondBPMRotational() +{ + if(fix_bond_history) modify->delete_fix("BOND_HISTORY_BPM_ROTATIONAL"); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(Kr); + memory->destroy(Ks); + memory->destroy(Kt); + memory->destroy(Kb); + memory->destroy(Fcr); + memory->destroy(Fcs); + memory->destroy(Tct); + memory->destroy(Tcb); + memory->destroy(gnorm); + memory->destroy(gslide); + memory->destroy(groll); + memory->destroy(gtwist); + } +} + +/* ---------------------------------------------------------------------- */ + +double BondBPMRotational::acos_limit(double c) +{ + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + return acos(c); +} + +/* ---------------------------------------------------------------------- + Store data for a single bond - if bond added after LAMMPS init (e.g. pour) +------------------------------------------------------------------------- */ + +double BondBPMRotational::store_bond(int n,int i,int j) +{ + int m,k; + double delx, dely, delz, r, rinv; + double **x = atom->x; + tagint *tag = atom->tag; + double **bondstore = fix_bond_history->bondstore; + + if (tag[i] < tag[j]) { + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + } else { + delx = x[j][0] - x[i][0]; + dely = x[j][1] - x[i][1]; + delz = x[j][2] - x[i][2]; + } + + r = sqrt(delx*delx + dely*dely + delz*delz); + rinv = 1.0/r; + + bondstore[n][0] = r; + bondstore[n][1] = delx*rinv; + bondstore[n][2] = dely*rinv; + bondstore[n][3] = delz*rinv; + + if (i < atom->nlocal) { + for (m = 0; m < atom->num_bond[i]; m ++) { + if (atom->bond_atom[i][m] == tag[j]) { + fix_bond_history->update_atom_value(i, m, 0, r); + fix_bond_history->update_atom_value(i, m, 1, delx*rinv); + fix_bond_history->update_atom_value(i, m, 2, dely*rinv); + fix_bond_history->update_atom_value(i, m, 3, delz*rinv); + } + } + } + + if (j < atom->nlocal) { + for (m = 0; m < atom->num_bond[j]; m ++) { + if (atom->bond_atom[j][m] == tag[i]) { + fix_bond_history->update_atom_value(j, m, 0, r); + fix_bond_history->update_atom_value(j, m, 1, delx*rinv); + fix_bond_history->update_atom_value(j, m, 2, dely*rinv); + fix_bond_history->update_atom_value(j, m, 3, delz*rinv); + } + } + } + + return r; +} + +/* ---------------------------------------------------------------------- + Store data for all bonds called once +------------------------------------------------------------------------- */ + +void BondBPMRotational::store_data() +{ + int i, j, m, type; + double delx, dely, delz, r, rinv; + double **x = atom->x; + int **bond_type = atom->bond_type; + tagint *tag = atom->tag; + + for (i = 0; i < atom->nlocal; i ++) { + for (m = 0; m < atom->num_bond[i]; m ++) { + type = bond_type[i][m]; + + //Skip if bond was turned off + if(type < 0) + continue; + + // map to find index n for tag + j = atom->map(atom->bond_atom[i][m]); + if(j == -1) error->one(FLERR, "Atom missing in BPM bond"); + + // Save orientation as pointing towards small tag + if(tag[i] < tag[j]){ + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + } else { + delx = x[j][0] - x[i][0]; + dely = x[j][1] - x[i][1]; + delz = x[j][2] - x[i][2]; + } + + // Get closest image in case bonded with ghost + domain->minimum_image(delx, dely, delz); + + r = sqrt(delx*delx + dely*dely + delz*delz); + rinv = 1.0/r; + fix_bond_history->update_atom_value(i, m, 0, r); + fix_bond_history->update_atom_value(i, m, 1, delx*rinv); + fix_bond_history->update_atom_value(i, m, 2, dely*rinv); + fix_bond_history->update_atom_value(i, m, 3, delz*rinv); + } + } + + fix_bond_history->post_neighbor(); +} + +/* ---------------------------------------------------------------------- */ + +void BondBPMRotational::compute(int eflag, int vflag) +{ + + if (! fix_bond_history->stored_flag) { + fix_bond_history->stored_flag = true; + store_data(); + } + + int i1,i2,itmp,m,n,type,itype,jtype; + double evdwl,fpair,rsq,ebond; + double q1[4], q2[4], r[3], r0[3]; + double r0_mag, r_mag, r_mag_inv, Fr_mag, Fs_mag; + double Tt_mag, Tb_mag; + double force1on2[3], torque1on2[3], torque2on1[3]; + double breaking, smooth, smooth_sq; + double rhat[3], wn1[3], wn2[3], wxn1[3], wxn2[3], vroll[3]; + double w1dotr, w2dotr, v1dotr, v2dotr; + double vn1[3], vn2[3], vt1[3], vt2[3], tmp[3], s1[3], s2[3], tdamp[3]; + double tor1, tor2, tor3, fs1, fs2, fs3; + + double q2inv[4], rb[3], rb_x_r0[3], s[3], t[3], Fs[3]; + double q21[4], qp21[4], Tbp[3], Ttp[3]; + double Tsp[3], Fsp[3], Tt[3], Tb[3], Ts[3], F_rot[3], T_rot[3]; + double mq[4], mqinv[4], Ttmp[3], Ftmp[3], qtmp[4]; + double r0_dot_rb, gamma, c, psi, theta, sin_phi, cos_phi, temp, mag_in_plane, mag_out_plane; + + ev_init(eflag,vflag); + + if (vflag_global == 2) + force->pair->vflag_either = force->pair->vflag_global = 1; + + double **cutsq = force->pair->cutsq; + double **x = atom->x; + double **v = atom->v; + double **omega = atom->omega; + double **f = atom->f; + double **torque = atom->torque; + double *radius = atom->radius; + double **quat = atom->quat; + tagint *tag = atom->tag; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + + double **bondstore = fix_bond_history->bondstore; + + for (n = 0; n < nbondlist; n++) { + + // skip bond if already broken + if (bondlist[n][2] <= 0) continue; + + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + type = bondlist[n][2]; + r0_mag = bondstore[n][0]; + + // Ensure pair is always ordered such that r0 points in + // a consistent direction and to ensure numerical operations + // are identical to minimize the possibility that a bond straddling + // an mpi grid (newton off) doesn't break on one proc but not the other + if (tag[i2] < tag[i1]) { + itmp = i1; + i1 = i2; + i2 = itmp; + } + + // If bond hasn't been set - should be initialized to zero + if (r0_mag < EPSILON || std::isnan(r0_mag)) + r0_mag = store_bond(n,i1,i2); + + r0[0] = bondstore[n][1]; + r0[1] = bondstore[n][2]; + r0[2] = bondstore[n][3]; + MathExtra::scale3(r0_mag, r0); + + q1[0] = quat[i1][0]; + q1[1] = quat[i1][1]; + q1[2] = quat[i1][2]; + q1[3] = quat[i1][3]; + + q2[0] = quat[i2][0]; + q2[1] = quat[i2][1]; + q2[2] = quat[i2][2]; + q2[3] = quat[i2][3]; + + // Note this is the reverse of Mora & Wang + MathExtra::sub3(x[i1], x[i2], r); + + rsq = MathExtra::lensq3(r); + r_mag = sqrt(rsq); + r_mag_inv = 1.0/r_mag; + MathExtra::scale3(r_mag_inv, r, rhat); + + // ------------------------------------------------------// + // Calculate forces using formulation in: + // 1) Y. Wang Acta Geotechnica 2009 + // 2) P. Mora & Y. Wang Advances in Geomcomputing 2009 + // ------------------------------------------------------// + + // Calculate normal forces, rb = bond vector in particle 1's frame + MathExtra::qconjugate(q2, q2inv); + MathExtra::quatrotvec(q2inv, r, rb); + Fr_mag = Kr[type]*(r_mag - r0_mag); + + MathExtra::scale3(Fr_mag*r_mag_inv, rb, F_rot); + + // Calculate forces due to tangential displacements (no rotation) + r0_dot_rb = dot3(r0, rb); + c = r0_dot_rb*r_mag_inv/r0_mag; + gamma = acos_limit(c); + + MathExtra::cross3(rb, r0, rb_x_r0); + MathExtra::cross3(rb, rb_x_r0, s); + MathExtra::norm3(s); + + MathExtra::scale3(Ks[type]*r_mag*gamma, s, Fs); + + // Calculate torque due to tangential displacements + MathExtra::cross3(r0, rb, t); + MathExtra::norm3(t); + + MathExtra::scale3(0.5*r_mag*Ks[type]*r_mag*gamma, t, Ts); + + // Relative rotation force/torque + // Use representation of X'Y'Z' rotations from Wang, Mora 2009 + temp = r_mag + rb[2]; + if (temp < 0.0) temp = 0.0; + mq[0] = sqrt(2)*0.5*sqrt(temp*r_mag_inv); + + temp = sqrt(rb[0]*rb[0]+rb[1]*rb[1]); + if (temp != 0.0) { + mq[1] = -sqrt(2)*0.5/temp; + temp = r_mag - rb[2]; + if (temp < 0.0) temp = 0.0; + mq[1] *= sqrt(temp*r_mag_inv); + mq[2] = -mq[1]; + mq[1] *= rb[1]; + mq[2] *= rb[0]; + } else { + // If aligned along z axis, x,y terms zero (r_mag-rb[2] = 0) + mq[1] = 0.0; + mq[2] = 0.0; + } + mq[3] = 0.0; + + // qp21 = opposite of r^\circ_21 in Wang + // q21 = opposite of r_21 in Wang + MathExtra::quatquat(q2inv, q1, qp21); + MathExtra::qconjugate(mq, mqinv); + MathExtra::quatquat(mqinv,qp21,qtmp); + MathExtra::quatquat(qtmp,mq,q21); + + temp = sqrt(q21[0]*q21[0] + q21[3]*q21[3]); + if (temp != 0.0) { + c = q21[0]/temp; + psi = 2.0*acos_limit(c); + } else { + c = 0.0; + psi = 0.0; + } + + // Map negative rotations + if (q21[3] < 0.0) // sin = q21[3]/temp + psi = -psi; + + if (q21[3] == 0.0) + psi = 0.0; + + c = q21[0]*q21[0] - q21[1]*q21[1] - q21[2]*q21[2] + q21[3]*q21[3]; + theta = acos_limit(c); + + // Separately calculte magnitude of quaternion in x-y and out of x-y planes + // to avoid dividing by zero + mag_out_plane = (q21[0]*q21[0] + q21[3]*q21[3]); + mag_in_plane = (q21[1]*q21[1] + q21[2]*q21[2]); + + if (mag_in_plane == 0.0) { + // No rotation => no bending/shear torque or extra shear force + // achieve by setting cos/sin = 0 + cos_phi = 0.0; + sin_phi = 0.0; + } else if (mag_out_plane == 0.0) { + // Calculate angle in plane + cos_phi = q21[2]/sqrt(mag_in_plane); + sin_phi = -q21[1]/sqrt(mag_in_plane); + } else { + // Default equations in Mora, Wang 2009 + cos_phi = q21[1]*q21[3] + q21[0]*q21[2]; + sin_phi = q21[2]*q21[3] - q21[0]*q21[1]; + + cos_phi /= sqrt(mag_out_plane*mag_in_plane); + sin_phi /= sqrt(mag_out_plane*mag_in_plane); + } + + Tbp[0] = -Kb[type]*theta*sin_phi; + Tbp[1] = Kb[type]*theta*cos_phi; + Tbp[2] = 0.0; + + Ttp[0] = 0.0; + Ttp[1] = 0.0; + Ttp[2] = Kt[type]*psi; + + Fsp[0] = -0.5*Ks[type]*r_mag*theta*cos_phi; + Fsp[1] = -0.5*Ks[type]*r_mag*theta*sin_phi; + Fsp[2] = 0.0; + + Tsp[0] = 0.25*Ks[type]*r_mag*r_mag*theta*sin_phi; + Tsp[1] = -0.25*Ks[type]*r_mag*r_mag*theta*cos_phi; + Tsp[2] = 0.0; + + // Rotate forces/torques back to 1st particle's frame + + MathExtra::quatrotvec(mq, Fsp, Ftmp); + MathExtra::quatrotvec(mq, Tsp, Ttmp); + for (m = 0; m < 3; m++) { + Fs[m] += Ftmp[m]; + Ts[m] += Ttmp[m]; + } + + MathExtra::quatrotvec(mq, Tbp, Tb); + MathExtra::quatrotvec(mq, Ttp, Tt); + + // Sum forces and calculate magnitudes + F_rot[0] += Fs[0]; + F_rot[1] += Fs[1]; + F_rot[2] += Fs[2]; + MathExtra::quatrotvec(q2, F_rot, force1on2); + + T_rot[0] = Ts[0] + Tt[0] + Tb[0]; + T_rot[1] = Ts[1] + Tt[1] + Tb[1]; + T_rot[2] = Ts[2] + Tt[2] + Tb[2]; + MathExtra::quatrotvec(q2, T_rot, torque1on2); + + T_rot[0] = Ts[0] - Tt[0] - Tb[0]; + T_rot[1] = Ts[1] - Tt[1] - Tb[1]; + T_rot[2] = Ts[2] - Tt[2] - Tb[2]; + MathExtra::quatrotvec(q2, T_rot, torque2on1); + + Fs_mag = MathExtra::len3(Fs); + Tt_mag = MathExtra::len3(Tt); + Tb_mag = MathExtra::len3(Tb); + + // ------------------------------------------------------// + // Check if bond breaks + // ------------------------------------------------------// + + if (r_mag < r0_mag) + breaking = Fs_mag/Fcs[type] + Tb_mag/Tcb[type] + Tt_mag/Tct[type]; + else + breaking = Fr_mag/Fcr[type] + Fs_mag/Fcs[type] + Tb_mag/Tcb[type] + Tt_mag/Tct[type]; + + if (breaking >= 1.0) { + bondlist[n][2] = 0; + process_broken(i1, i2); + continue; + } + + smooth = 1.0 - breaking; + smooth_sq = smooth*smooth; + + // Not actual energy, just a handy metric + if (eflag) ebond = -smooth_sq; + + // ------------------------------------------------------// + // Calculate damping using formulation in + // Y. Wang, F. Alonso-Marroquin, W. Guo 2015 + // ------------------------------------------------------// + // Note: n points towards 1 vs pointing towards 2 + + // Damp normal velocity difference + v1dotr = MathExtra::dot3(v[i1],rhat); + v2dotr = MathExtra::dot3(v[i2],rhat); + + MathExtra::scale3(v1dotr, rhat, vn1); + MathExtra::scale3(v2dotr, rhat, vn2); + + MathExtra::sub3(vn1, vn2, tmp); + MathExtra::scale3(gnorm[type], tmp); + MathExtra::add3(force1on2, tmp, force1on2); + + // Damp tangential objective velocities + MathExtra::sub3(v[i1], vn1, vt1); + MathExtra::sub3(v[i2], vn2, vt2); + + MathExtra::sub3(vt2, vt1, tmp); + MathExtra::scale3(-0.5, tmp); + + MathExtra::cross3(omega[i1], r, s1); + MathExtra::scale3(0.5, s1); + MathExtra::sub3(s1, tmp, s1); // Eq 12 + + MathExtra::cross3(omega[i2], r, s2); + MathExtra::scale3(-0.5,s2); + MathExtra::add3(s2, tmp, s2); // Eq 13 + MathExtra::scale3(-0.5,s2); + + MathExtra::sub3(s1, s2, tmp); + MathExtra::scale3(gslide[type], tmp); + MathExtra::add3(force1on2, tmp, force1on2); + + // Apply corresponding torque + MathExtra::cross3(r,tmp,tdamp); + MathExtra::scale3(-0.5, tdamp); // 0.5*r points from particle 2 to midpoint + MathExtra::add3(torque1on2, tdamp, torque1on2); + MathExtra::add3(torque2on1, tdamp, torque2on1); + + // Damp rolling + MathExtra::cross3(omega[i1], rhat, wxn1); + MathExtra::cross3(omega[i2], rhat, wxn2); + MathExtra::sub3(wxn1, wxn2, vroll); // Eq. 31 + MathExtra::cross3(r, vroll, tdamp); + + MathExtra::scale3(0.5*groll[type], tdamp); + MathExtra::add3(torque1on2, tdamp, torque1on2); + MathExtra::scale3(-1.0, tdamp); + MathExtra::add3(torque2on1, tdamp, torque2on1); + + // Damp twist + w1dotr = MathExtra::dot3(omega[i1],rhat); + w2dotr = MathExtra::dot3(omega[i2],rhat); + + MathExtra::scale3(w1dotr, rhat, wn1); + MathExtra::scale3(w2dotr, rhat, wn2); + + MathExtra::sub3(wn1, wn2, tdamp); // Eq. 38 + MathExtra::scale3(0.5*gtwist[type], tdamp); + MathExtra::add3(torque1on2, tdamp, torque1on2); + MathExtra::scale3(-1.0, tdamp); + MathExtra::add3(torque2on1, tdamp, torque2on1); + + // ------------------------------------------------------// + // Apply forces and torques to particles + // ------------------------------------------------------// + + if (newton_bond || i1 < nlocal) { + f[i1][0] -= force1on2[0]*smooth_sq; + f[i1][1] -= force1on2[1]*smooth_sq; + f[i1][2] -= force1on2[2]*smooth_sq; + + torque[i1][0] += torque2on1[0]*smooth_sq; + torque[i1][1] += torque2on1[1]*smooth_sq; + torque[i1][2] += torque2on1[2]*smooth_sq; + } + + if (newton_bond || i2 < nlocal) { + f[i2][0] += force1on2[0]*smooth_sq; + f[i2][1] += force1on2[1]*smooth_sq; + f[i2][2] += force1on2[2]*smooth_sq; + + torque[i2][0] += torque1on2[0]*smooth_sq; + torque[i2][1] += torque1on2[1]*smooth_sq; + torque[i2][2] += torque1on2[2]*smooth_sq; + } + + if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,Fr_mag,r[0],r[1],r[2]); + } +} + +/* ---------------------------------------------------------------------- */ + +void BondBPMRotational::allocate() +{ + allocated = 1; + int n = atom->nbondtypes; + + memory->create(Kr,n+1,"bond:Kr"); + memory->create(Ks,n+1,"bond:Ks"); + memory->create(Kt,n+1,"bond:Kt"); + memory->create(Kb,n+1,"bond:Kb"); + memory->create(Fcr,n+1,"bond:Fcr"); + memory->create(Fcs,n+1,"bond:Fcs"); + memory->create(Tct,n+1,"bond:Tct"); + memory->create(Tcb,n+1,"bond:Tcb"); + memory->create(gnorm,n+1,"bond:gnorm"); + memory->create(gslide,n+1,"bond:gslide"); + memory->create(groll,n+1,"bond:groll"); + memory->create(gtwist,n+1,"bond:gtwist"); + + memory->create(setflag,n+1,"bond:setflag"); + for (int i = 1; i <= n; i++) setflag[i] = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more types +------------------------------------------------------------------------- */ + +void BondBPMRotational::coeff(int narg, char **arg) +{ + if (narg != 13) error->all(FLERR,"Incorrect args for bond coefficients"); + if (!allocated) allocate(); + + int ilo,ihi; + utils::bounds(FLERR,arg[0],1,atom->nbondtypes,ilo,ihi,error); + + double Kr_one = utils::numeric(FLERR,arg[1],false,lmp); + double Ks_one = utils::numeric(FLERR,arg[2],false,lmp); + double Kt_one = utils::numeric(FLERR,arg[3],false,lmp); + double Kb_one = utils::numeric(FLERR,arg[4],false,lmp); + double Fcr_one = utils::numeric(FLERR,arg[5],false,lmp); + double Fcs_one = utils::numeric(FLERR,arg[6],false,lmp); + double Tct_one = utils::numeric(FLERR,arg[7],false,lmp); + double Tcb_one = utils::numeric(FLERR,arg[8],false,lmp); + double gnorm_one = utils::numeric(FLERR,arg[9],false,lmp); + double gslide_one = utils::numeric(FLERR,arg[10],false,lmp); + double groll_one = utils::numeric(FLERR,arg[11],false,lmp); + double gtwist_one = utils::numeric(FLERR,arg[12],false,lmp); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + Kr[i] = Kr_one; + Ks[i] = Ks_one; + Kt[i] = Kt_one; + Kb[i] = Kb_one; + Fcr[i] = Fcr_one; + Fcs[i] = Fcs_one; + Tct[i] = Tct_one; + Tcb[i] = Tcb_one; + gnorm[i] = gnorm_one; + gslide[i] = gslide_one; + groll[i] = groll_one; + gtwist[i] = gtwist_one; + setflag[i] = 1; + count++; + + if (Fcr[i]/Kr[i] > max_stretch) max_stretch = Fcr[i]/Kr[i]; + } + + if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients"); +} + +/* ---------------------------------------------------------------------- + check if pair defined and special_bond settings are valid +------------------------------------------------------------------------- */ + +void BondBPMRotational::init_style() +{ + BondBPM::init_style(); + + if (!atom->quat_flag || !atom->sphere_flag) + error->all(FLERR,"Bond bpm/rotational requires atom style sphere/bpm"); + if (comm->ghost_velocity == 0) + error->all(FLERR,"Bond bpm/rotational requires ghost atoms store velocity"); + + if(domain->dimension == 2) + error->warning(FLERR, "Bond style bpm/rotational not intended for 2d use"); + + if (!fix_bond_history) + fix_bond_history = (FixBondHistory *) modify->add_fix( + "BOND_HISTORY_BPM_ROTATIONAL all BOND_HISTORY 0 4"); +} + +/* ---------------------------------------------------------------------- + proc 0 writes out coeffs to restart file +------------------------------------------------------------------------- */ + +void BondBPMRotational::write_restart(FILE *fp) +{ + fwrite(&Kr[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Ks[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Kt[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Kb[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Fcr[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Fcs[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Tct[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&Tcb[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&gnorm[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&gslide[1],sizeof(double),atom->nbondtypes,fp); + fwrite(&groll[1],sizeof(double),atom->nbondtypes,fp); + fwrite(>wist[1],sizeof(double),atom->nbondtypes,fp); +} + +/* ---------------------------------------------------------------------- + proc 0 reads coeffs from restart file, bcasts them +------------------------------------------------------------------------- */ + +void BondBPMRotational::read_restart(FILE *fp) +{ + allocate(); + + if (comm->me == 0) { + utils::sfread(FLERR,&Kr[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Ks[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Kt[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Kb[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Fcr[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Fcs[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Tct[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&Tcb[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&gnorm[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&gslide[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,&groll[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + utils::sfread(FLERR,>wist[1],sizeof(double),atom->nbondtypes,fp,nullptr,error); + } + MPI_Bcast(&Kr[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Ks[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Kt[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Kb[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Fcr[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Fcs[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Tct[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&Tcb[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&gnorm[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&gslide[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(&groll[1],atom->nbondtypes,MPI_DOUBLE,0,world); + MPI_Bcast(>wist[1],atom->nbondtypes,MPI_DOUBLE,0,world); + + for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1; +} + +/* ---------------------------------------------------------------------- + proc 0 writes to data file +------------------------------------------------------------------------- */ + +void BondBPMRotational::write_data(FILE *fp) +{ + for (int i = 1; i <= atom->nbondtypes; i++) + fprintf(fp,"%d %g %g %g %g %g %g %g %g %g %g %g %g\n", + i,Kr[i],Ks[i],Kt[i],Kb[i],Fcr[i], Fcs[i], Tct[i], + Tcb[i], gnorm[i], gslide[i], groll[i], gtwist[i]); +} + +/* ---------------------------------------------------------------------- */ + +double BondBPMRotational::single(int type, double rsq, int i, int j, + double &fforce) +{ + // Not yet enabled + if (type <= 0) return 0.0; + + //double r0; + //for (int n = 0; n < atom->num_bond[i]; n ++) { + // if (atom->bond_atom[i][n] == atom->tag[j]) { + // r0 = fix_bond_history->get_atom_value(i, n, 0); + // } + //} +} diff --git a/src/BPM/bond_bpm_rotational.h b/src/BPM/bond_bpm_rotational.h new file mode 100644 index 0000000000..838226f0ff --- /dev/null +++ b/src/BPM/bond_bpm_rotational.h @@ -0,0 +1,77 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 BOND_CLASS +// clang-format off +BondStyle(bpm/rotational,BondBPMRotational) +// clang-format on +#else + +#ifndef LMP_BOND_BPM_ROTATIONAL_H +#define LMP_BOND_BPM_ROTATIONAL_H + +#include "bond_bpm.h" + +namespace LAMMPS_NS { + +class BondBPMRotational : public BondBPM { + public: + BondBPMRotational(class LAMMPS *); + virtual ~BondBPMRotational(); + virtual void compute(int, int); + void coeff(int, char **); + void init_style(); + void write_restart(FILE *); + void read_restart(FILE *); + void write_data(FILE *); + double single(int, double, int, int, double &); + + protected: + double *Kr, *Ks, *Kt, *Kb, *gnorm, *gslide, *groll, *gtwist; + double *Fcr, *Fcs, *Tct, *Tcb; + double acos_limit(double); + + void allocate(); + void store_data(); + double store_bond(int, int, int); + class FixBondHistory *fix_bond_history; +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Atom missing in BPM bond + +Bonded atom cannot be found + +E: Incorrect args for bond coefficients + +Self-explanatory. Check the input script or data file. + +E: Bond bpm/rotational requires atom style sphere/bpm + +Self-explanatory. + +E: Bond style bpm requires 1-3 and 1-4 special weights of 1.0 + +Self-explanatory. + +W: Bond style bpm/rotational not intended for 2d use, may be inefficient + +This bond style will perform a lot of unnecessary calculations in 2d + +*/ diff --git a/src/BPM/compute_nbond_atom.cpp b/src/BPM/compute_nbond_atom.cpp new file mode 100644 index 0000000000..ec39da9c5a --- /dev/null +++ b/src/BPM/compute_nbond_atom.cpp @@ -0,0 +1,147 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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. +------------------------------------------------------------------------- */ + +#include "compute_nbond_atom.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "update.h" + +#include + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeNBondAtom::ComputeNBondAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), + nbond(nullptr) +{ + if (narg < 3) error->all(FLERR,"Illegal compute nbond/atom command"); + + peratom_flag = 1; + size_peratom_cols = 0; + peatomflag = 1; + timeflag = 1; + comm_reverse = 1; + + nmax = 0; +} + +/* ---------------------------------------------------------------------- */ + +ComputeNBondAtom::~ComputeNBondAtom() +{ + memory->destroy(nbond); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeNBondAtom::compute_peratom() +{ + + invoked_peratom = update->ntimestep; + if (update->eflag_atom != invoked_peratom) + error->all(FLERR,"Per-atom nbond was not tallied on needed timestep"); + + // grow local nbond array if necessary + // needs to be atom->nmax in length + + if (atom->nmax > nmax) { + memory->destroy(nbond); + nmax = atom->nmax; + memory->create(nbond, nmax,"nbond/atom:nbond"); + vector_atom = nbond; + } + + // npair includes ghosts if either newton flag is set + // b/c some bonds/dihedrals call pair::ev_tally with pairwise info + // nbond includes ghosts if newton_bond is set + // ntotal includes ghosts if either newton flag is set + // KSpace includes ghosts if tip4pflag is set + + int nlocal = atom->nlocal; + tagint **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + + int ntotal = nlocal; + if (force->newton) ntotal += atom->nghost; + + // set local nbond array + int i, j, k; + int *num_bond = atom->num_bond; + int newton_bond = force->newton_bond; + + for (i = 0; i < ntotal; i++) nbond[i] = 0; + + for (i = 0; i < nlocal; i++) { + for (j = 0; j map(bond_atom[i][j]); + if (k < 0) continue; + + nbond[i] += 1; + if (newton_bond) nbond[k] += 1; + } + } + + // communicate ghost nbond between neighbor procs + if (force->newton) + comm->reverse_comm_compute(this); + + // zero nbond of atoms not in group + // only do this after comm since ghost contributions must be included + int *mask = atom->mask; + + for (i = 0; i < nlocal; i++) + if (!(mask[i] & groupbit)) nbond[i] = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +int ComputeNBondAtom::pack_reverse_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) buf[m++] = nbond[i]; + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeNBondAtom::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + nbond[j] += buf[m++]; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeNBondAtom::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/BPM/compute_nbond_atom.h b/src/BPM/compute_nbond_atom.h new file mode 100644 index 0000000000..3f92ef99c8 --- /dev/null +++ b/src/BPM/compute_nbond_atom.h @@ -0,0 +1,61 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 COMPUTE_CLASS +// clang-format off +ComputeStyle(nbond/atom,ComputeNBondAtom) +// clang-format on +#else + +#ifndef LMP_COMPUTE_NBOND_ATOM_H +#define LMP_COMPUTE_NBOND_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeNBondAtom : public Compute { + public: + ComputeNBondAtom(class LAMMPS *, int, char **); + ~ComputeNBondAtom(); + void init() {} + void compute_peratom(); + int pack_reverse_comm(int, int, double *); + void unpack_reverse_comm(int, int *, double *); + double memory_usage(); + + private: + int nmax; + double *nbond; +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Per-atom energy was not tallied on needed timestep + +You are using a thermo keyword that requires potentials to +have tallied energy, but they didn't on this timestep. See the +variable doc page for ideas on how to make this work. + +*/ diff --git a/src/BPM/fix_bond_history.cpp b/src/BPM/fix_bond_history.cpp new file mode 100644 index 0000000000..b0f046f43f --- /dev/null +++ b/src/BPM/fix_bond_history.cpp @@ -0,0 +1,294 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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. +------------------------------------------------------------------------- */ + +#include "fix_bond_history.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "group.h" +#include "memory.h" +#include "modify.h" +#include "neighbor.h" +#include "string.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define LB_FACTOR 1.5 +#define DELTA 10000 + +/* ---------------------------------------------------------------------- */ + +FixBondHistory::FixBondHistory(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg != 5) error->all(FLERR,"Illegal fix bond/history command"); + update_flag = utils::inumeric(FLERR,arg[3],false,lmp); + ndata = utils::inumeric(FLERR,arg[4],false,lmp); + nbond = atom->bond_per_atom; + + if (nbond == 0) + error->all(FLERR, "Cannot store bond variables without any bonds"); + + stored_flag = false; + restart_global = 1; + create_attribute = 1; + + bondstore = nullptr; + maxbond = 0; + allocate(); + + new_fix_id = nullptr; + array_id = nullptr; +} + +/* ---------------------------------------------------------------------- */ + +FixBondHistory::~FixBondHistory() +{ + if (new_fix_id && modify->nfix) modify->delete_fix(new_fix_id); + delete [] new_fix_id; + delete [] array_id; + + memory->destroy(bondstore); +} + +/* ---------------------------------------------------------------------- */ + +int FixBondHistory::setmask() +{ + int mask = 0; + mask |= PRE_EXCHANGE; + mask |= POST_NEIGHBOR; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixBondHistory::post_constructor() +{ + // Store saved bond quantities for each atom using fix property atom + char **newarg = new char*[5]; + int nvalue = 0; + int tmp1, tmp2; + + int nn = strlen(id) + strlen("_FIX_PROP_ATOM") + 1; + new_fix_id = new char[nn]; + strcpy(new_fix_id, "FIX_PROP_ATOM_"); + strcat(new_fix_id, id); + + nn = strlen(id) + 4 ; + array_id = new char[nn]; + strcpy(array_id, "d2_"); + strcat(array_id, id); + + char ncols[16]; + sprintf(ncols,"%d",nbond*ndata); + newarg[0] = new_fix_id; + newarg[1] = group->names[igroup]; + newarg[2] = (char *) "property/atom"; + newarg[3] = array_id; + newarg[4] = ncols; + + modify->add_fix(5,newarg); + index = atom->find_custom(&array_id[3],tmp1,tmp2); + + delete [] newarg; +} + +/* ---------------------------------------------------------------------- */ + +void FixBondHistory::update_atom_value(int i, int m, int idata, double value) +{ + if (idata >= ndata || m > nbond) error->all(FLERR, "Index exceeded in fix bond history"); + atom->darray[index][i][m*ndata+idata] = value; +} + +/* ---------------------------------------------------------------------- */ + +double FixBondHistory::get_atom_value(int i, int m, int idata) +{ + if (idata >= ndata || m > nbond) error->all(FLERR, "Index exceeded in fix bond history"); + return atom->darray[index][i][m*ndata+idata]; +} + +/* ---------------------------------------------------------------------- + If stored values are updated, need to copy to atom arrays before exchanging + Always called before neighborlist rebuilt + Also call prior to irregular communication in other fixes (e.g. deform) +------------------------------------------------------------------------- */ + +void FixBondHistory::pre_exchange() +{ + if (!update_flag) return; + + int i1, i2, n, m, idata; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + double **stored = atom->darray[index]; + + int nlocal = atom->nlocal; + tagint **bond_atom = atom->bond_atom; + int *num_bond = atom->num_bond; + tagint *tag = atom->tag; + + for (n = 0; n < nbondlist; n++) { + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + + // skip bond if already broken + if (bondlist[n][2] <= 0) { + continue; + } + + if (i1 < nlocal) { + for (m = 0; m < num_bond[i1]; m++) { + if (bond_atom[i1][m] == tag[i2]) { + for (idata = 0; idata < ndata; idata++) { + stored[i1][m*ndata+idata] = bondstore[n][idata]; + } + } + } + } + + if (i2 < nlocal) { + for (m = 0; m < num_bond[i2]; m++) { + if (bond_atom[i2][m] == tag[i1]) { + for (idata = 0; idata < ndata; idata++) { + stored[i1][m*ndata+idata] = bondstore[n][idata]; + } + } + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixBondHistory::allocate() +{ + //Ideally would just ask ntopo for maxbond, but protected + if (comm->nprocs == 1) maxbond = atom->nbonds; + else maxbond = static_cast (LB_FACTOR * atom->nbonds / comm->nprocs); + memory->create(bondstore,maxbond,ndata,"fix_bond_store:bondstore"); +} + +/* ---------------------------------------------------------------------- */ + +void FixBondHistory::setup_post_neighbor() +{ + post_neighbor(); +} + +/* ---------------------------------------------------------------------- + called after neighbor list is build + build array of stored bond quantities from fix property atom +------------------------------------------------------------------------- */ + +void FixBondHistory::post_neighbor() +{ + //Grow array if number of bonds has increased + while (neighbor->nbondlist >= maxbond) { + maxbond += DELTA; + memory->grow(bondstore,maxbond,ndata,"fix_bond_store:bondstore"); + } + + int i1, i2, n, m, idata; + int **bondlist = neighbor->bondlist; + int nbondlist = neighbor->nbondlist; + double **stored = atom->darray[index]; + + int nlocal = atom->nlocal; + tagint **bond_atom = atom->bond_atom; + int *num_bond = atom->num_bond; + tagint *tag = atom->tag; + + for (n = 0; n < nbondlist; n++) { + i1 = bondlist[n][0]; + i2 = bondlist[n][1]; + + // skip bond if already broken + if (bondlist[n][2] <= 0) { + continue; + } + + if (i1 < nlocal) { + for (m = 0; m < num_bond[i1]; m++) { + if (bond_atom[i1][m] == tag[i2]) { + for (idata = 0; idata < ndata; idata++) { + bondstore[n][idata] = stored[i1][m*ndata+idata]; + } + } + } + } + + if (i2 < nlocal){ + for (m = 0; m < num_bond[i2]; m++) { + if (bond_atom[i2][m] == tag[i1]) { + for (idata = 0; idata < ndata; idata++) { + bondstore[n][idata] = stored[i2][m*ndata+idata]; + } + } + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +double FixBondHistory::memory_usage() +{ + return maxbond * ndata * sizeof(double); +} + +/* ---------------------------------------------------------------------- */ + +void FixBondHistory::write_restart(FILE *fp) +{ + int n = 0; + double list[1]; + list[n++] = stored_flag; + + if (comm->me == 0) { + int size = n * sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(list,sizeof(double),n,fp); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixBondHistory::restart(char *buf) +{ + int n = 0; + double *list = (double *) buf; + stored_flag = static_cast (list[n++]); +} + +/* ---------------------------------------------------------------------- + initialize bond values to zero, called when atom is created +------------------------------------------------------------------------- */ + +void FixBondHistory::set_arrays(int i) +{ + double **stored = atom->darray[index]; + for (int m = 0; m < nbond; m++) + for (int idata = 0; idata < ndata; idata++) + stored[i][m*ndata+idata] = 0.0; +} \ No newline at end of file diff --git a/src/BPM/fix_bond_history.h b/src/BPM/fix_bond_history.h new file mode 100644 index 0000000000..0219e5de02 --- /dev/null +++ b/src/BPM/fix_bond_history.h @@ -0,0 +1,80 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 +// clang-format off +FixStyle(BOND_HISTORY,FixBondHistory) +// clang-format on +#else + +#ifndef LMP_FIX_BOND_HISTORY_H +#define LMP_FIX_BOND_HISTORY_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixBondHistory : public Fix { + public: + + FixBondHistory(class LAMMPS *, int, char **); + ~FixBondHistory(); + int setmask(); + void post_constructor(); + void setup_post_neighbor(); + void post_neighbor(); + void pre_exchange(); + double memory_usage(); + void write_restart(FILE *fp); + void restart(char *buf); + void set_arrays(int); + + void update_atom_value(int, int, int, double); + double get_atom_value(int, int, int); + double **bondstore; + int stored_flag; + + protected: + + void allocate(); + + int update_flag; + int nbond, maxbond, ndata; + int index; + char *new_fix_id; + char *array_id; +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +UNDOCUMENTED + +E: Index exceeded in fix bond history + +Bond requested non-existant data + +E: Cannot store bond variables without any bonds + +Atoms must have a nonzero number of bonds to store data + + + +*/ + diff --git a/src/BPM/fix_nve_sphere_bpm.cpp b/src/BPM/fix_nve_sphere_bpm.cpp new file mode 100644 index 0000000000..66801cc0b6 --- /dev/null +++ b/src/BPM/fix_nve_sphere_bpm.cpp @@ -0,0 +1,161 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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. +------------------------------------------------------------------------- */ + +#include "fix_nve_sphere_bpm.h" + +#include "atom.h" +#include "atom_vec.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "math_extra.h" + +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace MathExtra; + + +/* ---------------------------------------------------------------------- */ + +FixNVESphereBPM::FixNVESphereBPM(LAMMPS *lmp, int narg, char **arg) : + FixNVE(lmp, narg, arg) +{ + if (narg < 3) error->all(FLERR,"Illegal fix nve/sphere/bpm command"); + + time_integrate = 1; + + // process extra keywords + // inertia = moment of inertia prefactor for sphere or disc + + inertia = 0.4; + + int iarg = 3; + while (iarg < narg) { + if (strcmp(arg[iarg],"disc")==0) { + inertia = 0.5; + if (domain->dimension != 2) + error->all(FLERR,"Fix nve/sphere/bpm disc requires 2d simulation"); + iarg++; + } + else error->all(FLERR,"Illegal fix nve/sphere/bpm command"); + } + + inv_inertia = 1.0/inertia; + + // error checks + + if (!atom->quat_flag || !atom->sphere_flag) + error->all(FLERR,"Fix nve/sphere/bpm requires atom style sphere/bpm"); +} + +/* ---------------------------------------------------------------------- */ + +void FixNVESphereBPM::init() +{ + FixNVE::init(); + + // check that all particles are finite-size spheres + // no point particles allowed + + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) + if (radius[i] == 0.0) + error->one(FLERR,"Fix nve/sphere/bpm requires extended particles"); +} + +/* ---------------------------------------------------------------------- */ + +void FixNVESphereBPM::initial_integrate(int /*vflag*/) +{ + double dtq,dtfm,dtirotate,particle_inertia; + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **torque = atom->torque; + double **quat = atom->quat; + double *radius = atom->radius; + double *rmass = atom->rmass; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // set timestep here since dt may have changed or come via rRESPA + dtq = 0.5 * dtv; + + // update v,x,omega,quat for all particles + // d_omega/dt = torque / inertia + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + x[i][0] += dtv * v[i][0]; + x[i][1] += dtv * v[i][1]; + x[i][2] += dtv * v[i][2]; + + particle_inertia = inertia*(radius[i]*radius[i]*rmass[i]); + dtirotate = dtf / particle_inertia; + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + + MathExtra::richardson_sphere(quat[i],omega[i],dtq); + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNVESphereBPM::final_integrate() +{ + double dtfm,dtirotate,particle_inertia; + + double **v = atom->v; + double **f = atom->f; + double **omega = atom->omega; + double **torque = atom->torque; + double *rmass = atom->rmass; + double *radius = atom->radius; + int *mask = atom->mask; + int nlocal = atom->nlocal; + if (igroup == atom->firstgroup) nlocal = atom->nfirst; + + // update v,omega for all particles + // d_omega/dt = torque / inertia + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + dtfm = dtf / rmass[i]; + v[i][0] += dtfm * f[i][0]; + v[i][1] += dtfm * f[i][1]; + v[i][2] += dtfm * f[i][2]; + + particle_inertia = inertia*(radius[i]*radius[i]*rmass[i]); + dtirotate = dtf / particle_inertia; + omega[i][0] += dtirotate * torque[i][0]; + omega[i][1] += dtirotate * torque[i][1]; + omega[i][2] += dtirotate * torque[i][2]; + } +} diff --git a/src/BPM/fix_nve_sphere_bpm.h b/src/BPM/fix_nve_sphere_bpm.h new file mode 100644 index 0000000000..8e4d89ad67 --- /dev/null +++ b/src/BPM/fix_nve_sphere_bpm.h @@ -0,0 +1,71 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, 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 +// clang-format off +FixStyle(nve/sphere/bpm,FixNVESphereBPM) +// clang-format on +#else + +#ifndef LMP_FIX_NVE_SPHERE_BPM_H +#define LMP_FIX_NVE_SPHERE_BPM_H + +#include "fix_nve.h" + +namespace LAMMPS_NS { + +class FixNVESphereBPM : public FixNVE { + public: + FixNVESphereBPM(class LAMMPS *, int, char **); + virtual ~FixNVESphereBPM() {} + void init(); + virtual void initial_integrate(int); + virtual void final_integrate(); + + protected: + double inertia, inv_inertia; + int extra; + int dlm; +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Fix nve/sphere/bpm disc requires 2d simulation + +UNDOCUMENTED + +E: Fix nve/sphere/bpm requires atom style sphere/bpm + +Self-explanatory. + +E: Fix nve/sphere/bpm update dipole requires atom attribute mu + +An atom style with this attribute is needed. + +E: Fix nve/sphere/bpm requires extended particles + +This fix can only be used for particles of a finite size. + + +*/ diff --git a/src/BPM/fix_update_special_bonds.cpp b/src/BPM/fix_update_special_bonds.cpp new file mode 100644 index 0000000000..3fa9b84cd0 --- /dev/null +++ b/src/BPM/fix_update_special_bonds.cpp @@ -0,0 +1,235 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://lammps.sandia.gov/, Sandia National Laboratories + Steve Plimpton, sjplimp@sandia.gov + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#include "fix_update_special_bonds.h" + +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "pair.h" + +#include +#include +#include + +using namespace LAMMPS_NS; +using namespace FixConst; + +/* ---------------------------------------------------------------------- */ + +FixUpdateSpecialBonds::FixUpdateSpecialBonds(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg) +{ + if (narg != 3) error->all(FLERR,"Illegal fix update/special/bonds command"); + comm_forward = 1+atom->maxspecial; +} + +/* ---------------------------------------------------------------------- */ + +FixUpdateSpecialBonds::~FixUpdateSpecialBonds() +{ +} + +/* ---------------------------------------------------------------------- */ + +int FixUpdateSpecialBonds::setmask() +{ + int mask = 0; + mask |= PRE_EXCHANGE; + mask |= PRE_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixUpdateSpecialBonds::setup(int /*vflag*/) +{ + // Require atoms know about all of their bonds and if they break + if (force->newton_bond) + error->all(FLERR,"Fix update/special/bonds requires Newton bond off"); + + if (!atom->avec->bonds_allow) + error->all(FLERR,"Fix update/special/bonds requires atom bonds"); + + // special lj must be 0 1 1 to censor pair forces between bonded particles + // special coulomb must be 1 1 1 to ensure all pairs are included in the + // neighbor list and 1-3 and 1-4 special bond lists are skipped + if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || + force->special_lj[3] != 1.0) + error->all(FLERR,"Fix update/special/bonds requires special LJ weights = 0,1,1"); + if (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 || + force->special_coul[3] != 1.0) + error->all(FLERR,"Fix update/special/bonds requires special Coulomb weights = 1,1,1"); + + broken_pairs.clear(); +} + +/* ---------------------------------------------------------------------- + Update special bond list and atom bond arrays, empty broken bond list +------------------------------------------------------------------------- */ + +void FixUpdateSpecialBonds::pre_exchange() +{ + int i, j, key, m, n1, n3; + tagint min_tag, max_tag; + int nlocal = atom->nlocal; + + tagint *tag = atom->tag; + tagint *slist; + int **nspecial = atom->nspecial; + tagint **special = atom->special; + + for (auto const &key : broken_pairs) { + min_tag = key.first; + max_tag = key.second; + + i = atom->map(min_tag); + j = atom->map(max_tag); + + // remove i from special bond list for atom j and vice versa + if (i < nlocal) { + slist = special[i]; + n1 = nspecial[i][0]; + for (m = 0; m < n1; m++) + if (slist[m] == max_tag) break; + n3 = nspecial[i][2]; + for (; m < n3-1; m++) slist[m] = slist[m+1]; + nspecial[i][0]--; + nspecial[i][1]--; + nspecial[i][2]--; + } + + if (j < nlocal) { + slist = special[j]; + n1 = nspecial[j][0]; + for (m = 0; m < n1; m++) + if (slist[m] == min_tag) break; + n3 = nspecial[j][2]; + for (; m < n3-1; m++) slist[m] = slist[m+1]; + nspecial[j][0]--; + nspecial[j][1]--; + nspecial[j][2]--; + } + } + + // Forward updated special bond list + comm->forward_comm_fix(this); + + broken_pairs.clear(); +} + +/* ---------------------------------------------------------------------- + Loop neighbor list and update special bond lists for recently broken bonds +------------------------------------------------------------------------- */ + +void FixUpdateSpecialBonds::pre_force(int /*vflag*/) +{ + int i,j,n,m,ii,jj,inum,jnum; + int *ilist,*jlist,*numneigh,**firstneigh; + tagint min_tag, max_tag; + std::pair key; + + int **bond_type = atom->bond_type; + int *num_bond = atom->num_bond; + tagint **bond_atom = atom->bond_atom; + int nlocal = atom->nlocal; + + tagint *tag = atom->tag; + NeighList *list = force->pair->list; + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + min_tag = tag[i]; + max_tag = tag[j]; + if (max_tag < min_tag) { + min_tag = tag[j]; + max_tag = tag[i]; + } + key = std::make_pair(min_tag, max_tag); + + if (broken_pairs.find(key) != broken_pairs.end()) + jlist[jj] = j; // Clear special bond bits + } + } +} + +/* ---------------------------------------------------------------------- */ + +int FixUpdateSpecialBonds::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i,j,k,m,ns; + int **nspecial = atom->nspecial; + tagint **special = atom->special; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + ns = nspecial[j][0]; + buf[m++] = ubuf(ns).d; + for (k = 0; k < ns; k++) + buf[m++] = ubuf(special[j][k]).d; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void FixUpdateSpecialBonds::unpack_forward_comm(int n, int first, double *buf) +{ + int i,j,m,ns,last; + int **nspecial = atom->nspecial; + tagint **special = atom->special; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + ns = (int) ubuf(buf[m++]).i; + nspecial[i][0] = ns; + for (j = 0; j < ns; j++) + special[i][j] = (tagint) ubuf(buf[m++]).i; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixUpdateSpecialBonds::add_broken_bond(int i, int j) +{ + tagint *tag = atom->tag; + + tagint min_tag = tag[i]; + tagint max_tag = tag[j]; + if (max_tag < min_tag) { + min_tag = tag[j]; + max_tag = tag[i]; + } + std::pair key = std::make_pair(min_tag, max_tag); + + broken_pairs.insert(key); +} diff --git a/src/BPM/fix_update_special_bonds.h b/src/BPM/fix_update_special_bonds.h new file mode 100644 index 0000000000..ff47dd4d9b --- /dev/null +++ b/src/BPM/fix_update_special_bonds.h @@ -0,0 +1,76 @@ +/* -*- 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 +// clang-format off +FixStyle(update/special/bonds,FixUpdateSpecialBonds) +// clang-format on +#else + +#ifndef LMP_FIX_UPDATE_SPECIAL_BONDS_H +#define LMP_FIX_UPDATE_SPECIAL_BONDS_H + +#include "fix.h" + +#include +#include + +namespace LAMMPS_NS { + +class FixUpdateSpecialBonds : public Fix { + public: + FixUpdateSpecialBonds(class LAMMPS *, int, char **); + ~FixUpdateSpecialBonds(); + int setmask(); + void setup(int); + void pre_exchange(); + void pre_force(int); + int pack_forward_comm(int, int *, double *, int, int *); + void unpack_forward_comm(int, int, double *); + void add_broken_bond(int,int); + + protected: + std::set > broken_pairs; + inline int sbmask(int j) const { + return j >> SBBITS & 3; + } +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal fix censor/bonded/pairs command + +Self-explanatory. + +E: Fix censor/bonded/pairs requires Newton bond off + +Self-explanatory. + +E: Fix censor/bonded/pairs requires atom bonds + +Self-explanatory. + +E: Fix censor/bonded/pairs must be used without special bonds + +Self-explanatory. Look at the atom modify special command. + +E: Fix censor/bonded/pairs requires special_bonds = 0,0,0 + +Self-explanatory. + +*/ diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index 2fa0f5acb8..c45b535d6f 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -43,7 +43,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; - double radi,radj,radsum,rsq,r,rinv,rsqinv; + double radi,radj,radsum,rsq,r,rinv,rsqinv,factor_lj; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; @@ -112,8 +112,11 @@ void PairGranHertzHistory::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; + if (factor_lj == 0) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; @@ -244,6 +247,9 @@ void PairGranHertzHistory::compute(int eflag, int vflag) fx = delx*ccel + fs1; fy = dely*ccel + fs2; fz = delz*ccel + fs3; + fx *= factor_lj; + fy *= factor_lj; + fz *= factor_lj; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; @@ -251,6 +257,9 @@ void PairGranHertzHistory::compute(int eflag, int vflag) tor1 = rinv * (dely*fs3 - delz*fs2); tor2 = rinv * (delz*fs1 - delx*fs3); tor3 = rinv * (delx*fs2 - dely*fs1); + tor1 *= factor_lj; + tor2 *= factor_lj; + tor3 *= factor_lj; torque[i][0] -= radi*tor1; torque[i][1] -= radi*tor2; torque[i][2] -= radi*tor3; diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 1ef69e041d..11c867cd81 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -42,7 +42,7 @@ void PairGranHooke::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; - double radi,radj,radsum,rsq,r,rinv,rsqinv; + double radi,radj,radsum,rsq,r,rinv,rsqinv,factor_lj; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; @@ -101,8 +101,11 @@ void PairGranHooke::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; + if (factor_lj == 0) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; @@ -187,6 +190,9 @@ void PairGranHooke::compute(int eflag, int vflag) fx = delx*ccel + fs1; fy = dely*ccel + fs2; fz = delz*ccel + fs3; + fx *= factor_lj; + fy *= factor_lj; + fz *= factor_lj; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; @@ -194,6 +200,9 @@ void PairGranHooke::compute(int eflag, int vflag) tor1 = rinv * (dely*fs3 - delz*fs2); tor2 = rinv * (delz*fs1 - delx*fs3); tor3 = rinv * (delx*fs2 - dely*fs1); + tor1 *= factor_lj; + tor2 *= factor_lj; + tor3 *= factor_lj; torque[i][0] -= radi*tor1; torque[i][1] -= radi*tor2; torque[i][2] -= radi*tor3; diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index e6013b9940..9cd5f1085e 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -101,7 +101,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz; - double radi,radj,radsum,rsq,r,rinv,rsqinv; + double radi,radj,radsum,rsq,r,rinv,rsqinv,factor_lj; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; @@ -170,8 +170,11 @@ void PairGranHookeHistory::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; + if (factor_lj == 0) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; @@ -301,6 +304,9 @@ void PairGranHookeHistory::compute(int eflag, int vflag) fx = delx*ccel + fs1; fy = dely*ccel + fs2; fz = delz*ccel + fs3; + fx *= factor_lj; + fy *= factor_lj; + fz *= factor_lj; f[i][0] += fx; f[i][1] += fy; f[i][2] += fz; @@ -308,6 +314,9 @@ void PairGranHookeHistory::compute(int eflag, int vflag) tor1 = rinv * (dely*fs3 - delz*fs2); tor2 = rinv * (delz*fs1 - delx*fs3); tor3 = rinv * (delx*fs2 - dely*fs1); + tor1 *= factor_lj; + tor2 *= factor_lj; + tor3 *= factor_lj; torque[i][0] -= radi*tor1; torque[i][1] -= radi*tor2; torque[i][2] -= radi*tor3; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 2846403e4c..d82423b05e 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -150,7 +150,7 @@ void PairGranular::compute(int eflag, int vflag) { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz; - double radi,radj,radsum,rsq,r,rinv; + double radi,radj,radsum,rsq,r,rinv,factor_lj; double Reff, delta, dR, dR2, dist_to_contact; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; @@ -249,8 +249,11 @@ void PairGranular::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; j &= NEIGHMASK; + if (factor_lj == 0) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; @@ -651,6 +654,9 @@ void PairGranular::compute(int eflag, int vflag) fx = nx*Fntot + fs1; fy = ny*Fntot + fs2; fz = nz*Fntot + fs3; + fx *= factor_lj; + fy *= factor_lj; + fz *= factor_lj; f[i][0] += fx; f[i][1] += fy; @@ -659,6 +665,9 @@ void PairGranular::compute(int eflag, int vflag) tor1 = ny*fs3 - nz*fs2; tor2 = nz*fs1 - nx*fs3; tor3 = nx*fs2 - ny*fs1; + tor1 *= factor_lj; + tor2 *= factor_lj; + tor3 *= factor_lj; dist_to_contact = radi-0.5*delta; torque[i][0] -= dist_to_contact*tor1; @@ -666,9 +675,9 @@ void PairGranular::compute(int eflag, int vflag) torque[i][2] -= dist_to_contact*tor3; if (twist_model[itype][jtype] != TWIST_NONE) { - tortwist1 = magtortwist * nx; - tortwist2 = magtortwist * ny; - tortwist3 = magtortwist * nz; + tortwist1 = magtortwist * nx * factor_lj; + tortwist2 = magtortwist * ny * factor_lj; + tortwist3 = magtortwist * nz * factor_lj; torque[i][0] += tortwist1; torque[i][1] += tortwist2; @@ -676,9 +685,9 @@ void PairGranular::compute(int eflag, int vflag) } if (roll_model[itype][jtype] != ROLL_NONE) { - torroll1 = Reff*(ny*fr3 - nz*fr2); // n cross fr - torroll2 = Reff*(nz*fr1 - nx*fr3); - torroll3 = Reff*(nx*fr2 - ny*fr1); + torroll1 = Reff*(ny*fr3 - nz*fr2) * factor_lj; // n cross fr + torroll2 = Reff*(nz*fr1 - nx*fr3) * factor_lj; + torroll3 = Reff*(nx*fr2 - ny*fr1) * factor_lj; torque[i][0] += torroll1; torque[i][1] += torroll2; diff --git a/src/MISC/fix_pair_tracker.cpp b/src/MISC/fix_pair_tracker.cpp deleted file mode 100644 index 755025baed..0000000000 --- a/src/MISC/fix_pair_tracker.cpp +++ /dev/null @@ -1,369 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/ 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. -------------------------------------------------------------------------- */ - -#include "fix_pair_tracker.h" -#include "atom.h" -#include "atom_vec.h" -#include "error.h" -#include "group.h" -#include "memory.h" -#include "modify.h" -#include "tokenizer.h" -#include "update.h" - -#include - -using namespace LAMMPS_NS; -using namespace FixConst; - -#define DELTA 1000 - -/* ---------------------------------------------------------------------- */ - -FixPairTracker::FixPairTracker(LAMMPS *lmp, int narg, char **arg) : - Fix(lmp, narg, arg), nvalues(0), vector(NULL), array(NULL), pack_choice(NULL) -{ - if (narg < 3) error->all(FLERR, "Illegal fix pair/tracker command"); - local_flag = 1; - - nevery = utils::inumeric(FLERR, arg[3], false, lmp); - if (nevery <= 0) error->all(FLERR, "Illegal fix pair/tracker command"); - local_freq = nevery; - - // If optional arguments included, this will be oversized - nvalues = narg - 4; - pack_choice = new FnPtrPack[nvalues]; - - tmin = -1; - type_filter = nullptr; - int iarg = 4; - nvalues = 0; - while (iarg < narg) { - if (strcmp(arg[iarg], "id1") == 0) { - pack_choice[nvalues++] = &FixPairTracker::pack_id1; - } else if (strcmp(arg[iarg], "id2") == 0) { - pack_choice[nvalues++] = &FixPairTracker::pack_id2; - - } else if (strcmp(arg[iarg], "time/created") == 0) { - pack_choice[nvalues++] = &FixPairTracker::pack_time_created; - } else if (strcmp(arg[iarg], "time/broken") == 0) { - pack_choice[nvalues++] = &FixPairTracker::pack_time_broken; - } else if (strcmp(arg[iarg], "time/total") == 0) { - pack_choice[nvalues++] = &FixPairTracker::pack_time_total; - - } else if (strcmp(arg[iarg], "x") == 0) { - pack_choice[nvalues++] = &FixPairTracker::pack_x; - } else if (strcmp(arg[iarg], "y") == 0) { - pack_choice[nvalues++] = &FixPairTracker::pack_y; - } else if (strcmp(arg[iarg], "z") == 0) { - pack_choice[nvalues++] = &FixPairTracker::pack_z; - - } else if (strcmp(arg[iarg], "r/min") == 0) { - pack_choice[nvalues++] = &FixPairTracker::pack_rmin; - } else if (strcmp(arg[iarg], "r/ave") == 0) { - pack_choice[nvalues++] = &FixPairTracker::pack_rave; - - } else if (strcmp(arg[iarg], "time/min") == 0) { - if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in fix pair/tracker command"); - tmin = utils::numeric(FLERR, arg[iarg + 1], false, lmp); - iarg++; - - } else if (strcmp(arg[iarg], "type/include") == 0) { - if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in fix pair/tracker command"); - int ntypes = atom->ntypes; - - int i, j, itype, jtype, in, jn, infield, jnfield; - int inlo, inhi, jnlo, jnhi; - char *istr, *jstr; - if (!type_filter) { - memory->create(type_filter, ntypes + 1, ntypes + 1, "fix/pair/tracker:type_filter"); - - for (i = 0; i <= ntypes; i++) { - for (j = 0; j <= ntypes; j++) { type_filter[i][j] = 0; } - } - } - - in = strlen(arg[iarg + 1]) + 1; - istr = new char[in]; - strcpy(istr, arg[iarg + 1]); - std::vector iwords = Tokenizer(istr, ",").as_vector(); - infield = iwords.size(); - - jn = strlen(arg[iarg + 2]) + 1; - jstr = new char[jn]; - strcpy(jstr, arg[iarg + 2]); - std::vector jwords = Tokenizer(jstr, ",").as_vector(); - jnfield = jwords.size(); - - for (i = 0; i < infield; i++) { - const char *ifield = iwords[i].c_str(); - utils::bounds(FLERR, ifield, 1, ntypes, inlo, inhi, error); - - for (j = 0; j < jnfield; j++) { - const char *jfield = jwords[j].c_str(); - utils::bounds(FLERR, jfield, 1, ntypes, jnlo, jnhi, error); - - for (itype = inlo; itype <= inhi; itype++) { - for (jtype = jnlo; jtype <= jnhi; jtype++) { - type_filter[itype][jtype] = 1; - type_filter[jtype][itype] = 1; - } - } - } - } - - delete[] istr; - delete[] jstr; - - iarg += 2; - - } else - error->all(FLERR, "Invalid keyword in fix pair/tracker command"); - - iarg++; - } - - if (nvalues == 1) - size_local_cols = 0; - else - size_local_cols = nvalues; - - nmax = 0; - ncount = 0; - vector = NULL; - array = NULL; -} - -/* ---------------------------------------------------------------------- */ - -FixPairTracker::~FixPairTracker() -{ - delete[] pack_choice; - - memory->destroy(vector); - memory->destroy(array); - memory->destroy(type_filter); -} - -/* ---------------------------------------------------------------------- */ - -int FixPairTracker::setmask() -{ - int mask = 0; - mask |= POST_FORCE; - return mask; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::init() -{ - // Set size of array/vector - ncount = 0; - - if (ncount > nmax) reallocate(ncount); - - size_local_rows = ncount; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::lost_contact(int i, int j, double time_tmp, double nstep_tmp, double rsum_tmp, - double rmin_tmp) -{ - - double time = update->atime + (update->ntimestep - update->atimestep) * update->dt; - if ((time - time_tmp) < tmin) return; - - if (type_filter) { - int *type = atom->type; - if (type_filter[type[i]][type[j]] == 0) return; - } - - int *mask = atom->mask; - if (!(mask[i] & groupbit)) return; - if (!(mask[j] & groupbit)) return; - - if (ncount == nmax) reallocate(ncount); - - index_i = i; - index_j = j; - - rmin = rmin_tmp; - rsum = rsum_tmp; - time_initial = time_tmp; - nstep_initial = nstep_tmp; - - // fill vector or array with local values - if (nvalues == 1) { - (this->*pack_choice[0])(0); - } else { - for (int k = 0; k < nvalues; k++) { (this->*pack_choice[k])(k); } - } - - ncount += 1; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::post_force(int /*vflag*/) -{ - if (update->ntimestep % nevery == 0) { - size_local_rows = ncount; - ncount = 0; - } -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::reallocate(int n) -{ - // grow vector or array - while (nmax <= n) nmax += DELTA; - - if (nvalues == 1) { - memory->grow(vector, nmax, "fix_pair_tracker:vector"); - vector_local = vector; - } else { - memory->grow(array, nmax, nvalues, "fix_pair_tracker:array"); - array_local = array; - } -} - -/* ---------------------------------------------------------------------- - memory usage of local data -------------------------------------------------------------------------- */ - -double FixPairTracker::memory_usage() -{ - double bytes = nmax * (double) nvalues * sizeof(double); - bytes += nmax * 2 * sizeof(int); - return bytes; -} - -/* ---------------------------------------------------------------------- - one method for every keyword fix pair/tracker can output - the atom property is packed into a local vector or array -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::pack_time_created(int n) -{ - if (nvalues == 1) - vector[ncount] = time_initial; - else - array[ncount][n] = time_initial; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::pack_time_broken(int n) -{ - double time = update->atime + (update->ntimestep - update->atimestep) * update->dt; - if (nvalues == 1) - vector[ncount] = time; - else - array[ncount][n] = time; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::pack_time_total(int n) -{ - double time = update->atime + (update->ntimestep - update->atimestep) * update->dt; - if (nvalues == 1) - vector[ncount] = time - time_initial; - else - array[ncount][n] = time - time_initial; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::pack_id1(int n) -{ - tagint *tag = atom->tag; - - if (nvalues == 1) - vector[ncount] = tag[index_i]; - else - array[ncount][n] = tag[index_i]; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::pack_id2(int n) -{ - tagint *tag = atom->tag; - - if (nvalues == 1) - vector[ncount] = tag[index_j]; - else - array[ncount][n] = tag[index_j]; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::pack_x(int n) -{ - double **x = atom->x; - - if (nvalues == 1) - vector[ncount] = (x[index_i][0] + x[index_j][0]) / 2; - else - array[ncount][n] = (x[index_i][0] + x[index_j][0]) / 2; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::pack_y(int n) -{ - double **x = atom->x; - - if (nvalues == 1) - vector[ncount] = (x[index_i][1] + x[index_j][1]) / 2; - else - array[ncount][n] = (x[index_i][1] + x[index_j][1]) / 2; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::pack_z(int n) -{ - double **x = atom->x; - - if (nvalues == 1) - vector[ncount] = (x[index_i][2] + x[index_j][2]) / 2; - else - array[ncount][n] = (x[index_i][2] + x[index_j][2]) / 2; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::pack_rmin(int n) -{ - if (nvalues == 1) - vector[ncount] = rmin; - else - array[ncount][n] = rmin; -} - -/* ---------------------------------------------------------------------- */ - -void FixPairTracker::pack_rave(int n) -{ - if (nvalues == 1) - vector[ncount] = rsum / (update->ntimestep - nstep_initial); - else - array[ncount][n] = rsum / (update->ntimestep - nstep_initial); -} diff --git a/src/MISC/pair_tracker.cpp b/src/MISC/pair_tracker.cpp index f2fdb71081..d98d462ab6 100644 --- a/src/MISC/pair_tracker.cpp +++ b/src/MISC/pair_tracker.cpp @@ -19,20 +19,23 @@ #include "fix.h" #include "fix_dummy.h" #include "fix_neigh_history.h" -#include "fix_pair_tracker.h" +#include "fix_store_local.h" #include "force.h" #include "memory.h" #include "modify.h" #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" +#include "tokenizer.h" +#include "utils.h" #include "update.h" using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp) +PairTracker::PairTracker(LAMMPS *lmp) : + Pair(lmp), pack_choice(NULL) { single_enable = 1; no_virial_fdotr_compute = 1; @@ -43,6 +46,7 @@ PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp) nondefault_history_transfer = 1; finitecutflag = 0; + tmin = -1; // create dummy fix as placeholder for FixNeighHistory // this is so final order of Modify:fix will conform to input script @@ -50,6 +54,12 @@ PairTracker::PairTracker(LAMMPS *lmp) : Pair(lmp) fix_history = nullptr; modify->add_fix("NEIGH_HISTORY_TRACK_DUMMY all DUMMY"); fix_dummy = (FixDummy *) modify->fix[modify->nfix - 1]; + + fix_store_local = nullptr; + + output_data = nullptr; + pack_choice = nullptr; + type_filter = nullptr; } /* ---------------------------------------------------------------------- */ @@ -71,6 +81,12 @@ PairTracker::~PairTracker() delete[] maxrad_dynamic; delete[] maxrad_frozen; } + + delete[] pack_choice; + + delete [] id_fix_store_local; + memory->destroy(output_data); + memory->destroy(type_filter); } /* ---------------------------------------------------------------------- */ @@ -130,8 +146,9 @@ void PairTracker::compute(int eflag, int vflag) data = &alldata[size_history * jj]; if (touch[jj] == 1) { - fix_pair_tracker->lost_contact(i, j, data[0], data[1], data[2], data[3]); + process_data(i, j, data); } + touch[jj] = 0; data[0] = 0.0; // initial time data[1] = 0.0; // initial timestep @@ -159,7 +176,7 @@ void PairTracker::compute(int eflag, int vflag) data = &alldata[size_history * jj]; if (touch[jj] == 1) { - fix_pair_tracker->lost_contact(i, j, data[0], data[1], data[2], data[3]); + process_data(i, j, data); } touch[jj] = 0; @@ -216,14 +233,98 @@ void PairTracker::allocate() void PairTracker::settings(int narg, char **arg) { - if (narg != 0 && narg != 1) error->all(FLERR, "Illegal pair_style command"); + if (narg < 1) error->all(FLERR, "Illegal pair_style command"); - if (narg == 1) { - if (strcmp(arg[0], "finite") == 0) + id_fix_store_local = utils::strdup(arg[0]); + + // If optional arguments included, this will be oversized + pack_choice = new FnPtrPack[narg - 1]; + + nvalues = 0; + int iarg = 1; + while (iarg < narg) { + if (strcmp(arg[iarg], "finite") == 0) { finitecutflag = 1; - else + } else if (strcmp(arg[iarg], "id1") == 0) { + pack_choice[nvalues++] = &PairTracker::pack_id1; + } else if (strcmp(arg[iarg], "id2") == 0) { + pack_choice[nvalues++] = &PairTracker::pack_id2; + } else if (strcmp(arg[iarg], "time/created") == 0) { + pack_choice[nvalues++] = &PairTracker::pack_time_created; + } else if (strcmp(arg[iarg], "time/broken") == 0) { + pack_choice[nvalues++] = &PairTracker::pack_time_broken; + } else if (strcmp(arg[iarg], "time/total") == 0) { + pack_choice[nvalues++] = &PairTracker::pack_time_total; + } else if (strcmp(arg[iarg], "x") == 0) { + pack_choice[nvalues++] = &PairTracker::pack_x; + } else if (strcmp(arg[iarg], "y") == 0) { + pack_choice[nvalues++] = &PairTracker::pack_y; + } else if (strcmp(arg[iarg], "z") == 0) { + pack_choice[nvalues++] = &PairTracker::pack_z; + } else if (strcmp(arg[iarg], "r/min") == 0) { + pack_choice[nvalues++] = &PairTracker::pack_rmin; + } else if (strcmp(arg[iarg], "r/ave") == 0) { + pack_choice[nvalues++] = &PairTracker::pack_rave; + } else if (strcmp(arg[iarg], "time/min") == 0) { + if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in pair tracker command"); + tmin = utils::numeric(FLERR, arg[iarg + 1], false, lmp); + iarg++; + + } else if (strcmp(arg[iarg], "type/include") == 0) { + if (iarg + 1 >= narg) error->all(FLERR, "Invalid keyword in pair tracker command"); + int ntypes = atom->ntypes; + + int i, j, itype, jtype, in, jn, infield, jnfield; + int inlo, inhi, jnlo, jnhi; + char *istr, *jstr; + if (!type_filter) { + memory->create(type_filter, ntypes + 1, ntypes + 1, "pair/tracker:type_filter"); + + for (i = 0; i <= ntypes; i++) { + for (j = 0; j <= ntypes; j++) type_filter[i][j] = 0; + } + } + + in = strlen(arg[iarg + 1]) + 1; + istr = new char[in]; + strcpy(istr, arg[iarg + 1]); + std::vector iwords = Tokenizer(istr, ",").as_vector(); + infield = iwords.size(); + + jn = strlen(arg[iarg + 2]) + 1; + jstr = new char[jn]; + strcpy(jstr, arg[iarg + 2]); + std::vector jwords = Tokenizer(jstr, ",").as_vector(); + jnfield = jwords.size(); + + for (i = 0; i < infield; i++) { + const char *ifield = iwords[i].c_str(); + utils::bounds(FLERR, ifield, 1, ntypes, inlo, inhi, error); + + for (j = 0; j < jnfield; j++) { + const char *jfield = jwords[j].c_str(); + utils::bounds(FLERR, jfield, 1, ntypes, jnlo, jnhi, error); + + for (itype = inlo; itype <= inhi; itype++) { + for (jtype = jnlo; jtype <= jnhi; jtype++) { + type_filter[itype][jtype] = 1; + type_filter[jtype][itype] = 1; + } + } + } + } + delete[] istr; + delete[] jstr; + iarg += 2; + + } else { error->all(FLERR, "Illegal pair_style command"); + } + iarg ++; } + + if (nvalues == 0) error->all(FLERR, "Must request at least one value to output"); + memory->create(output_data, nvalues, "pair/tracker:output_data"); } /* ---------------------------------------------------------------------- @@ -262,7 +363,6 @@ void PairTracker::coeff(int narg, char **arg) void PairTracker::init_style() { int i; - // error and warning checks if (!atom->radius_flag && finitecutflag) @@ -275,7 +375,7 @@ void PairTracker::init_style() neighbor->requests[irequest]->size = 1; neighbor->requests[irequest]->history = 1; // history flag won't affect results, but match granular pairstyles - // so neighborlist can be copied to reduce overhead + // to copy neighbor list and reduce overhead } // if history is stored and first init, create Fix to store history @@ -284,7 +384,7 @@ void PairTracker::init_style() if (fix_history == nullptr) { modify->replace_fix("NEIGH_HISTORY_TRACK_DUMMY", - fmt::format("NEIGH_HISTORY_TRACK all NEIGH_HISTORY {}", size_history), 1); + fmt::format("NEIGH_HISTORY_TRACK all NEIGH_HISTORY {}", size_history), 1); int ifix = modify->find_fix("NEIGH_HISTORY_TRACK"); fix_history = (FixNeighHistory *) modify->fix[ifix]; fix_history->pair = this; @@ -295,9 +395,14 @@ void PairTracker::init_style() if (force->pair->beyond_contact) error->all(FLERR, - "Pair tracker incompatible with granular pairstyles that extend beyond contact"); - // check for FixPour and FixDeposit so can extract particle radii + "Pair tracker incompatible with granular pairstyles that extend beyond contact"); + // check for FixFreeze and set freeze_group_bit + int ifix = modify->find_fix_by_style("^freeze"); + if (ifix < 0) freeze_group_bit = 0; + else freeze_group_bit = modify->fix[ifix]->groupbit; + + // check for FixPour and FixDeposit so can extract particle radii int ipour; for (ipour = 0; ipour < modify->nfix; ipour++) if (strcmp(modify->fix[ipour]->style, "pour") == 0) break; @@ -310,7 +415,6 @@ void PairTracker::init_style() // set maxrad_dynamic and maxrad_frozen for each type // include future FixPour and FixDeposit particles as dynamic - int itype; for (i = 1; i <= atom->ntypes; i++) { onerad_dynamic[i] = onerad_frozen[i] = 0.0; @@ -343,9 +447,13 @@ void PairTracker::init_style() if (ifix < 0) error->all(FLERR, "Could not find pair fix neigh history ID"); fix_history = (FixNeighHistory *) modify->fix[ifix]; - ifix = modify->find_fix_by_style("pair/tracker"); - if (ifix < 0) error->all(FLERR, "Cannot use pair tracker without fix pair/tracker"); - fix_pair_tracker = (FixPairTracker *) modify->fix[ifix]; + ifix = modify->find_fix(id_fix_store_local); + if (ifix < 0) error->all(FLERR, "Cannot find fix store/local"); + if (strcmp(modify->fix[ifix]->style, "store/local") != 0) + error->all(FLERR, "Incorrect fix style matched, not store/local"); + fix_store_local = (FixStoreLocal *) modify->fix[ifix]; + if (fix_store_local->nvalues != nvalues) + error->all(FLERR, "Inconsistent number of output variables in fix store/local"); } /* ---------------------------------------------------------------------- @@ -357,7 +465,6 @@ double PairTracker::init_one(int i, int j) if (!allocated) allocate(); // always mix prefactors geometrically - if (setflag[i][j] == 0) { cut[i][j] = mix_distance(cut[i][i], cut[j][j]); } cut[j][i] = cut[i][j]; @@ -469,3 +576,105 @@ double PairTracker::radii2cut(double r1, double r2) double cut = r1 + r2; return cut; } + +/* ---------------------------------------------------------------------- */ + +void PairTracker::process_data(int i, int j, double * input_data) +{ + double time = update->atime + (update->ntimestep - update->atimestep) * update->dt; + int time_initial = (int) input_data[0]; + if ((time - time_initial) < tmin) return; + + if (type_filter) { + int *type = atom->type; + if (type_filter[type[i]][type[j]] == 0) return; + } + + for (int k = 0; k < nvalues; k++) (this->*pack_choice[k])(k, i, j, input_data); + fix_store_local->add_data(output_data, i, j); +} + +/* ---------------------------------------------------------------------- + one method for every keyword fix pair/tracker can output + the atom property is packed into a local vector or array +------------------------------------------------------------------------- */ + +void PairTracker::pack_time_created(int n, int i, int j, double * data) +{ + double time_initial = data[0]; + output_data[n] = time_initial; +} + +/* ---------------------------------------------------------------------- */ + +void PairTracker::pack_time_broken(int n, int i, int j, double * data) +{ + double time = update->atime + (update->ntimestep - update->atimestep) * update->dt; + output_data[n] = time; +} + +/* ---------------------------------------------------------------------- */ + +void PairTracker::pack_time_total(int n, int i, int j, double * data) +{ + double time = update->atime + (update->ntimestep - update->atimestep) * update->dt; + double time_initial = data[0]; + output_data[n] = time - time_initial; +} + +/* ---------------------------------------------------------------------- */ + +void PairTracker::pack_id1(int n, int i, int j, double * data) +{ + tagint *tag = atom->tag; + output_data[n] = tag[i]; +} + +/* ---------------------------------------------------------------------- */ + +void PairTracker::pack_id2(int n, int i, int j, double * data) +{ + tagint *tag = atom->tag; + output_data[n] = tag[j]; +} + +/* ---------------------------------------------------------------------- */ + +void PairTracker::pack_x(int n, int i, int j, double * data) +{ + double **x = atom->x; + output_data[n] = (x[i][0] + x[j][0])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void PairTracker::pack_y(int n, int i, int j, double * data) +{ + double **x = atom->x; + output_data[n] = (x[i][1] + x[j][1])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void PairTracker::pack_z(int n, int i, int j, double * data) +{ + double **x = atom->x; + output_data[n] = (x[i][2] + x[j][2])*0.5; +} + +/* ---------------------------------------------------------------------- */ + +void PairTracker::pack_rmin(int n, int i, int j, double * data) +{ + double rmin = data[3]; + output_data[n] = rmin; +} + +/* ---------------------------------------------------------------------- */ + +void PairTracker::pack_rave(int n, int i, int j, double * data) +{ + double rsum = data[2]; + double nstep_initial = data[1]; + output_data[n] = rsum / (update->ntimestep - nstep_initial); +} \ No newline at end of file diff --git a/src/MISC/pair_tracker.h b/src/MISC/pair_tracker.h index c6825c410e..3cb9344212 100644 --- a/src/MISC/pair_tracker.h +++ b/src/MISC/pair_tracker.h @@ -51,11 +51,33 @@ class PairTracker : public Pair { double *maxrad_dynamic, *maxrad_frozen; int freeze_group_bit; + char *id_fix_store_local; class FixDummy *fix_dummy; class FixNeighHistory *fix_history; - class FixPairTracker *fix_pair_tracker; + class FixStoreLocal *fix_store_local; + int **type_filter; + double tmin; + + int nvalues, ncount; + double *output_data; + typedef void (PairTracker::*FnPtrPack)(int, int, int, double *); + FnPtrPack *pack_choice; // ptrs to pack functions + + void pack_id1(int, int, int, double *); + void pack_id2(int, int, int, double *); + void pack_time_created(int, int, int, double *); + void pack_time_broken(int, int, int, double *); + void pack_time_total(int, int, int, double *); + void pack_x(int, int, int, double *); + void pack_y(int, int, int, double *); + void pack_z(int, int, int, double *); + void pack_rmin(int, int, int, double *); + void pack_rave(int, int, int, double *); + + void transfer_data(int, int); void transfer_history(double *, double *); + void process_data(int, int, double *); void allocate(); }; @@ -72,20 +94,37 @@ Self-explanatory. Check the input script syntax and compare to the documentation for the command. You can use -echo screen as a command-line option when running LAMMPS to see the offending line. +E: Invalid keyword in pair tracker command + +Self-explanatory. + E: Incorrect args for pair coefficients Self-explanatory. Check the input script or data file. +E: Must request at least one value to output + +Must include at least one bond property to store in fix store/local + E: Pair tracker requires atom attribute radius for finite cutoffs The atom style defined does not have these attributes. +E: Pair tracker incompatible with granular pairstyles that extend beyond contact + +Self-explanatory. + E: Could not find pair fix neigh history ID The associated fix neigh/history is missing -E: Cannot use pair tracker without fix pair/tracker +E: Cannot use pair tracker without fix store/local -This pairstyle requires one to define a pair/tracker fix +The associated fix store/local does not exist + +E: Inconsistent number of output variables in fix store/local + +The number of values specified in fix store/local disagrees with +the number of values requested in pair tracker */ diff --git a/src/MOLECULE/bond_quartic.cpp b/src/MOLECULE/bond_quartic.cpp index b3eea79064..40e2583e63 100644 --- a/src/MOLECULE/bond_quartic.cpp +++ b/src/MOLECULE/bond_quartic.cpp @@ -34,6 +34,7 @@ using namespace LAMMPS_NS; BondQuartic::BondQuartic(LAMMPS *lmp) : Bond(lmp) { + partial_flag = 1; TWO_1_3 = pow(2.0,(1.0/3.0)); } diff --git a/src/Makefile b/src/Makefile index 1527adfd07..015a96775a 100644 --- a/src/Makefile +++ b/src/Makefile @@ -53,6 +53,7 @@ PACKAGE = \ awpmd \ bocs \ body \ + bpm \ brownian \ cg-dna \ cg-sdk \ @@ -147,7 +148,8 @@ PACKMOST = \ asphere \ bocs \ body \ - brownian \ + bpm + brownian \ cg-dna \ cg-sdk \ class2 \ diff --git a/src/OPENMP/fix_neigh_history_omp.cpp b/src/OPENMP/fix_neigh_history_omp.cpp index 50599b2669..086ddb1c4f 100644 --- a/src/OPENMP/fix_neigh_history_omp.cpp +++ b/src/OPENMP/fix_neigh_history_omp.cpp @@ -576,7 +576,7 @@ void FixNeighHistoryOMP::post_neighbor() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - rflag = sbmask(j); + rflag = histmask(j); j &= NEIGHMASK; jlist[jj] = j; diff --git a/src/OPENMP/npair_half_size_bin_newtoff_omp.cpp b/src/OPENMP/npair_half_size_bin_newtoff_omp.cpp index 445c2ff286..07082868f6 100644 --- a/src/OPENMP/npair_half_size_bin_newtoff_omp.cpp +++ b/src/OPENMP/npair_half_size_bin_newtoff_omp.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_bin_newtoff_omp.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neigh_list.h" #include "npair_omp.h" @@ -40,8 +43,10 @@ NPairHalfSizeBinNewtoffOmp::NPairHalfSizeBinNewtoffOmp(LAMMPS *lmp) : void NPairHalfSizeBinNewtoffOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; @@ -50,7 +55,8 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,ibin; + int i,j,jh,k,n,ibin,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -61,7 +67,14 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -81,6 +94,11 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list) ztmp = x[i][2]; radi = radius[i]; ibin = atom2bin[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in surrounding bins in stencil including self // only store pair if i < j @@ -100,10 +118,23 @@ void NPairHalfSizeBinNewtoffOmp::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/OPENMP/npair_half_size_bin_newton_omp.cpp b/src/OPENMP/npair_half_size_bin_newton_omp.cpp index be110f07e0..42a74f13aa 100644 --- a/src/OPENMP/npair_half_size_bin_newton_omp.cpp +++ b/src/OPENMP/npair_half_size_bin_newton_omp.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_bin_newton_omp.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neigh_list.h" #include "npair_omp.h" @@ -39,8 +42,10 @@ NPairHalfSizeBinNewtonOmp::NPairHalfSizeBinNewtonOmp(LAMMPS *lmp) : void NPairHalfSizeBinNewtonOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; @@ -49,7 +54,8 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,ibin; + int i,j,jh,k,n,ibin,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -58,7 +64,14 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -77,6 +90,11 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list @@ -101,10 +119,23 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } } @@ -123,10 +154,23 @@ void NPairHalfSizeBinNewtonOmp::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/OPENMP/npair_half_size_bin_newton_tri_omp.cpp b/src/OPENMP/npair_half_size_bin_newton_tri_omp.cpp index bc3bda1181..5a77de0796 100644 --- a/src/OPENMP/npair_half_size_bin_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_size_bin_newton_tri_omp.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_bin_newton_tri_omp.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neigh_list.h" #include "npair_omp.h" @@ -39,8 +42,10 @@ NPairHalfSizeBinNewtonTriOmp::NPairHalfSizeBinNewtonTriOmp(LAMMPS *lmp) : void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -48,7 +53,8 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,ibin; + int i,j,jh,k,n,ibin,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -59,7 +65,14 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -78,6 +91,11 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in bins in stencil // pairs for atoms j "below" i are excluded @@ -107,10 +125,23 @@ void NPairHalfSizeBinNewtonTriOmp::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/OPENMP/npair_half_size_multi_newtoff_omp.cpp b/src/OPENMP/npair_half_size_multi_newtoff_omp.cpp index da2d7a7590..d58d250209 100644 --- a/src/OPENMP/npair_half_size_multi_newtoff_omp.cpp +++ b/src/OPENMP/npair_half_size_multi_newtoff_omp.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_multi_newtoff_omp.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neighbor.h" #include "neigh_list.h" @@ -41,8 +44,10 @@ NPairHalfSizeMultiNewtoffOmp::NPairHalfSizeMultiNewtoffOmp(LAMMPS *lmp) : NPair( void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -50,7 +55,9 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns; + int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns; + int which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -63,7 +70,14 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -84,6 +98,11 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } ibin = atom2bin[i]; @@ -104,27 +123,40 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list) ns = nstencil_multi[icollection][jcollection]; for (k = 0; k < ns; k++) { - js = binhead_multi[jcollection][jbin + s[k]]; - for (j = js; j >=0; j = bins[j]) { - if (j <= i) continue; + js = binhead_multi[jcollection][jbin + s[k]]; + for (j = js; j >=0; j = bins[j]) { + if (j <= i) continue; jtype = type[j]; - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; - } - } + if (rsq <= cutdistsq) { + jh = j; + if (history && rsq < radsum*radsum) + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; + } + } } } diff --git a/src/OPENMP/npair_half_size_multi_newton_omp.cpp b/src/OPENMP/npair_half_size_multi_newton_omp.cpp index cdc68e1b42..df3d7abc9f 100644 --- a/src/OPENMP/npair_half_size_multi_newton_omp.cpp +++ b/src/OPENMP/npair_half_size_multi_newton_omp.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_multi_newton_omp.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neighbor.h" #include "neigh_list.h" @@ -40,8 +43,10 @@ NPairHalfSizeMultiNewtonOmp::NPairHalfSizeMultiNewtonOmp(LAMMPS *lmp) : NPair(lm void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -49,7 +54,9 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns; + int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns; + int which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -62,7 +69,14 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -83,6 +97,11 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } ibin = atom2bin[i]; @@ -107,7 +126,7 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) // if j is owned atom, store it if j > i // if j is ghost, only store if j coords are "above and to the right" of i - for (j = js; j >= 0; j = bins[j]) { + for (j = js; j >= 0; j = bins[j]) { if(icollection != jcollection and j < i) continue; if (j >= nlocal) { @@ -121,19 +140,32 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; - } + if (rsq <= cutdistsq) { + jh = j; + if (history && rsq < radsum*radsum) + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; + } } } @@ -142,31 +174,43 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) // stencil is half if i same size as j // stencil is full if i smaller than j - s = stencil_multi[icollection][jcollection]; - ns = nstencil_multi[icollection][jcollection]; + s = stencil_multi[icollection][jcollection]; + ns = nstencil_multi[icollection][jcollection]; - for (k = 0; k < ns; k++) { - js = binhead_multi[jcollection][jbin + s[k]]; - for (j = js; j >= 0; j = bins[j]) { + for (k = 0; k < ns; k++) { + js = binhead_multi[jcollection][jbin + s[k]]; + for (j = js; j >= 0; j = bins[j]) { jtype = type[j]; - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; - } - } + if (rsq <= cutdistsq) { + if (history && rsq < radsum*radsum) + j = j ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } + } + } } ilist[i] = i; diff --git a/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp b/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp index d38e52fe18..717e023f4a 100644 --- a/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_size_multi_newton_tri_omp.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_multi_newton_tri_omp.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neighbor.h" #include "neigh_list.h" @@ -41,8 +44,10 @@ NPairHalfSizeMultiNewtonTriOmp::NPairHalfSizeMultiNewtonTriOmp(LAMMPS *lmp) : void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -50,7 +55,9 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns; + int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns; + int which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -63,7 +70,14 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -84,6 +98,11 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } ibin = atom2bin[i]; @@ -104,12 +123,12 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) // (equal zyx and j <= i) // latter excludes self-self interaction but allows superposed atoms - s = stencil_multi[icollection][jcollection]; - ns = nstencil_multi[icollection][jcollection]; + s = stencil_multi[icollection][jcollection]; + ns = nstencil_multi[icollection][jcollection]; - for (k = 0; k < ns; k++) { - js = binhead_multi[jcollection][jbin + s[k]]; - for (j = js; j >= 0; j = bins[j]) { + for (k = 0; k < ns; k++) { + js = binhead_multi[jcollection][jbin + s[k]]; + for (j = js; j >= 0; j = bins[j]) { // if same size (same collection), use half stencil if(cutcollectionsq[icollection][icollection] == cutcollectionsq[jcollection][jcollection]){ @@ -126,21 +145,34 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; - } - } + if (rsq <= cutdistsq) { + jh = j; + if (history && rsq < radsum*radsum) + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } + } + } } ilist[i] = i; diff --git a/src/OPENMP/npair_half_size_multi_old_newtoff_omp.cpp b/src/OPENMP/npair_half_size_multi_old_newtoff_omp.cpp index eadcde2154..39210f1231 100644 --- a/src/OPENMP/npair_half_size_multi_old_newtoff_omp.cpp +++ b/src/OPENMP/npair_half_size_multi_old_newtoff_omp.cpp @@ -18,6 +18,8 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" +#include "domain.h" +#include "molecule.h" #include "my_page.h" #include "error.h" @@ -40,8 +42,10 @@ NPairHalfSizeMultiOldNewtoffOmp::NPairHalfSizeMultiOldNewtoffOmp(LAMMPS *lmp) : void NPairHalfSizeMultiOldNewtoffOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -49,7 +53,8 @@ void NPairHalfSizeMultiOldNewtoffOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,ns; + int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -59,7 +64,14 @@ void NPairHalfSizeMultiOldNewtoffOmp::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -79,6 +91,11 @@ void NPairHalfSizeMultiOldNewtoffOmp::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in other bins in stencil including self // only store pair if i < j @@ -107,10 +124,23 @@ void NPairHalfSizeMultiOldNewtoffOmp::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/OPENMP/npair_half_size_multi_old_newton_omp.cpp b/src/OPENMP/npair_half_size_multi_old_newton_omp.cpp index eb154dce3a..11a2441990 100644 --- a/src/OPENMP/npair_half_size_multi_old_newton_omp.cpp +++ b/src/OPENMP/npair_half_size_multi_old_newton_omp.cpp @@ -18,6 +18,8 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" +#include "domain.h" +#include "molecule.h" #include "my_page.h" #include "error.h" @@ -39,8 +41,10 @@ NPairHalfSizeMultiOldNewtonOmp::NPairHalfSizeMultiOldNewtonOmp(LAMMPS *lmp) : void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -48,7 +52,8 @@ void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,ns; + int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -58,7 +63,14 @@ void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -78,6 +90,11 @@ void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list @@ -102,10 +119,23 @@ void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } } @@ -132,10 +162,23 @@ void NPairHalfSizeMultiOldNewtonOmp::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp b/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp index e67c1463f7..03f9f6f8eb 100644 --- a/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp +++ b/src/OPENMP/npair_half_size_multi_old_newton_tri_omp.cpp @@ -18,6 +18,8 @@ #include "neigh_list.h" #include "atom.h" #include "atom_vec.h" +#include "domain.h" +#include "molecule.h" #include "my_page.h" #include "error.h" @@ -39,8 +41,10 @@ NPairHalfSizeMultiOldNewtonTriOmp::NPairHalfSizeMultiOldNewtonTriOmp(LAMMPS *lmp void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; #if defined(_OPENMP) @@ -48,7 +52,8 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,ns; + int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -58,7 +63,14 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -78,6 +90,11 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in bins, including self, in stencil // skip if i,j neighbor cutoff is less than bin distance @@ -115,10 +132,23 @@ void NPairHalfSizeMultiOldNewtonTriOmp::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/OPENMP/npair_half_size_nsq_newtoff_omp.cpp b/src/OPENMP/npair_half_size_nsq_newtoff_omp.cpp index f182492efe..7f4b22e729 100644 --- a/src/OPENMP/npair_half_size_nsq_newtoff_omp.cpp +++ b/src/OPENMP/npair_half_size_nsq_newtoff_omp.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_nsq_newtoff_omp.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "group.h" #include "my_page.h" #include "neigh_list.h" @@ -42,8 +45,10 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; @@ -52,7 +57,8 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,n; + int i,j,jh,n,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -61,9 +67,17 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + int nall = atom->nlocal + atom->nghost; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; @@ -81,6 +95,11 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over remaining atoms, owned and ghost @@ -96,10 +115,23 @@ void NPairHalfSizeNsqNewtoffOmp::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } } diff --git a/src/OPENMP/npair_half_size_nsq_newton_omp.cpp b/src/OPENMP/npair_half_size_nsq_newton_omp.cpp index fdea8f4468..9262b869af 100644 --- a/src/OPENMP/npair_half_size_nsq_newton_omp.cpp +++ b/src/OPENMP/npair_half_size_nsq_newton_omp.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_nsq_newton_omp.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "group.h" #include "my_page.h" #include "neigh_list.h" @@ -42,9 +45,11 @@ NPairHalfSizeNsqNewtonOmp::NPairHalfSizeNsqNewtonOmp(LAMMPS *lmp) : void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; - const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0;; + const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + const int molecular = atom->molecular; + const int moltemplate = (molecular == Atom::TEMPLATE) ? 1 : 0; const int history = list->history; - const int mask_history = 3 << SBBITS; + const int mask_history = 1 << HISTBITS; NPAIR_OMP_INIT; @@ -53,7 +58,8 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,n,itag,jtag; + int i,j,jh,n,itag,jtag,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -64,6 +70,13 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) int *type = atom->type; int *mask = atom->mask; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; @@ -84,6 +97,11 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over remaining atoms, owned and ghost @@ -115,10 +133,23 @@ void NPairHalfSizeNsqNewtonOmp::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + } else neighptr[n++] = jh; } } diff --git a/src/atom.cpp b/src/atom.cpp index 4135298673..ed7e95fe64 100644 --- a/src/atom.cpp +++ b/src/atom.cpp @@ -122,6 +122,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) omega = angmom = torque = nullptr; radius = rmass = nullptr; ellipsoid = line = tri = body = nullptr; + quat = nullptr; // molecular systems @@ -419,6 +420,8 @@ void Atom::peratom_create() add_peratom("tri",&tri,INT,0); add_peratom("body",&body,INT,0); + add_peratom("quat",&quat,DOUBLE,4); + // MOLECULE package add_peratom("molecule",&molecule,tagintsize,0); @@ -639,6 +642,7 @@ void Atom::set_atomflag_defaults() // identical list as 2nd customization in atom.h sphere_flag = ellipsoid_flag = line_flag = tri_flag = body_flag = 0; + quat_flag = 0; peri_flag = electron_flag = 0; wavepacket_flag = sph_flag = 0; molecule_flag = molindex_flag = molatom_flag = 0; @@ -2660,6 +2664,10 @@ length of the data area, and a short description. - int - 1 - 1 if the particle is a body particle, 0 if not + * - quat + - double + - 4 + - four quaternion components of the particles * - i_name - int - 1 @@ -2715,6 +2723,7 @@ void *Atom::extract(const char *name) if (strcmp(name,"line") == 0) return (void *) line; if (strcmp(name,"tri") == 0) return (void *) tri; if (strcmp(name,"body") == 0) return (void *) body; + if (strcmp(name,"quat") == 0) return (void *) quat; if (strcmp(name,"vfrac") == 0) return (void *) vfrac; if (strcmp(name,"s0") == 0) return (void *) s0; @@ -2837,6 +2846,7 @@ int Atom::extract_datatype(const char *name) if (strcmp(name,"line") == 0) return LAMMPS_INT; if (strcmp(name,"tri") == 0) return LAMMPS_INT; if (strcmp(name,"body") == 0) return LAMMPS_INT; + if (strcmp(name,"quat") == 0) return LAMMPS_DOUBLE_2D; if (strcmp(name,"vfrac") == 0) return LAMMPS_DOUBLE; if (strcmp(name,"s0") == 0) return LAMMPS_DOUBLE; diff --git a/src/atom.h b/src/atom.h index a16f1e2752..7e80151ef4 100644 --- a/src/atom.h +++ b/src/atom.h @@ -79,6 +79,7 @@ class Atom : protected Pointers { double *radius; double **omega, **angmom, **torque; int *ellipsoid, *line, *tri, *body; + double **quat; // molecular systems @@ -180,7 +181,7 @@ class Atom : protected Pointers { int molecule_flag, molindex_flag, molatom_flag; int q_flag, mu_flag; - int rmass_flag, radius_flag, omega_flag, torque_flag, angmom_flag; + int rmass_flag, radius_flag, omega_flag, torque_flag, angmom_flag, quat_flag; int vfrac_flag, spin_flag, eradius_flag, ervel_flag, erforce_flag; int cs_flag, csforce_flag, vforce_flag, ervelforce_flag, etag_flag; int rho_flag, esph_flag, cv_flag, vest_flag; diff --git a/src/atom_vec_sphere.cpp b/src/atom_vec_sphere.cpp index 5a69944501..02db644a2b 100644 --- a/src/atom_vec_sphere.cpp +++ b/src/atom_vec_sphere.cpp @@ -6,7 +6,7 @@ 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 distributead under + 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. diff --git a/src/bond.cpp b/src/bond.cpp index cc13ad6b86..e07e184cc9 100644 --- a/src/bond.cpp +++ b/src/bond.cpp @@ -42,8 +42,11 @@ Bond::Bond(LAMMPS *lmp) : Pointers(lmp) virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0; writedata = 1; + comm_forward = comm_reverse = comm_reverse_off = 0; + allocated = 0; suffix_flag = Suffix::NONE; + partial_flag = 0; maxeatom = maxvatom = 0; eatom = nullptr; diff --git a/src/bond.h b/src/bond.h index b2b2008b9e..ea213da476 100644 --- a/src/bond.h +++ b/src/bond.h @@ -25,11 +25,16 @@ class Bond : protected Pointers { public: int allocated; int *setflag; + int partial_flag; // 1 if bond type can be set to 0 and deleted int writedata; // 1 if writes coeffs to data file double energy; // accumulated energies double virial[6]; // accumulated virial: xx,yy,zz,xy,xz,yz double *eatom, **vatom; // accumulated per-atom energy/virial + int comm_forward; // size of forward communication (0 if none) + int comm_reverse; // size of reverse communication (0 if none) + int comm_reverse_off; // size of reverse comm even if newton off + int reinitflag; // 1 if compatible with fix adapt and alike // KOKKOS host/device flag and data masks @@ -55,6 +60,10 @@ class Bond : protected Pointers { virtual double memory_usage(); virtual void *extract(const char *, int &) { return nullptr; } virtual void reinit(); + virtual int pack_forward_comm(int, int *, double *, int, int *) {return 0;} + virtual void unpack_forward_comm(int, int, double *) {} + virtual int pack_reverse_comm(int, int, double *) {return 0;} + virtual void unpack_reverse_comm(int, int *, double *) {} void write_file(int, char **); diff --git a/src/comm.cpp b/src/comm.cpp index 2dddba11c2..16127b37e3 100644 --- a/src/comm.cpp +++ b/src/comm.cpp @@ -217,6 +217,9 @@ void Comm::init() if (force->pair) maxforward = MAX(maxforward,force->pair->comm_forward); if (force->pair) maxreverse = MAX(maxreverse,force->pair->comm_reverse); + if (force->bond) maxforward = MAX(maxforward,force->bond->comm_forward); + if (force->bond) maxreverse = MAX(maxreverse,force->bond->comm_reverse); + for (int i = 0; i < modify->nfix; i++) { maxforward = MAX(maxforward,modify->fix[i]->comm_forward); maxreverse = MAX(maxreverse,modify->fix[i]->comm_reverse); @@ -234,6 +237,7 @@ void Comm::init() if (force->newton == 0) maxreverse = 0; if (force->pair) maxreverse = MAX(maxreverse,force->pair->comm_reverse_off); + if (force->bond) maxreverse = MAX(maxreverse,force->bond->comm_reverse_off); // maxexchange_atom = size of an exchanged atom, set by AtomVec // only needs to be set if size > BUFEXTRA diff --git a/src/comm.h b/src/comm.h index bc5faa49f4..14253603b5 100644 --- a/src/comm.h +++ b/src/comm.h @@ -84,6 +84,8 @@ class Comm : protected Pointers { virtual void forward_comm_pair(class Pair *) = 0; virtual void reverse_comm_pair(class Pair *) = 0; + virtual void forward_comm_bond(class Bond *) = 0; + virtual void reverse_comm_bond(class Bond *) = 0; virtual void forward_comm_fix(class Fix *, int size = 0) = 0; virtual void reverse_comm_fix(class Fix *, int size = 0) = 0; virtual void reverse_comm_fix_variable(class Fix *) = 0; diff --git a/src/comm_brick.cpp b/src/comm_brick.cpp index ab1b27f1b2..8dacdfb29b 100644 --- a/src/comm_brick.cpp +++ b/src/comm_brick.cpp @@ -20,6 +20,7 @@ #include "atom.h" #include "atom_vec.h" +#include "bond.h" #include "compute.h" #include "domain.h" #include "dump.h" @@ -1074,6 +1075,79 @@ void CommBrick::reverse_comm_pair(Pair *pair) } } +/* ---------------------------------------------------------------------- + forward communication invoked by a Bond + nsize used only to set recv buffer limit +------------------------------------------------------------------------- */ + +void CommBrick::forward_comm_bond(Bond *bond) +{ + int iswap,n; + double *buf; + MPI_Request request; + + int nsize = bond->comm_forward; + + for (iswap = 0; iswap < nswap; iswap++) { + + // pack buffer + + n = bond->pack_forward_comm(sendnum[iswap],sendlist[iswap], + buf_send,pbc_flag[iswap],pbc[iswap]); + + // exchange with another proc + // if self, set recv buffer to send buffer + + if (sendproc[iswap] != me) { + if (recvnum[iswap]) + MPI_Irecv(buf_recv,nsize*recvnum[iswap],MPI_DOUBLE, + recvproc[iswap],0,world,&request); + if (sendnum[iswap]) + MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap],0,world); + if (recvnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE); + buf = buf_recv; + } else buf = buf_send; + + // unpack buffer + + bond->unpack_forward_comm(recvnum[iswap],firstrecv[iswap],buf); + } +} + +/* ---------------------------------------------------------------------- + reverse communication invoked by a Bond + nsize used only to set recv buffer limit +------------------------------------------------------------------------- */ + +void CommBrick::reverse_comm_bond(Bond *bond) +{ + int iswap,n; + double *buf; + MPI_Request request; + + int nsize = MAX(bond->comm_reverse,bond->comm_reverse_off); + + for (iswap = nswap-1; iswap >= 0; iswap--) { + + // pack buffer + + n = bond->pack_reverse_comm(recvnum[iswap],firstrecv[iswap],buf_send); + + // exchange with another proc + // if self, set recv buffer to send buffer + + if (sendproc[iswap] != me) { + if (sendnum[iswap]) + MPI_Irecv(buf_recv,nsize*sendnum[iswap],MPI_DOUBLE,sendproc[iswap],0, + world,&request); + if (recvnum[iswap]) + MPI_Send(buf_send,n,MPI_DOUBLE,recvproc[iswap],0,world); + if (sendnum[iswap]) MPI_Wait(&request,MPI_STATUS_IGNORE); + buf = buf_recv; + } else buf = buf_send; + } +} + /* ---------------------------------------------------------------------- forward communication invoked by a Fix size/nsize used only to set recv buffer limit diff --git a/src/comm_brick.h b/src/comm_brick.h index 9ff177837d..48863864d9 100644 --- a/src/comm_brick.h +++ b/src/comm_brick.h @@ -33,6 +33,8 @@ class CommBrick : public Comm { virtual void forward_comm_pair(class Pair *); // forward comm from a Pair virtual void reverse_comm_pair(class Pair *); // reverse comm from a Pair + virtual void forward_comm_bond(class Bond *); // forward comm from a Bond + virtual void reverse_comm_bond(class Bond *); // reverse comm from a Bond virtual void forward_comm_fix(class Fix *, int size = 0); // forward comm from a Fix virtual void reverse_comm_fix(class Fix *, int size = 0); diff --git a/src/comm_tiled.cpp b/src/comm_tiled.cpp index e9175094f8..9267ac3c69 100644 --- a/src/comm_tiled.cpp +++ b/src/comm_tiled.cpp @@ -21,6 +21,7 @@ #include "atom.h" #include "atom_vec.h" +#include "bond.h" #include "compute.h" #include "domain.h" #include "dump.h" @@ -1464,6 +1465,99 @@ void CommTiled::reverse_comm_pair(Pair *pair) } } +/* ---------------------------------------------------------------------- + forward communication invoked by a Bond + nsize used only to set recv buffer limit +------------------------------------------------------------------------- */ + +void CommTiled::forward_comm_bond(Bond *bond) +{ + int i,irecv,n,nsend,nrecv; + + int nsize = bond->comm_forward; + + for (int iswap = 0; iswap < nswap; iswap++) { + nsend = nsendproc[iswap] - sendself[iswap]; + nrecv = nrecvproc[iswap] - sendself[iswap]; + + if (recvother[iswap]) { + for (i = 0; i < nrecv; i++) + MPI_Irecv(&buf_recv[nsize*forward_recv_offset[iswap][i]], + nsize*recvnum[iswap][i], + MPI_DOUBLE,recvproc[iswap][i],0,world,&requests[i]); + } + + if (sendother[iswap]) { + for (i = 0; i < nsend; i++) { + n = bond->pack_forward_comm(sendnum[iswap][i],sendlist[iswap][i], + buf_send,pbc_flag[iswap][i],pbc[iswap][i]); + MPI_Send(buf_send,n,MPI_DOUBLE,sendproc[iswap][i],0,world); + } + } + + if (sendself[iswap]) { + bond->pack_forward_comm(sendnum[iswap][nsend],sendlist[iswap][nsend], + buf_send,pbc_flag[iswap][nsend], + pbc[iswap][nsend]); + bond->unpack_forward_comm(recvnum[iswap][nrecv],firstrecv[iswap][nrecv], + buf_send); + } + if (recvother[iswap]) { + for (i = 0; i < nrecv; i++) { + MPI_Waitany(nrecv,requests,&irecv,MPI_STATUS_IGNORE); + bond->unpack_forward_comm(recvnum[iswap][irecv],firstrecv[iswap][irecv], + &buf_recv[nsize* + forward_recv_offset[iswap][irecv]]); + } + } + } +} + +/* ---------------------------------------------------------------------- + reverse communication invoked by a Bond + nsize used only to set recv buffer limit +------------------------------------------------------------------------- */ + +void CommTiled::reverse_comm_bond(Bond *bond) +{ + int i,irecv,n,nsend,nrecv; + + int nsize = MAX(bond->comm_reverse,bond->comm_reverse_off); + + for (int iswap = nswap-1; iswap >= 0; iswap--) { + nsend = nsendproc[iswap] - sendself[iswap]; + nrecv = nrecvproc[iswap] - sendself[iswap]; + + if (sendother[iswap]) { + for (i = 0; i < nsend; i++) + MPI_Irecv(&buf_recv[nsize*reverse_recv_offset[iswap][i]], + nsize*sendnum[iswap][i],MPI_DOUBLE, + sendproc[iswap][i],0,world,&requests[i]); + } + if (recvother[iswap]) { + for (i = 0; i < nrecv; i++) { + n = bond->pack_reverse_comm(recvnum[iswap][i],firstrecv[iswap][i], + buf_send); + MPI_Send(buf_send,n,MPI_DOUBLE,recvproc[iswap][i],0,world); + } + } + if (sendself[iswap]) { + bond->pack_reverse_comm(recvnum[iswap][nrecv],firstrecv[iswap][nrecv], + buf_send); + bond->unpack_reverse_comm(sendnum[iswap][nsend],sendlist[iswap][nsend], + buf_send); + } + if (sendother[iswap]) { + for (i = 0; i < nsend; i++) { + MPI_Waitany(nsend,requests,&irecv,MPI_STATUS_IGNORE); + bond->unpack_reverse_comm(sendnum[iswap][irecv],sendlist[iswap][irecv], + &buf_recv[nsize* + reverse_recv_offset[iswap][irecv]]); + } + } + } +} + /* ---------------------------------------------------------------------- forward communication invoked by a Fix size/nsize used only to set recv buffer limit diff --git a/src/comm_tiled.h b/src/comm_tiled.h index 75109e097d..ec43b50787 100644 --- a/src/comm_tiled.h +++ b/src/comm_tiled.h @@ -33,6 +33,8 @@ class CommTiled : public Comm { virtual void forward_comm_pair(class Pair *); // forward comm from a Pair virtual void reverse_comm_pair(class Pair *); // reverse comm from a Pair + virtual void forward_comm_bond(class Bond *); // forward comm from a Bond + virtual void reverse_comm_bond(class Bond *); // reverse comm from a Bond virtual void forward_comm_fix(class Fix *, int size = 0); // forward comm from a Fix virtual void reverse_comm_fix(class Fix *, int size = 0); diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp index d4cec70fe4..e7c9b9a40c 100644 --- a/src/compute_property_atom.cpp +++ b/src/compute_property_atom.cpp @@ -231,25 +231,25 @@ ComputePropertyAtom::ComputePropertyAtom(LAMMPS *lmp, int narg, char **arg) : } else if (strcmp(arg[iarg],"quatw") == 0) { avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); avec_body = (AtomVecBody *) atom->style_match("body"); - if (!avec_ellipsoid && !avec_body) + if (!avec_ellipsoid && !avec_body && !atom->quat_flag) error->all(FLERR,"Compute property/atom for atom property that isn't allocated"); pack_choice[i] = &ComputePropertyAtom::pack_quatw; } else if (strcmp(arg[iarg],"quati") == 0) { avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); avec_body = (AtomVecBody *) atom->style_match("body"); - if (!avec_ellipsoid && !avec_body) + if (!avec_ellipsoid && !avec_body && !atom->quat_flag) error->all(FLERR,"Compute property/atom for atom property that isn't allocated"); pack_choice[i] = &ComputePropertyAtom::pack_quati; } else if (strcmp(arg[iarg],"quatj") == 0) { avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); avec_body = (AtomVecBody *) atom->style_match("body"); - if (!avec_ellipsoid && !avec_body) + if (!avec_ellipsoid && !avec_body && !atom->quat_flag) error->all(FLERR,"Compute property/atom for atom property that isn't allocated"); pack_choice[i] = &ComputePropertyAtom::pack_quatj; } else if (strcmp(arg[iarg],"quatk") == 0) { avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid"); avec_body = (AtomVecBody *) atom->style_match("body"); - if (!avec_ellipsoid && !avec_body) + if (!avec_ellipsoid && !avec_body && !atom->quat_flag) error->all(FLERR,"Compute property/atom for atom property that isn't allocated"); pack_choice[i] = &ComputePropertyAtom::pack_quatk; @@ -1334,7 +1334,7 @@ void ComputePropertyAtom::pack_quatw(int n) n += nvalues; } - } else { + } else if (avec_body) { AtomVecBody::Bonus *bonus = avec_body->bonus; int *body = atom->body; int *mask = atom->mask; @@ -1346,6 +1346,17 @@ void ComputePropertyAtom::pack_quatw(int n) else buf[n] = 0.0; n += nvalues; } + } else { + double **quat = atom->quat; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = quat[i][0]; + else buf[n] = 0.0; + n += nvalues; + } } } @@ -1366,7 +1377,7 @@ void ComputePropertyAtom::pack_quati(int n) n += nvalues; } - } else { + } else if (avec_body) { AtomVecBody::Bonus *bonus = avec_body->bonus; int *body = atom->body; int *mask = atom->mask; @@ -1378,6 +1389,17 @@ void ComputePropertyAtom::pack_quati(int n) else buf[n] = 0.0; n += nvalues; } + } else { + double **quat = atom->quat; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = quat[i][1]; + else buf[n] = 0.0; + n += nvalues; + } } } @@ -1398,7 +1420,7 @@ void ComputePropertyAtom::pack_quatj(int n) n += nvalues; } - } else { + } else if (avec_body) { AtomVecBody::Bonus *bonus = avec_body->bonus; int *body = atom->body; int *mask = atom->mask; @@ -1410,6 +1432,17 @@ void ComputePropertyAtom::pack_quatj(int n) else buf[n] = 0.0; n += nvalues; } + } else { + double **quat = atom->quat; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = quat[i][2]; + else buf[n] = 0.0; + n += nvalues; + } } } @@ -1430,7 +1463,7 @@ void ComputePropertyAtom::pack_quatk(int n) n += nvalues; } - } else { + } else if (avec_body) { AtomVecBody::Bonus *bonus = avec_body->bonus; int *body = atom->body; int *mask = atom->mask; @@ -1442,6 +1475,17 @@ void ComputePropertyAtom::pack_quatk(int n) else buf[n] = 0.0; n += nvalues; } + } else { + double **quat = atom->quat; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) + buf[n] = quat[i][3]; + else buf[n] = 0.0; + n += nvalues; + } } } diff --git a/src/fix_neigh_history.cpp b/src/fix_neigh_history.cpp index 4ff7e7841b..2e5c15d7ba 100644 --- a/src/fix_neigh_history.cpp +++ b/src/fix_neigh_history.cpp @@ -627,7 +627,7 @@ void FixNeighHistory::post_neighbor() j = jlist[jj]; if (use_bit_flag) { - rflag = sbmask(j) | pair->beyond_contact; + rflag = histmask(j) | pair->beyond_contact; j &= NEIGHMASK; jlist[jj] = j; } else { diff --git a/src/fix_neigh_history.h b/src/fix_neigh_history.h index 1405c05731..f382207541 100644 --- a/src/fix_neigh_history.h +++ b/src/fix_neigh_history.h @@ -93,7 +93,8 @@ class FixNeighHistory : public Fix { virtual void pre_exchange_no_newton(); void allocate_pages(); - inline int sbmask(int j) const { return j >> SBBITS & 3; } + // Shift by HISTBITS and check the first bit + inline int histmask(int j) const { return j >> HISTBITS & 1; } }; } // namespace LAMMPS_NS diff --git a/src/fix_property_atom.cpp b/src/fix_property_atom.cpp index 1e583a4ec9..558e44de48 100644 --- a/src/fix_property_atom.cpp +++ b/src/fix_property_atom.cpp @@ -550,18 +550,24 @@ void FixPropertyAtom::copy_arrays(int i, int j, int /*delflag*/) atom->q[j] = atom->q[i]; else if (style[nv] == RMASS) atom->rmass[j] = atom->rmass[i]; - else if (style[nv] == IVEC) + else if (style[nv] == IVEC) { atom->ivector[index[nv]][j] = atom->ivector[index[nv]][i]; - else if (style[nv] == DVEC) + atom->ivector[index[nv]][i] = 0; + } else if (style[nv] == DVEC) { atom->dvector[index[nv]][j] = atom->dvector[index[nv]][i]; - else if (style[nv] == IARRAY) { + atom->dvector[index[nv]][i] = 0.0; + } else if (style[nv] == IARRAY) { ncol = cols[nv]; - for (k = 0; k < ncol; k++) - atom->iarray[index[nv]][j][k] = atom->iarray[index[nv]][i][k]; + for (k = 0; k < ncol; k++) { + atom->iarray[index[nv]][j][k] = atom->iarray[index[nv]][i][k]; + atom->iarray[index[nv]][i][k] = 0; + } } else if (style[nv] == DARRAY) { ncol = cols[nv]; - for (k = 0; k < ncol; k++) - atom->darray[index[nv]][j][k] = atom->darray[index[nv]][i][k]; + for (k = 0; k < ncol; k++) { + atom->darray[index[nv]][j][k] = atom->darray[index[nv]][i][k]; + atom->darray[index[nv]][i][k] = 0.0; + } } } } @@ -696,16 +702,24 @@ int FixPropertyAtom::pack_exchange(int i, double *buf) if (style[nv] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d; else if (style[nv] == CHARGE) buf[m++] = atom->q[i]; else if (style[nv] == RMASS) buf[m++] = atom->rmass[i]; - else if (style[nv] == IVEC) buf[m++] = ubuf(atom->ivector[index[nv]][i]).d; - else if (style[nv] == DVEC) buf[m++] = atom->dvector[index[nv]][i]; - else if (style[nv] == IARRAY) { + else if (style[nv] == IVEC) { + buf[m++] = ubuf(atom->ivector[index[nv]][i]).d; + atom->ivector[index[nv]][i] = 0; + } else if (style[nv] == DVEC) { + buf[m++] = atom->dvector[index[nv]][i]; + atom->dvector[index[nv]][i] = 0; + } else if (style[nv] == IARRAY) { ncol = cols[nv]; - for (k = 0; k < ncol; k++) - buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d; + for (k = 0; k < ncol; k++) { + buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d; + atom->iarray[index[nv]][i][k] = 0; + } } else if (style[nv] == DARRAY) { ncol = cols[nv]; - for (k = 0; k < ncol; k++) - buf[m++] = atom->darray[index[nv]][i][k]; + for (k = 0; k < ncol; k++) { + buf[m++] = atom->darray[index[nv]][i][k]; + atom->darray[index[nv]][i][k] = 0.0; + } } } diff --git a/src/fix_store_local.cpp b/src/fix_store_local.cpp new file mode 100644 index 0000000000..7067d83ed1 --- /dev/null +++ b/src/fix_store_local.cpp @@ -0,0 +1,200 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ 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. +------------------------------------------------------------------------- */ + +#include "fix_store_local.h" +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "memory.h" +#include "update.h" + +using namespace LAMMPS_NS; +using namespace FixConst; + +#define DELTA 1000 + +/* ---------------------------------------------------------------------- */ + +FixStoreLocal::FixStoreLocal(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), nvalues(0), vector(nullptr), array(nullptr) +{ + if (narg != 5) error->all(FLERR, "Illegal fix store/local command"); + local_flag = 1; + + nevery = utils::inumeric(FLERR, arg[3], false, lmp); + if (nevery <= 0) error->all(FLERR, "Illegal fix store/local command"); + local_freq = nevery; + + nvalues = utils::inumeric(FLERR, arg[4], false, lmp); + + if (nvalues <= 0) error->all(FLERR, "Illegal fix store/local command"); + if (nvalues == 1) + size_local_cols = 0; + else + size_local_cols = nvalues; + size_local_rows = 0; + + vector = nullptr; + array = nullptr; + nmax = 0; + ncount = 0; +} + +/* ---------------------------------------------------------------------- */ + +FixStoreLocal::~FixStoreLocal() +{ + memory->destroy(vector); + memory->destroy(array); +} + +/* ---------------------------------------------------------------------- */ + +int FixStoreLocal::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixStoreLocal::add_data(double *input_data, int i, int j) +{ + int *mask = atom->mask; + if (!(mask[i] & groupbit)) return; + if (!(mask[j] & groupbit)) return; + + if (ncount >= nmax) reallocate(ncount); + + // fill vector or array with local values + if (nvalues == 1) { + vector[ncount] = input_data[0]; + } else { + for (int i = 0; i < nvalues; i++) array[ncount][i] = input_data[i]; + } + + ncount += 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixStoreLocal::post_force(int /*vflag*/) +{ + if (update->ntimestep % nevery == 0) { + size_local_rows = ncount; + ncount = 0; + } +} + +/* ---------------------------------------------------------------------- */ + +void FixStoreLocal::reallocate(int n) +{ + // grow vector or array + while (nmax <= n) nmax += DELTA; + + if (nvalues == 1) { + memory->grow(vector, nmax, "fix_store_local:vector"); + vector_local = vector; + } else { + memory->grow(array, nmax, nvalues, "fix_store_local:array"); + array_local = array; + } +} + +/* ---------------------------------------------------------------------- + write global array to restart file +------------------------------------------------------------------------- */ + +void FixStoreLocal::write_restart(FILE *fp) +{ + // fill rbuf with size and vec/array values + + rbuf[0] = nmax; + rbuf[1] = nvalues; + if (nvalues == 1) memcpy(&rbuf[2],vector,sizeof(double)*nmax); + else memcpy(&rbuf[2],&array[0][0],sizeof(double)*nmax*nvalues); + + int n = nmax*nvalues + 2; + if (comm->me == 0) { + int size = n * sizeof(double); + fwrite(&size,sizeof(int),1,fp); + fwrite(rbuf,sizeof(double),n,fp); + } +} + +/* ---------------------------------------------------------------------- + use global array from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixStoreLocal::restart(char *buf) +{ + // first 2 values in buf are vec/array sizes + + double *dbuf = (double *) buf; + int nrow_restart = dbuf[0]; + int ncol_restart = dbuf[1]; + + // if size of vec/array has changed, + // means the restart file is setting size of vec or array and doing init + // because caller did not know size at time this fix was instantiated + // reallocate vector or array accordingly + + if (nmax != nrow_restart || nvalues != ncol_restart) { + memory->destroy(vector); + memory->destroy(array); + memory->destroy(rbuf); + vector = nullptr; + array = nullptr; + + nmax = nrow_restart; + nvalues = ncol_restart; + if (nvalues == 1) memory->create(vector,nmax,"fix/store/local:vector"); + else memory->create(array,nmax,nvalues,"fix/store/local:array"); + memory->create(rbuf,nmax*nvalues+2,"fix/store:rbuf"); + } + + int n = nmax*nvalues; + if (nvalues == 1) memcpy(vector,&dbuf[2],n*sizeof(double)); + else memcpy(&array[0][0],&dbuf[2],n*sizeof(double)); +} + +/* ---------------------------------------------------------------------- + maxsize of any atom's restart data +------------------------------------------------------------------------- */ + +int FixStoreLocal::maxsize_restart() +{ + return nvalues+1; +} + +/* ---------------------------------------------------------------------- + size of atom nlocal's restart data +------------------------------------------------------------------------- */ + +int FixStoreLocal::size_restart(int /*nlocal*/) +{ + return nvalues+1; +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double FixStoreLocal::memory_usage() +{ + double bytes = (double) nmax * (double) nvalues * sizeof(double); + return bytes; +} + diff --git a/src/MISC/fix_pair_tracker.h b/src/fix_store_local.h similarity index 60% rename from src/MISC/fix_pair_tracker.h rename to src/fix_store_local.h index 7c2e3ff322..4102c88f37 100644 --- a/src/MISC/fix_pair_tracker.h +++ b/src/fix_store_local.h @@ -1,85 +1,75 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - https://www.lammps.org/ 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 -// clang-format off -FixStyle(pair/tracker,FixPairTracker); -// clang-format on -#else - -#ifndef LMP_FIX_PAIR_TRACKING_H -#define LMP_FIX_PAIR_TRACKING_H - -#include "fix.h" - -namespace LAMMPS_NS { - -class FixPairTracker : public Fix { - public: - FixPairTracker(class LAMMPS *, int, char **); - ~FixPairTracker(); - int setmask(); - void init(); - void post_force(int); - double memory_usage(); - void lost_contact(int, int, double, double, double, double); - - private: - int nvalues, nmax; - int index_i, index_j; - double tmin, rmin, rsum, time_initial, nstep_initial; - - double *vector; - double **array; - int **type_filter; - - int ncount; - - void reallocate(int); - - typedef void (FixPairTracker::*FnPtrPack)(int); - FnPtrPack *pack_choice; // ptrs to pack functions - - void pack_id1(int); - void pack_id2(int); - - void pack_time_created(int); - void pack_time_broken(int); - void pack_time_total(int); - - void pack_x(int); - void pack_y(int); - void pack_z(int); - - void pack_rmin(int); - void pack_rave(int); -}; - -} // namespace LAMMPS_NS - -#endif -#endif - -/* ERROR/WARNING messages: - -E: Illegal ... command - -Self-explanatory. Check the input script syntax and compare to the -documentation for the command. You can use -echo screen as a -command-line option when running LAMMPS to see the offending line. - -E: Invalid keyword in fix pair/tracker command - -Self-explanatory. - -*/ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/ 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 +// clang-format off +FixStyle(store/local,FixStoreLocal); +// clang-format on +#else + +#ifndef LMP_FIX_STORE_LOCAL_H +#define LMP_FIX_STORE_LOCAL_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixStoreLocal : public Fix { + public: + FixStoreLocal(class LAMMPS *, int, char **); + ~FixStoreLocal(); + int setmask(); + void post_force(int); + void write_restart(FILE *); + void restart(char *); + int size_restart(int); + int maxsize_restart(); + double memory_usage(); + void add_data(double *, int, int); + int nvalues; + + private: + int nmax; + + double *vector; + double **array; + + int ncount; + + void reallocate(int); + double *rbuf; // restart buffer for GLOBAL vec/array +}; + +} // namespace LAMMPS_NS + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Invalid keyword in fix store/local command + +Self-explanatory. + +E: Unused instance of fix store/local + +Instance of fix store/local is not associated with any other LAMMPS +class such as a bond style, pair style, etc. + +*/ diff --git a/src/lmptype.h b/src/lmptype.h index 871bf5ff6c..b188e4b7e7 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -56,10 +56,12 @@ namespace LAMMPS_NS { // reserve 2 hi bits in molecular system neigh list for special bonds flag -// max local + ghost atoms per processor = 2^30 - 1 +// reserve 3rd last bit in neigh list for fix neigh/history flag +// max local + ghost atoms per processor = 2^29 - 1 #define SBBITS 30 -#define NEIGHMASK 0x3FFFFFFF +#define HISTBITS 29 +#define NEIGHMASK 0x1FFFFFFF // default to 32-bit smallint and other ints, 64-bit bigint diff --git a/src/math_extra.cpp b/src/math_extra.cpp index 54f6204182..15d6e472c9 100644 --- a/src/math_extra.cpp +++ b/src/math_extra.cpp @@ -143,6 +143,58 @@ void richardson(double *q, double *m, double *w, double *moments, double dtq) MathExtra::qnormalize(q); } +/* ---------------------------------------------------------------------- + Richardson iteration to update quaternion from angular velocity + return new normalized quaternion q + also returns updated omega at 1/2 step + Assumes spherical particles - no need to rotate to match moments +------------------------------------------------------------------------- */ + +void richardson_sphere(double *q, double *w, double dtq) +{ + // full update from dq/dt = 1/2 w q + + double wq[4]; + MathExtra::vecquat(w,q,wq); + + double qfull[4]; + qfull[0] = q[0] + dtq * wq[0]; + qfull[1] = q[1] + dtq * wq[1]; + qfull[2] = q[2] + dtq * wq[2]; + qfull[3] = q[3] + dtq * wq[3]; + MathExtra::qnormalize(qfull); + + // 1st half update from dq/dt = 1/2 w q + + double qhalf[4]; + qhalf[0] = q[0] + 0.5*dtq * wq[0]; + qhalf[1] = q[1] + 0.5*dtq * wq[1]; + qhalf[2] = q[2] + 0.5*dtq * wq[2]; + qhalf[3] = q[3] + 0.5*dtq * wq[3]; + MathExtra::qnormalize(qhalf); + + // re-compute q at 1/2 step + // recompute wq + + MathExtra::vecquat(w,qhalf,wq); + + // 2nd half update from dq/dt = 1/2 w q + + qhalf[0] += 0.5*dtq * wq[0]; + qhalf[1] += 0.5*dtq * wq[1]; + qhalf[2] += 0.5*dtq * wq[2]; + qhalf[3] += 0.5*dtq * wq[3]; + MathExtra::qnormalize(qhalf); + + // corrected Richardson update + + q[0] = 2.0*qhalf[0] - qfull[0]; + q[1] = 2.0*qhalf[1] - qfull[1]; + q[2] = 2.0*qhalf[2] - qfull[2]; + q[3] = 2.0*qhalf[3] - qfull[3]; + MathExtra::qnormalize(q); +} + /* ---------------------------------------------------------------------- apply evolution operators to quat, quat momentum Miller et al., J Chem Phys. 116, 8649-8659 (2002) diff --git a/src/math_extra.h b/src/math_extra.h index 7c71e4c11b..55d680b0f2 100644 --- a/src/math_extra.h +++ b/src/math_extra.h @@ -76,6 +76,7 @@ void write3(const double mat[3][3]); int mldivide3(const double mat[3][3], const double *vec, double *ans); void rotate(double matrix[3][3], int i, int j, int k, int l, double s, double tau); void richardson(double *q, double *m, double *w, double *moments, double dtq); +void richardson_sphere(double *q, double *w, double dtq); void no_squish_rotate(int k, double *p, double *q, double *inertia, double dt); // shape matrix operations @@ -91,6 +92,7 @@ inline void vecquat(double *a, double *b, double *c); inline void quatvec(double *a, double *b, double *c); inline void quatquat(double *a, double *b, double *c); inline void invquatvec(double *a, double *b, double *c); +inline void quatrotvec(double *a, double *b, double *c); inline void axisangle_to_quat(const double *v, const double angle, double *quat); void angmom_to_omega(double *m, double *ex, double *ey, double *ez, double *idiag, double *w); @@ -651,6 +653,29 @@ inline void MathExtra::invquatvec(double *a, double *b, double *c) c[2] = -a[3] * b[0] + a[2] * b[1] - a[1] * b[2] + a[0] * b[3]; } +/* ---------------------------------------------------------------------- + quaternion rotation of vector: c = a*b*conj(a) + a is a quaternion + b is a three component vector + c is a three component vector +------------------------------------------------------------------------- */ + +inline void MathExtra::quatrotvec(double *a, double *b, double *c) +{ + double temp[4]; + + // temp = a*b + temp[0] = -a[1]*b[0] - a[2]*b[1] - a[3]*b[2]; + temp[1] = a[0]*b[0] + a[2]*b[2] - a[3]*b[1]; + temp[2] = a[0]*b[1] + a[3]*b[0] - a[1]*b[2]; + temp[3] = a[0]*b[2] + a[1]*b[1] - a[2]*b[0]; + + // c = temp*conj(a) + c[0] = -a[1]*temp[0] + a[0]*temp[1] - a[3]*temp[2] + a[2]*temp[3]; + c[1] = -a[2]*temp[0] + a[3]*temp[1] + a[0]*temp[2] - a[1]*temp[3]; + c[2] = -a[3]*temp[0] - a[2]*temp[1] + a[1]*temp[2] + a[0]*temp[3]; +} + /* ---------------------------------------------------------------------- compute quaternion from axis-angle rotation v MUST be a unit vector diff --git a/src/neighbor.cpp b/src/neighbor.cpp index b7217926ca..d2fab58c24 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -21,6 +21,7 @@ #include "atom.h" #include "atom_vec.h" +#include "bond.h" #include "citeme.h" #include "comm.h" #include "compute.h" @@ -1444,7 +1445,7 @@ void Neighbor::init_topology() // bonds,etc can only be broken for atom->molecular = Atom::MOLECULAR, not Atom::TEMPLATE // SHAKE sets bonds and angles negative // gcmc sets all bonds, angles, etc negative - // bond_quartic sets bonds to 0 + // partial_flag sets bonds to 0 // delete_bonds sets all interactions negative int bond_off = 0; @@ -1453,7 +1454,9 @@ void Neighbor::init_topology() if (utils::strmatch(modify->fix[i]->style,"^shake") || utils::strmatch(modify->fix[i]->style,"^rattle")) bond_off = angle_off = 1; - if (force->bond && force->bond_match("quartic")) bond_off = 1; + if (force->bond) + if (force->bond->partial_flag) + bond_off = 1; if (atom->avec->bonds_allow && atom->molecular == Atom::MOLECULAR) { for (i = 0; i < atom->nlocal; i++) { diff --git a/src/npair_half_size_bin_newtoff.cpp b/src/npair_half_size_bin_newtoff.cpp index e836503190..001d1281c9 100644 --- a/src/npair_half_size_bin_newtoff.cpp +++ b/src/npair_half_size_bin_newtoff.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_bin_newtoff.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neigh_list.h" @@ -35,7 +38,8 @@ NPairHalfSizeBinNewtoff::NPairHalfSizeBinNewtoff(LAMMPS *lmp) : NPair(lmp) {} void NPairHalfSizeBinNewtoff::build(NeighList *list) { - int i,j,k,n,ibin; + int i,j,jh,k,n,ibin,which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -44,17 +48,26 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -68,6 +81,11 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list) ztmp = x[i][2]; radi = radius[i]; ibin = atom2bin[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in surrounding bins in stencil including self // only store pair if i < j @@ -87,10 +105,24 @@ void NPairHalfSizeBinNewtoff::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/npair_half_size_bin_newton.cpp b/src/npair_half_size_bin_newton.cpp index b6a786b0cf..b3b1f3216d 100644 --- a/src/npair_half_size_bin_newton.cpp +++ b/src/npair_half_size_bin_newton.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_bin_newton.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neigh_list.h" @@ -34,7 +37,8 @@ NPairHalfSizeBinNewton::NPairHalfSizeBinNewton(LAMMPS *lmp) : NPair(lmp) {} void NPairHalfSizeBinNewton::build(NeighList *list) { - int i,j,k,n,ibin; + int i,j,jh,k,n,ibin,which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -43,17 +47,26 @@ void NPairHalfSizeBinNewton::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -66,6 +79,11 @@ void NPairHalfSizeBinNewton::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list @@ -90,10 +108,24 @@ void NPairHalfSizeBinNewton::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; } } @@ -112,10 +144,24 @@ void NPairHalfSizeBinNewton::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/npair_half_size_bin_newton_tri.cpp b/src/npair_half_size_bin_newton_tri.cpp index a86feb6f40..2704c809ba 100644 --- a/src/npair_half_size_bin_newton_tri.cpp +++ b/src/npair_half_size_bin_newton_tri.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_bin_newton_tri.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neigh_list.h" @@ -35,7 +38,8 @@ NPairHalfSizeBinNewtonTri::NPairHalfSizeBinNewtonTri(LAMMPS *lmp) : void NPairHalfSizeBinNewtonTri::build(NeighList *list) { - int i,j,k,n,ibin; + int i,j,jh,k,n,ibin,which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -44,17 +48,26 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -67,6 +80,11 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in bins in stencil // pairs for atoms j "below" i are excluded @@ -96,10 +114,24 @@ void NPairHalfSizeBinNewtonTri::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/npair_half_size_multi_newtoff.cpp b/src/npair_half_size_multi_newtoff.cpp index cb06ddfb6f..ab57a5b76e 100644 --- a/src/npair_half_size_multi_newtoff.cpp +++ b/src/npair_half_size_multi_newtoff.cpp @@ -16,7 +16,10 @@ es certain rights in this software. This software is distributed under #include "npair_half_size_multi_newtoff.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neighbor.h" #include "neigh_list.h" @@ -38,7 +41,9 @@ NPairHalfSizeMultiNewtoff::NPairHalfSizeMultiNewtoff(LAMMPS *lmp) : NPair(lmp) { void NPairHalfSizeMultiNewtoff::build(NeighList *list) { - int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns; + int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns; + int which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -49,17 +54,26 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -73,6 +87,11 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } ibin = atom2bin[i]; @@ -93,27 +112,41 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list) ns = nstencil_multi[icollection][jcollection]; for (k = 0; k < ns; k++) { - js = binhead_multi[jcollection][jbin + s[k]]; - for (j = js; j >= 0; j = bins[j]) { - if (j <= i) continue; + js = binhead_multi[jcollection][jbin + s[k]]; + for (j = js; j >= 0; j = bins[j]) { + if (j <= i) continue; jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; - } - } + if (rsq <= cutdistsq) { + jh = j; + if (history && rsq < radsum*radsum) + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; + } + } } } diff --git a/src/npair_half_size_multi_newton.cpp b/src/npair_half_size_multi_newton.cpp index 8af60aa435..6a5491247a 100644 --- a/src/npair_half_size_multi_newton.cpp +++ b/src/npair_half_size_multi_newton.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_multi_newton.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neighbor.h" #include "neigh_list.h" @@ -36,7 +39,9 @@ NPairHalfSizeMultiNewton::NPairHalfSizeMultiNewton(LAMMPS *lmp) : NPair(lmp) {} void NPairHalfSizeMultiNewton::build(NeighList *list) { - int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js; + int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js; + int which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -46,17 +51,26 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -70,6 +84,11 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } ibin = atom2bin[i]; @@ -94,33 +113,47 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) // if j is owned atom, store it if j > i // if j is ghost, only store if j coords are "above and to the right" of i - for (j = js; j >= 0; j = bins[j]) { + for (j = js; j >= 0; j = bins[j]) { if ((icollection != jcollection) && (j < i)) continue; - if (j >= nlocal) { - if (x[j][2] < ztmp) continue; - if (x[j][2] == ztmp) { - if (x[j][1] < ytmp) continue; - if (x[j][1] == ytmp && x[j][0] < xtmp) continue; - } - } + if (j >= nlocal) { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; - } + if (rsq <= cutdistsq) { + jh = j; + if (history && rsq < radsum*radsum) + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; + } } } @@ -129,30 +162,44 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) // stencil is half if i same size as j // stencil is full if i smaller than j - s = stencil_multi[icollection][jcollection]; - ns = nstencil_multi[icollection][jcollection]; + s = stencil_multi[icollection][jcollection]; + ns = nstencil_multi[icollection][jcollection]; - for (k = 0; k < ns; k++) { - js = binhead_multi[jcollection][jbin + s[k]]; - for (j = js; j >= 0; j = bins[j]) { + for (k = 0; k < ns; k++) { + js = binhead_multi[jcollection][jbin + s[k]]; + for (j = js; j >= 0; j = bins[j]) { jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - radsum = radi + radius[j]; - cutdistsq = (radsum+skin) * (radsum+skin); + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx*delx + dely*dely + delz*delz; + radsum = radi + radius[j]; + cutdistsq = (radsum+skin) * (radsum+skin); - if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; - } - } + if (rsq <= cutdistsq) { + jh = j; + if (history && rsq < radsum*radsum) + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; + } + } } } diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp index 20d4d8b421..ecb0036e02 100644 --- a/src/npair_half_size_multi_newton_tri.cpp +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_multi_newton_tri.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "my_page.h" #include "neighbor.h" #include "neigh_list.h" @@ -36,7 +39,9 @@ NPairHalfSizeMultiNewtonTri::NPairHalfSizeMultiNewtonTri(LAMMPS *lmp) : NPair(lm void NPairHalfSizeMultiNewtonTri::build(NeighList *list) { - int i,j,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js; + int i,j,jh,k,n,itype,jtype,icollection,jcollection,ibin,jbin,ns,js; + int which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -46,17 +51,26 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -70,6 +84,11 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } ibin = atom2bin[i]; @@ -119,10 +138,24 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/npair_half_size_multi_old_newtoff.cpp b/src/npair_half_size_multi_old_newtoff.cpp index 26bde10b0e..77781b7266 100644 --- a/src/npair_half_size_multi_old_newtoff.cpp +++ b/src/npair_half_size_multi_old_newtoff.cpp @@ -15,8 +15,11 @@ #include "npair_half_size_multi_old_newtoff.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" #include "molecule.h" +#include "molecule.h" #include "my_page.h" #include "neigh_list.h" @@ -37,7 +40,8 @@ NPairHalfSizeMultiOldNewtoff::NPairHalfSizeMultiOldNewtoff(LAMMPS *lmp) : NPair( void NPairHalfSizeMultiOldNewtoff::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,ns; + int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -47,17 +51,26 @@ void NPairHalfSizeMultiOldNewtoff::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -71,6 +84,11 @@ void NPairHalfSizeMultiOldNewtoff::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in other bins in stencil including self // only store pair if i < j @@ -99,10 +117,24 @@ void NPairHalfSizeMultiOldNewtoff::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/npair_half_size_multi_old_newton.cpp b/src/npair_half_size_multi_old_newton.cpp index 7f63679bc1..413090448c 100644 --- a/src/npair_half_size_multi_old_newton.cpp +++ b/src/npair_half_size_multi_old_newton.cpp @@ -15,6 +15,8 @@ #include "npair_half_size_multi_old_newton.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" #include "molecule.h" #include "my_page.h" @@ -36,7 +38,8 @@ NPairHalfSizeMultiOldNewton::NPairHalfSizeMultiOldNewton(LAMMPS *lmp) : NPair(lm void NPairHalfSizeMultiOldNewton::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,ns; + int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -46,17 +49,26 @@ void NPairHalfSizeMultiOldNewton::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -70,6 +82,11 @@ void NPairHalfSizeMultiOldNewton::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list @@ -95,10 +112,24 @@ void NPairHalfSizeMultiOldNewton::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; } } @@ -125,10 +156,23 @@ void NPairHalfSizeMultiOldNewton::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { - if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else + if (history && rsq < radsum*radsum) + j = j ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = j; + else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; + else if (which > 0) neighptr[n++] = j ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } diff --git a/src/npair_half_size_multi_old_newton_tri.cpp b/src/npair_half_size_multi_old_newton_tri.cpp index f71a470479..9117e1449d 100644 --- a/src/npair_half_size_multi_old_newton_tri.cpp +++ b/src/npair_half_size_multi_old_newton_tri.cpp @@ -15,6 +15,8 @@ #include "npair_half_size_multi_old_newton_tri.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" #include "molecule.h" #include "my_page.h" @@ -35,7 +37,8 @@ NPairHalfSizeMultiOldNewtonTri::NPairHalfSizeMultiOldNewtonTri(LAMMPS *lmp) : NP void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,ns; + int i,j,jh,k,n,itype,jtype,ibin,ns,which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -45,17 +48,26 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; if (includegroup) nlocal = atom->nfirst; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -69,7 +81,11 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; - + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in bins, including self, in stencil // skip if i,j neighbor cutoff is less than bin distance @@ -107,10 +123,24 @@ void NPairHalfSizeMultiOldNewtonTri::build(NeighList *list) cutdistsq = (radsum+skin) * (radsum+skin); if (rsq <= cutdistsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; } } } diff --git a/src/npair_half_size_nsq_newtoff.cpp b/src/npair_half_size_nsq_newtoff.cpp index e8eac6f54b..0a127f3886 100644 --- a/src/npair_half_size_nsq_newtoff.cpp +++ b/src/npair_half_size_nsq_newtoff.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_nsq_newtoff.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "group.h" #include "my_page.h" #include "neigh_list.h" @@ -35,7 +38,8 @@ NPairHalfSizeNsqNewtoff::NPairHalfSizeNsqNewtoff(LAMMPS *lmp) : NPair(lmp) {} void NPairHalfSizeNsqNewtoff::build(NeighList *list) { - int i,j,n,bitmask; + int i,j,jh,n,bitmask,which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -44,7 +48,10 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list) double *radius = atom->radius; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; if (includegroup) { @@ -52,13 +59,19 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list) bitmask = group->bitmask[includegroup]; } + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -71,6 +84,11 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over remaining atoms, owned and ghost @@ -86,10 +104,24 @@ void NPairHalfSizeNsqNewtoff::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; } } diff --git a/src/npair_half_size_nsq_newton.cpp b/src/npair_half_size_nsq_newton.cpp index 7f90596755..866165c12e 100644 --- a/src/npair_half_size_nsq_newton.cpp +++ b/src/npair_half_size_nsq_newton.cpp @@ -15,7 +15,10 @@ #include "npair_half_size_nsq_newton.h" #include "atom.h" +#include "atom_vec.h" +#include "domain.h" #include "error.h" +#include "molecule.h" #include "group.h" #include "my_page.h" #include "neigh_list.h" @@ -36,7 +39,8 @@ NPairHalfSizeNsqNewton::NPairHalfSizeNsqNewton(LAMMPS *lmp) : NPair(lmp) {} void NPairHalfSizeNsqNewton::build(NeighList *list) { - int i,j,n,itag,jtag,bitmask; + int i,j,jh,n,itag,jtag,bitmask,which,imol,iatom,moltemplate; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutsq; int *neighptr; @@ -47,6 +51,8 @@ void NPairHalfSizeNsqNewton::build(NeighList *list) int *type = atom->type; int *mask = atom->mask; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; int nlocal = atom->nlocal; int nall = nlocal + atom->nghost; if (includegroup) { @@ -54,13 +60,19 @@ void NPairHalfSizeNsqNewton::build(NeighList *list) bitmask = group->bitmask[includegroup]; } + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + if (molecular == Atom::TEMPLATE) moltemplate = 1; + else moltemplate = 0; + int history = list->history; int *ilist = list->ilist; int *numneigh = list->numneigh; int **firstneigh = list->firstneigh; MyPage *ipage = list->ipage; - int mask_history = 3 << SBBITS; + int mask_history = 1 << HISTBITS; int inum = 0; ipage->reset(); @@ -74,6 +86,11 @@ void NPairHalfSizeNsqNewton::build(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over remaining atoms, owned and ghost @@ -105,10 +122,24 @@ void NPairHalfSizeNsqNewton::build(NeighList *list) cutsq = (radsum+skin) * (radsum+skin); if (rsq <= cutsq) { + jh = j; if (history && rsq < radsum*radsum) - neighptr[n++] = j ^ mask_history; - else - neighptr[n++] = j; + jh = jh ^ mask_history; + + if (molecular != Atom::ATOMIC) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >= 0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; + if (which == 0) neighptr[n++] = jh; + else if (domain->minimum_image_check(delx,dely,delz)) + neighptr[n++] = jh; + else if (which > 0) neighptr[n++] = jh ^ (which << SBBITS); + // OLD: if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = jh; } } diff --git a/src/set.cpp b/src/set.cpp index 2966106345..5f0d1dbec1 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -303,7 +303,7 @@ void Set::command(int narg, char **arg) else zvalue = utils::numeric(FLERR,arg[iarg+3],false,lmp); if (utils::strmatch(arg[iarg+4],"^v_")) varparse(arg[iarg+4],4); else wvalue = utils::numeric(FLERR,arg[iarg+4],false,lmp); - if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag) + if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag && !atom->quat_flag) error->all(FLERR,"Cannot set this attribute for this atom style"); set(QUAT); iarg += 5; @@ -311,7 +311,7 @@ void Set::command(int narg, char **arg) } else if (strcmp(arg[iarg],"quat/random") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal set command"); ivalue = utils::inumeric(FLERR,arg[iarg+1],false,lmp); - if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag) + if (!atom->ellipsoid_flag && !atom->tri_flag && !atom->body_flag && !atom->quat_flag) error->all(FLERR,"Cannot set this attribute for this atom style"); if (ivalue <= 0) error->all(FLERR,"Invalid random number seed in set command"); @@ -937,18 +937,20 @@ void Set::set(int keyword) sp[i][3] = dvalue; } - // set quaternion orientation of ellipsoid or tri or body particle - // set quaternion orientation of ellipsoid or tri or body particle + // set quaternion orientation of ellipsoid or tri or body particle or sphere/bpm // enforce quat rotation vector in z dir for 2d systems else if (keyword == QUAT) { double *quat = nullptr; + double **quat2 = nullptr; if (avec_ellipsoid && atom->ellipsoid[i] >= 0) quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat; else if (avec_tri && atom->tri[i] >= 0) quat = avec_tri->bonus[atom->tri[i]].quat; else if (avec_body && atom->body[i] >= 0) quat = avec_body->bonus[atom->body[i]].quat; + else if (atom->quat_flag) + quat2 = atom->quat; else error->one(FLERR,"Cannot set quaternion for atom that has none"); if (domain->dimension == 2 && (xvalue != 0.0 || yvalue != 0.0)) @@ -957,11 +959,24 @@ void Set::set(int keyword) double theta2 = MY_PI2 * wvalue/180.0; double sintheta2 = sin(theta2); - quat[0] = cos(theta2); - quat[1] = xvalue * sintheta2; - quat[2] = yvalue * sintheta2; - quat[3] = zvalue * sintheta2; - MathExtra::qnormalize(quat); + if (atom->quat_flag) { + double temp[4]; + temp[0] = cos(theta2); + temp[1] = xvalue * sintheta2; + temp[2] = yvalue * sintheta2; + temp[3] = zvalue * sintheta2; + MathExtra::qnormalize(temp); + quat2[i][0] = temp[0]; + quat2[i][1] = temp[1]; + quat2[i][2] = temp[2]; + quat2[i][3] = temp[3]; + } else { + quat[0] = cos(theta2); + quat[1] = xvalue * sintheta2; + quat[2] = yvalue * sintheta2; + quat[3] = zvalue * sintheta2; + MathExtra::qnormalize(quat); + } } // set theta of line particle @@ -1223,6 +1238,7 @@ void Set::setrandom(int keyword) } else if (keyword == QUAT_RANDOM) { int nlocal = atom->nlocal; double *quat; + double **quat2; if (domain->dimension == 3) { double s,t1,t2,theta1,theta2; @@ -1234,6 +1250,8 @@ void Set::setrandom(int keyword) quat = avec_tri->bonus[atom->tri[i]].quat; else if (avec_body && atom->body[i] >= 0) quat = avec_body->bonus[atom->body[i]].quat; + else if (atom->quat_flag) + quat2 = atom->quat; else error->one(FLERR,"Cannot set quaternion for atom that has none"); @@ -1243,10 +1261,17 @@ void Set::setrandom(int keyword) t2 = sqrt(s); theta1 = 2.0*MY_PI*ranpark->uniform(); theta2 = 2.0*MY_PI*ranpark->uniform(); - quat[0] = cos(theta2)*t2; - quat[1] = sin(theta1)*t1; - quat[2] = cos(theta1)*t1; - quat[3] = sin(theta2)*t2; + if (atom->quat_flag) { + quat2[i][0] = cos(theta2)*t2; + quat2[i][1] = sin(theta1)*t1; + quat2[i][2] = cos(theta1)*t1; + quat2[i][3] = sin(theta2)*t2; + } else { + quat[0] = cos(theta2)*t2; + quat[1] = sin(theta1)*t1; + quat[2] = cos(theta1)*t1; + quat[3] = sin(theta2)*t2; + } count++; } @@ -1258,15 +1283,24 @@ void Set::setrandom(int keyword) quat = avec_ellipsoid->bonus[atom->ellipsoid[i]].quat; else if (avec_body && atom->body[i] >= 0) quat = avec_body->bonus[atom->body[i]].quat; + else if (atom->quat_flag) + quat2 = atom->quat; else error->one(FLERR,"Cannot set quaternion for atom that has none"); ranpark->reset(seed,x[i]); theta2 = MY_PI*ranpark->uniform(); - quat[0] = cos(theta2); - quat[1] = 0.0; - quat[2] = 0.0; - quat[3] = sin(theta2); + if (atom->quat_flag) { + quat2[i][0] = cos(theta2); + quat2[i][1] = 0.0; + quat2[i][2] = 0.0; + quat2[i][3] = sin(theta2); + } else { + quat[0] = cos(theta2); + quat[1] = 0.0; + quat[2] = 0.0; + quat[3] = sin(theta2); + } count++; } }