ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3421 from phankl/mesocnt_stable

Browse Source 
Major update to mesocnt styles
This commit is contained in:
Axel Kohlmeyer
2022-08-31 11:54:22 -04:00
committed by GitHub
parent dd4260da77 0c23ba71df
commit 22f7255e1d
28 changed files with 241044 additions and 2219 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
doc/src/Commands_bond.rst
View File

@ -44,6 +44,7 @@ OPT.
* :doc:`harmonic (iko) <bond_harmonic>`
* :doc:`harmonic/shift (o) <bond_harmonic_shift>`
* :doc:`harmonic/shift/cut (o) <bond_harmonic_shift_cut>`
* :doc:`mesocnt <bond_mesocnt>`
* :doc:`mm3 <bond_mm3>`
* :doc:`morse (o) <bond_morse>`
* :doc:`nonlinear (o) <bond_nonlinear>`
@ -92,6 +93,7 @@ OPT.
* :doc:`fourier/simple (o) <angle_fourier_simple>`
* :doc:`gaussian <angle_gaussian>`
* :doc:`harmonic (iko) <angle_harmonic>`
* :doc:`mesocnt <angle_mesocnt>`
* :doc:`mm3 <angle_mm3>`
* :doc:`quartic (o) <angle_quartic>`
* :doc:`spica (o) <angle_spica>`

1
doc/src/Commands_pair.rst
View File

@ -201,6 +201,7 @@ OPT.
* :doc:`meam/spline (o) <pair_meam_spline>`
* :doc:`meam/sw/spline <pair_meam_sw_spline>`
* :doc:`mesocnt <pair_mesocnt>`
* :doc:`mesocnt/viscous <pair_mesocnt>`
* :doc:`mesont/tpm <pair_mesont_tpm>`
* :doc:`mgpt <pair_mgpt>`
* :doc:`mie/cut (g) <pair_mie>`

39
doc/src/Packages_details.rst
View File

@ -1556,31 +1556,40 @@ MESONT package
**Contents:**
MESONT is a LAMMPS package for simulation of nanomechanics of
nanotubes (NTs). The model is based on a coarse-grained representation
of NTs as "flexible cylinders" consisting of a variable number of
MESONT is a LAMMPS package for simulation of nanomechanics of nanotubes
(NTs). The model is based on a coarse-grained representation of NTs as
"flexible cylinders" consisting of a variable number of
segments. Internal interactions within a NT and the van der Waals
interaction between the tubes are described by a mesoscopic force field
designed and parameterized based on the results of atomic-level
molecular dynamics simulations. The description of the force field is
provided in the papers listed below. This package contains two
independent implementations of this model: :doc:`pair_style mesocnt
<pair_mesocnt>` is a (minimal) C++ implementation, and :doc:`pair_style
mesont/tpm <pair_mesont_tpm>` is a more general and feature rich
implementation based on a Fortran library in the ``lib/mesont`` folder.
provided in the papers listed below.
This package contains two independent implementations of this model:
:doc:`pair_style mesont/tpm <pair_mesont_tpm>` is the original
implementation of the model based on a Fortran library in the
``lib/mesont`` folder. The second implementation is provided by the
mesocnt styles (:doc:`bond_style mesocnt <bond_mesocnt>`,
:doc:`angle_style mesocnt <angle_mesocnt>` and :doc:`pair_style mesocnt
<pair_mesocnt>`). The mesocnt implementation has the same features as
the original implementation with the addition of friction, but is
directly implemented in C++, interfaces more cleanly with general LAMMPS
functionality, and is typically faster. It also does not require its own
atom style and can be installed without any external libraries.
**Download of potential files:**
The potential files for these pair styles are *very* large and thus
are not included in the regular downloaded packages of LAMMPS or the
git repositories. Instead, they will be automatically downloaded
from a web server when the package is installed for the first time.
The potential files for these pair styles are *very* large and thus are
not included in the regular downloaded packages of LAMMPS or the git
repositories. Instead, they will be automatically downloaded from a web
server when the package is installed for the first time.
**Authors of the *mesont* styles:**
Maxim V. Shugaev (University of Virginia), Alexey N. Volkov (University of Alabama), Leonid V. Zhigilei (University of Virginia)
Maxim V. Shugaev (University of Virginia), Alexey N. Volkov (University
of Alabama), Leonid V. Zhigilei (University of Virginia)
**Author of the *mesocnt* pair style:**
**Author of the *mesocnt* styles:**
Philipp Kloza (U Cambridge)
**Supporting info:**
@ -1590,6 +1599,8 @@ Philipp Kloza (U Cambridge)
* :doc:`atom_style mesont <atom_style>`
* :doc:`pair_style mesont/tpm <pair_mesont_tpm>`
* :doc:`compute mesont <compute_mesont>`
* :doc:`bond_style mesocnt <bond_mesocnt>`
* :doc:`angle_style mesocnt <angle_mesocnt>`
* :doc:`pair_style mesocnt <pair_mesocnt>`
* examples/PACKAGES/mesont
* tools/mesont

146
doc/src/angle_mesocnt.rst Normal file
View File

@ -0,0 +1,146 @@
.. index:: angle_style mesocnt
angle_style mesocnt command
===========================
Syntax
""""""
.. code-block:: LAMMPS
angle_style mesocnt
Examples
""""""""
.. code-block:: LAMMPS
angle_style mesocnt
angle_coeff 1 buckling C 10 10 20.0
angle_coeff 4 harmonic C 8 4 10.0
angle_coeff 2 buckling custom 400.0 50.0 5.0
angle_coeff 1 harmonic custom 300.0
Description
"""""""""""
.. versionadded:: TBD
The *mesocnt* angle style uses the potential
.. math::
E = K_\text{H} \Delta \theta^2, \qquad |\Delta \theta| < \Delta
\theta_\text{B} \\
E = K_\text{H} \Delta \theta_\text{B}^2 +
K_\text{B} (\Delta \theta - \Delta \theta_\text{B}), \qquad |\Delta
\theta| \geq \Delta \theta_\text{B}
where :math:`\Delta \theta = \theta - \pi` is the bending angle of the
nanotube, :math:`K_\text{H}` and :math:`K_\text{B}` are prefactors for
the harmonic and linear regime respectively and :math:`\Delta
\theta_\text{B}` is the buckling angle. Note that the usual 1/2 factor
for the harmonic potential is included in :math:`K_\text{H}`.
The style implements parameterization presets of :math:`K_\text{H}`,
:math:`K_\text{B}` and :math:`\Delta \theta_\text{B}` for mesoscopic
simulations of carbon nanotubes based on the atomistic simulations of
:ref:`(Srivastava) <Srivastava_2>` and buckling considerations of
:ref:`(Zhigilei) <Zhigilei1_1>`.
The following coefficients must be defined for each angle type via the
:doc:`angle_coeff <angle_coeff>` command as in the examples above, or
in the data file or restart files read by the :doc:`read_data
<read_data>` or :doc:`read_restart <read_restart>` commands:
* mode = *buckling* or *harmonic*
* preset = *C* or *custom*
* additional parameters depending on preset
If mode *harmonic* is chosen, the potential is simply harmonic and
does not switch to the linear term when the buckling angle is
reached. In *buckling* mode, the full piecewise potential is used.
Preset *C* is for carbon nanotubes, and the additional parameters are:
* chiral index :math:`n` (unitless)
* chiral index :math:`m` (unitless)
* :math:`r_0` (distance)
Here, :math:`r_0` is the equilibrium distance of the bonds included in
the angle, see :doc:`bond_style mesocnt <bond_mesocnt>`.
In harmonic mode with preset *custom*, the additional parameter is:
* :math:`K_\text{H}` (energy)
Hence, this setting is simply a wrapper for :doc:`bond_style harmonic
<bond_harmonic>` with an equilibrium angle of 180 degrees.
In harmonic mode with preset *custom*, the additional parameters are:
* :math:`K_\text{H}` (energy)
* :math:`K_\text{B}` (energy)
* :math:`\Delta \theta_\text{B}` (degrees)
:math:`\Delta \theta_\text{B}` is specified in degrees, but LAMMPS
converts it to radians internally; hence :math:`K_\text{H}` is
effectively energy per radian\^2 and :math:`K_\text{B}` is energy per
radian.
----------
In *buckling* mode, this angle style adds the *buckled* property to
all atoms in the simulation, which is an integer flag indicating
whether the bending angle at a given atom has exceeded :math:`\Delta
\theta_\text{B}`. It can be accessed as an atomic variable, e.g. for
custom dump commands, as *i_buckled*.
.. note::
If the initial state of the simulation contains buckled nanotubes
and :doc:`pair_style mesocnt <pair_mesocnt>` is used, the
*i_buckled* atomic variable needs to be initialized before the
pair_style is defined by doing a *run 0* command straight after the
angle_style command. See below for an example.
If CNTs are already buckled at the start of the simulation, this
script will correctly initialize *i_buckled*:
.. code-block:: LAMMPS
angle_style mesocnt
angle_coeff 1 buckling C 10 10 20.0
run 0
pair_style mesocnt 60.0
pair_coeff * * C_10_10.mesocnt 1
Restrictions
""""""""""""
This angle style can only be used if LAMMPS was built with the
MOLECULE and MESONT packages. See the :doc:`Build package
<Build_package>` doc page for more info.
Related commands
""""""""""""""""
:doc:`angle_coeff <angle_coeff>`
Default
"""""""
none
----------
.. _Srivastava_2:
**(Srivastava)** Zhigilei, Wei, Srivastava, Phys. Rev. B 71, 165417
(2005).
.. _Zhigilei1_1:
**(Zhigilei)** Volkov and Zhigilei, ACS Nano 4, 6187 (2010).

1
doc/src/angle_style.rst
View File

@ -90,6 +90,7 @@ of (g,i,k,o,t) to indicate which accelerated styles exist.
* :doc:`fourier/simple <angle_fourier_simple>` - angle with a single cosine term
* :doc:`gaussian <angle_gaussian>` - multi-centered Gaussian-based angle potential
* :doc:`harmonic <angle_harmonic>` - harmonic angle
* :doc:`mesocnt <angle_mesocnt>` - piecewise harmonic and linear angle for bending-buckling of nanotubes
* :doc:`mm3 <angle_mm3>` - anharmonic angle
* :doc:`quartic <angle_quartic>` - angle with cubic and quartic terms
* :doc:`spica <angle_spica>` - harmonic angle with repulsive SPICA pair style between 1-3 atoms

84
doc/src/bond_mesocnt.rst Normal file
View File

@ -0,0 +1,84 @@
.. index:: bond_style mesocnt
bond_style mesocnt command
===========================
Syntax
""""""
.. code-block:: LAMMPS
bond_style mesocnt
Examples
""""""""
.. code-block:: LAMMPS
bond_style mesocnt
bond_coeff 1 C 10 10 20.0
bond_coeff 4 custom 800.0 10.0
Description
"""""""""""
.. versionadded:: TBD
The *mesocnt* bond style is a wrapper for the :doc:`harmonic
<bond_harmonic>` style, and uses the potential
.. math::
E = K (r - r_0)^2
where :math:`r_0` is the equilibrium bond distance. Note that the
usual 1/2 factor is included in :math:`K`. The style implements
parameterization presets of :math:`K` for mesoscopic simulations of
carbon nanotubes based on the atomistic simulations of
:ref:`(Srivastava) <Srivastava_1>`.
Other presets can be readily implemented in the future.
The following coefficients must be defined for each bond type via the
:doc:`bond_coeff <bond_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data
<read_data>` or :doc:`read_restart <read_restart>` commands:
* preset = *C* or *custom*
* additional parameters depending on preset
Preset *C* is for carbon nanotubes, and the additional parameters are:
* chiral index :math:`n` (unitless)
* chiral index :math:`m` (unitless)
* :math:`r_0` (distance)
Preset *custom* is simply a direct wrapper for the :doc:`harmonic
<bond_harmonic>` style, and the additional parameters are:
* :math:`K` (energy/distance\^2)
* :math:`r_0` (distance)
Restrictions
""""""""""""
This bond style can only be used if LAMMPS was built with the MOLECULE
and MESONT packages. See the :doc:`Build package <Build_package>`
page for more info.
Related commands
""""""""""""""""
:doc:`bond_coeff <bond_coeff>`, :doc:`delete_bonds <delete_bonds>`
Default
"""""""
none
----------
.. _Srivastava_1:
**(Srivastava)** Zhigilei, Wei and Srivastava, Phys. Rev. B 71, 165417
(2005).

1
doc/src/bond_style.rst
View File

@ -95,6 +95,7 @@ accelerated styles exist.
* :doc:`harmonic <bond_harmonic>` - harmonic bond
* :doc:`harmonic/shift <bond_harmonic_shift>` - shifted harmonic bond
* :doc:`harmonic/shift/cut <bond_harmonic_shift_cut>` - shifted harmonic bond with a cutoff
* :doc:`mesocnt <bond_mesocnt>` - Harmonic bond wrapper with parameterization presets for nanotubes
* :doc:`mm3 <bond_mm3>` - MM3 anharmonic bond
* :doc:`morse <bond_morse>` - Morse bond
* :doc:`nonlinear <bond_nonlinear>` - nonlinear bond

242
doc/src/pair_mesocnt.rst
View File

@ -1,68 +1,153 @@
.. index:: pair_style mesocnt
.. index:: pair_style mesocnt/viscous
pair_style mesocnt command
==========================
pair_style mesocnt/viscous command
==================================
Syntax
""""""
.. code-block:: LAMMPS
pair_style mesocnt
pair_style style neigh_cutoff mode neigh_mode
* style = *mesocnt* or *mesocnt/viscous*
* neigh_cutoff = neighbor list cutoff (distance units)
* mode = *chain* or *segment* (optional)
* neigh_mode = *id* or *topology* (optional)
Examples
""""""""
.. code-block:: LAMMPS
pair_style mesocnt
pair_coeff * * 10_10.cnt
pair_style mesocnt 30.0
pair_coeff * * C_10_10.mesocnt 2
pair_style mesocnt/viscous 60.0 chain topology
pair_coeff * * C_10_10.mesocnt 0.001 20.0 0.2 2 4
Description
"""""""""""
Style *mesocnt* implements a mesoscopic potential
for the interaction of carbon nanotubes (CNTs). In this potential,
CNTs are modelled as chains of cylindrical segments in which
each infinitesimal surface element interacts with all other
CNT surface elements with the Lennard-Jones (LJ) term adopted from
the :doc:`airebo <pair_airebo>` style. The interaction energy
is then computed by integrating over the surfaces of all interacting
CNTs.
Style *mesocnt* implements a mesoscopic potential for the interaction
of carbon nanotubes (CNTs), or other quasi-1D objects such as other
kinds of nanotubes or nanowires. In this potential, CNTs are modelled
as chains of cylindrical segments in which each infinitesimal surface
element interacts with all other CNT surface elements with the
Lennard-Jones (LJ) term adopted from the :doc:`airebo <pair_airebo>`
style. The interaction energy is then computed by integrating over the
surfaces of all interacting CNTs.
The potential is based on interactions between one cylindrical
segment and infinitely or semi-infinitely long CNTs as described
in :ref:`(Volkov1) <Volkov1>`. Chains of segments are
converted to these (semi-)infinite CNTs bases on an approximate
chain approach outlined in :ref:`(Volkov2) <Volkov2>`.
This allows to simplify the computation of the interactions
significantly and reduces the computational times to the
same order of magnitude as for regular bead spring models
where beads interact with the standard :doc:`pair_lj/cut <pair_lj>`
potential.
In LAMMPS, cylindrical segments are represented by bonds. Each segment
is defined by its two end points ("nodes") which correspond to atoms
in LAMMPS. For the exact functional form of the potential and
implementation details, the reader is referred to the original papers
:ref:`(Volkov1) <Volkov1>` and :ref:`(Volkov2) <Volkov2>`.
In LAMMPS, cylindrical segments are represented by bonds. Each
segment is defined by its two end points ("nodes") which correspond
to atoms in LAMMPS. For the exact functional form of the potential
and implementation details, the reader is referred to the
original papers :ref:`(Volkov1) <Volkov1>` and
:ref:`(Volkov2) <Volkov2>`.
.. versionchanged:: TBD
The potential requires tabulated data provided in a single ASCII
text file specified in the :doc:`pair_coeff <pair_coeff>` command.
The first line of the file provides a time stamp and
general information. The second line lists four integers giving
the number of data points provided in the subsequent four
data tables. The third line lists four floating point numbers:
the CNT radius R, the LJ parameter sigma and two numerical
parameters delta1 and delta2. These four parameters are given
in Angstroms. This is followed by four data tables each separated
by a single empty line. The first two tables have two columns
and list the parameters uInfParallel and Gamma respectively.
The last two tables have three columns giving data on a quadratic
array and list the parameters Phi and uSemiParallel respectively.
uInfParallel and uSemiParallel are given in eV/Angstrom, Phi is
given in eV and Gamma is unitless.
The potential supports two modes, *segment* and *chain*. By default,
*chain* mode is enabled. In *segment* mode, interactions are
pair-wise between all neighboring segments based on a segment-segment
approach (keyword *segment* in pair_style command). In *chain* mode,
interactions are calculated between each segment and infinitely or
semi-infinitely long CNTs as described in :ref:`(Volkov1) <Volkov1>`.
Chains of segments are converted to these (semi-)infinite CNTs bases
on an approximate chain approach outlined in :ref:`(Volkov2)
<Volkov2>`. Hence, interactions are calculated on a segment-chain
basis (keyword *chain* in the pair_style command). Using *chain* mode
allows to simplify the computation of the interactions significantly
and reduces the computational times to the same order of magnitude as
for regular bead spring models where beads interact with the standard
:doc:`pair_lj/cut <pair_lj>` potential. However, this method is only
valid when the curvature of the CNTs in the system is small. When
CNTs are buckled (see :doc:`angle_mesocnt <angle_mesocnt>`), local
curvature can be very high and the pair_style automatically switches
to *segment* mode for interactions involving buckled CNTs.
The potential further implements two different neighbor list
construction modes. Mode *id* uses atom and mol IDs to construct
neighbor lists while *topology* modes uses only the bond topology of
the system. While *id* mode requires bonded atoms to have consecutive
LAMMPS atom IDs and atoms in different CNTs to have different LAMMPS
molecule IDs, *topology* mode has no such requirement. Using *id* mode
is faster and is enabled by default.
.. note::
Neighbor *id* mode requires all CNTs in the system to have distinct
LAMMPS molecule IDs and bonded atoms to have consecutive LAMMPS atom
IDs. If this is not possible (e.g. in simulations of CNT rings),
*topology* mode needs to be enabled in the pair_style command.
.. versionadded:: TBD
In addition to the LJ interactions described above, style
*mesocnt/viscous* explicitly models friction between neighboring
segments. Friction forces are a function of the relative velocity
between a segment and its neighboring approximate chain (even in
*segment* mode) and only act along the axes of the interacting segment
and chain. In this potential, friction forces acting per unit length
of a nanotube segent are modelled as a shifted logistic function:
.. math::
F^{\text{FRICTION}}(v) / L = \frac{F^{\text{max}}}{1 +
\exp(-k(v-v_0))} - \frac{F^{\text{max}}}{1 + \exp(k v_0)}
----------
In the pair_style command, the modes described above can be toggled
using the *segment* or *chain* keywords. The neighbor list cutoff
defines the cutoff within which atoms are included in the neighbor
list for constructing neighboring CNT chains. This is different from
the potential cutoff, which is directly calculated from parameters
specified in the potential file. We recommend using a neighbor list
cutoff of at least 3 times the maximum segment length used in the
simulation to ensure proper neighbor chain construction.
.. note::
CNT ends are treated differently by all *mesocnt* styles. Atoms on
CNT ends need to be assigned different LAMMPS atom types than atoms
not on CNT ends.
Style *mesocnt* requires tabulated data provided in a single ASCII
text file, as well as a list of integers corresponding to all LAMMPS
atom types representing CNT ends:
* filename
* :math:`N` CNT end atom types
For example, if your LAMMPS simulation of (10, 10) nanotubes has 4
atom types where atom types 1 and 3 are assigned to 'inner' nodes and
atom types 2 and 4 are assigned to CNT end nodes, the pair_coeff
command would be:
.. code-block:: LAMMPS
pair_coeff * * C_10_10.mesocnt 2 4
Likewise, style *mesocnt/viscous* also requires the same information
as style *mesocnt*, with the addition of 3 parameters for the viscous
friction forces as listed above:
* filename
* :math:`F^{\text{max}}`
* :math:`k`
* :math:`v_0`
* :math:`N` CNT end atom types
Using the same example system as with style *mesocnt* with the
addition of friction, the pair_coeff command is:
.. code-block:: LAMMPS
pair_coeff * * C_10_10.mesocnt 0.03 20.0 0.20 2 4
Potential files for CNTs can be readily generated using the freely
available code provided on
@ -71,45 +156,51 @@ available code provided on
https://github.com/phankl/cntpot
Using the same approach, it should also be possible to
generate potential files for other 1D systems such as
boron nitride nanotubes.
Using the same approach, it should also be possible to generate
potential files for other 1D systems mentioned above.
.. note::
Because of their size, *mesocnt* style potential files
are not bundled with LAMMPS. When compiling LAMMPS from
source code, the file ``C_10_10.mesocnt`` should be downloaded
transparently from `https://download.lammps.org/potentials/C_10_10.mesocnt <https://download.lammps.org/potentials/C_10_10.mesocnt>`_
This file has as number of data points per table 1001.
This is sufficient for NVT simulations. For proper energy
conservation, we recommend using a potential file where
the resolution for Phi is at least 2001 data points.
Because of their size, *mesocnt* style potential files are not
bundled with LAMMPS. When compiling LAMMPS from source code, the
file ``C_10_10.mesocnt`` should be downloaded separately from
`https://download.lammps.org/potentials/C_10_10.mesocnt
<https://download.lammps.org/potentials/C_10_10.mesocnt>`_
.. note::
The first line of the potential file provides a time stamp and
general information. The second line lists four integers giving the
number of data points provided in the subsequent four data
tables. The third line lists four floating point numbers: the CNT
radius R, the LJ parameter sigma and two numerical parameters
delta1 and delta2. These four parameters are given in
Angstroms. This is followed by four data tables each separated by a
single empty line. The first two tables have two columns and list
the parameters uInfParallel and Gamma respectively. The last two
tables have three columns giving data on a quadratic array and list
the parameters Phi and uSemiParallel respectively. uInfParallel
and uSemiParallel are given in eV/Angstrom, Phi is given in eV and
Gamma is unitless.
The *mesocnt* style requires CNTs to be represented
as a chain of atoms connected by bonds. Atoms need
to be numbered consecutively within one chain.
Atoms belonging to different CNTs need to be assigned
different molecule IDs.
If a simulation produces many warnings about segment-chain
interactions falling outside the interpolation range, we recommend
generating a potential file with lower values of delta1 and delta2.
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This pair style does not support mixing.
These pair styles does not support mixing.
This pair style does not support the :doc:`pair_modify <pair_modify>`
shift, table, and tail options.
These pair styles does not support the :doc:`pair_modify
<pair_modify>` shift, table, and tail options.
The *mesocnt* pair style do not write their information to :doc:`binary restart files <restart>`,
since it is stored in tabulated potential files.
Thus, you need to re-specify the pair_style and pair_coeff commands in
an input script that reads a restart file.
These pair styles do not write their information to :doc:`binary
restart files <restart>`, since it is stored in tabulated potential
files. Thus, you need to re-specify the pair_style and pair_coeff
commands in an input script that reads a restart file.
This pair style can only be used via the *pair* keyword of the
These pair styles can only be used via the *pair* keyword of the
:doc:`run_style respa <run_style>` command. They do not support the
*inner*, *middle*, *outer* keywords.
@ -118,12 +209,21 @@ This pair style can only be used via the *pair* keyword of the
Restrictions
""""""""""""
This style is part of the MESONT package. It is only
enabled if LAMMPS was built with that package. See the :doc:`Build package <Build_package>` page for more info.
These styles are part of the MESONT package. They are only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
This pair potential requires the :doc:`newton <newton>` setting to be
These pair styles require the :doc:`newton <newton>` setting to be
"on" for pair interactions.
These pair styles require all 3 :doc:`special_bonds lj
<special_bonds>` settings to be non-zero for proper neighbor list
construction.
Pair style *mesocnt/viscous* requires you to use the :doc:`comm_modify
vel yes <comm_modify>` command so that velocities are stored by ghost
atoms.
Related commands
""""""""""""""""
@ -132,7 +232,7 @@ Related commands
Default
"""""""
none
mode = chain, neigh_mode = id
----------

3
doc/src/pair_style.rst
View File

@ -278,7 +278,8 @@ accelerated styles exist.
* :doc:`meam <pair_meam>` - modified embedded atom method (MEAM) in C
* :doc:`meam/spline <pair_meam_spline>` - splined version of MEAM
* :doc:`meam/sw/spline <pair_meam_sw_spline>` - splined version of MEAM with a Stillinger-Weber term
* :doc:`mesocnt <pair_mesocnt>` - mesoscale model for (carbon) nanotubes
* :doc:`mesocnt <pair_mesocnt>` - mesoscopic vdW potential for (carbon) nanotubes
* :doc:`mesocnt/viscous <pair_mesocnt>` - mesoscopic vdW potential for (carbon) nanotubes with friction
* :doc:`mgpt <pair_mgpt>` - simplified model generalized pseudopotential theory (MGPT) potential
* :doc:`mesont/tpm <pair_mesont_tpm>` - nanotubes mesoscopic force field
* :doc:`mie/cut <pair_mie>` - Mie potential

1
doc/utils/sphinx-config/false_positives.txt
View File

@ -2250,6 +2250,7 @@ nanoparticles
nanotube
Nanotube
nanotubes
nanowires
Narulkar
nasa
nasr

1521
examples/PACKAGES/mesont/cnt.data
View File

File diff suppressed because it is too large Load Diff

237628
examples/PACKAGES/mesont/data.film_mesocnt Normal file
View File

File diff suppressed because it is too large Load Diff

38
examples/PACKAGES/mesont/in.cnt
View File

@ -1,38 +0,0 @@
#Initialisation
units nano
dimension 3
boundary p p p
atom_style full
comm_modify cutoff 11.0
neighbor 7.80 bin
newton on
#Read data
read_data cnt.data
replicate 1 2 2
#Force field
bond_style harmonic
bond_coeff 1 268896.77 2.0
angle_style harmonic
angle_coeff 1 46562.17 180.0
pair_style mesocnt
pair_coeff * * C_10_10.mesocnt
#Output
thermo 1000
dump xyz all xyz 1000 cnt.xyz
#Simulation setup
timestep 1.0e-05
#Nose-Hoover thermostat
fix nvt all nvt temp 300 300 0.001
run 10000

46
examples/PACKAGES/mesont/in.film_mesocnt Normal file
View File

@ -0,0 +1,46 @@
# initialisation
units metal
dimension 3
boundary p p s
atom_style full
special_bonds lj 1 1 1
neigh_modify every 5 delay 0 check yes
newton on
read_data data.film_mesocnt
# force field
bond_style mesocnt
bond_coeff 1 C 10 10 10
angle_style mesocnt
angle_coeff 1 buckling C 10 10 10
pair_style mesocnt 40 chain
pair_coeff * * C_10_10.mesocnt 1
# output
compute epair all pe pair
compute ebond all pe bond
compute eangle all pe angle
compute epair_atom all pe/atom pair
compute ebond_atom all pe/atom bond
compute eangle_atom all pe/atom angle
fix angle_mesocnt_buckled all property/atom i_buckled ghost yes
thermo_style custom step temp etotal ke pe c_ebond c_eangle c_epair
thermo 10
#dump custom all custom 100 film_mesocnt.lmp id mol type x y z c_ebond_atom c_eangle_atom c_epair_atom i_buckled
# simulation setup
velocity all create 600.0 2022 loop geom
timestep 0.01
fix nvt all nvt temp 300.0 300.0 1
run 100

126
examples/PACKAGES/mesont/log.31Aug2022.film_mesocnt.g++.1 Normal file
View File

@ -0,0 +1,126 @@
LAMMPS (3 Aug 2022)
# initialisation
units metal
dimension 3
boundary p p s
atom_style full
special_bonds lj 1 1 1
neigh_modify every 5 delay 0 check yes
newton on
read_data data.film_mesocnt
Reading data file ...
orthogonal box = (-2500 -2500 -300) to (2500 2500 402.42)
1 by 1 by 1 MPI processor grid
reading atoms ...
79596 atoms
scanning bonds ...
1 = max bonds/atom
scanning angles ...
1 = max angles/atom
reading bonds ...
79200 bonds
reading angles ...
78804 angles
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 1 1 1
special bond factors coul: 0 0 0
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.012 seconds
read_data CPU = 0.160 seconds
# force field
bond_style mesocnt
bond_coeff 1 C 10 10 10
angle_style mesocnt
angle_coeff 1 buckling C 10 10 10
pair_style mesocnt 40 chain
pair_coeff * * C_10_10.mesocnt 1
Reading mesocnt potential file C_10_10.mesocnt with DATE: 2022-07-12
# output
compute epair all pe pair
compute ebond all pe bond
compute eangle all pe angle
compute epair_atom all pe/atom pair
compute ebond_atom all pe/atom bond
compute eangle_atom all pe/atom angle
fix angle_mesocnt_buckled all property/atom i_buckled ghost yes
thermo_style custom step temp etotal ke pe c_ebond c_eangle c_epair
thermo 10
#dump custom all custom 100 film_mesocnt.lmp id mol type x y z c_ebond_atom c_eangle_atom c_epair_atom i_buckled
# simulation setup
velocity all create 600.0 2022 loop geom
timestep 0.01
fix nvt all nvt temp 300.0 300.0 1
run 100
Neighbor list info ...
update: every = 5 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 42
ghost atom cutoff = 42
binsize = 21, bins = 239 239 4
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair mesocnt, perpetual
attributes: full, newton on
pair build: full/bin
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 49.89 | 49.89 | 49.89 Mbytes
Step Temp TotEng KinEng PotEng c_ebond c_eangle c_epair
0 600 1355.8735 6173.0767 -4817.2032 28.668731 21.29199 -4867.1639
10 389.92732 1373.1839 4011.7521 -2638.5683 848.77485 1387.6166 -4874.9597
20 311.7468 1388.0161 3207.3949 -1819.3788 1211.534 1868.6817 -4899.5945
30 291.75586 1385.2114 3001.7189 -1616.5075 1184.1423 2133.2088 -4933.8586
40 320.42607 1364.2644 3296.6912 -1932.4267 951.62534 2088.33 -4972.382
50 341.37701 1346.0547 3512.2441 -2166.1894 956.59158 1891.9196 -5014.7005
60 333.15461 1337.5518 3427.6483 -2090.0965 1137.7331 1825.9946 -5053.8242
70 324.47061 1328.7153 3338.3033 -2009.588 1133.3639 1953.8288 -5096.7807
80 324.39487 1315.4948 3337.5241 -2022.0293 979.16213 2141.4697 -5142.6611
90 323.39973 1302.3471 3327.2856 -2024.9385 972.70286 2185.7686 -5183.41
100 322.73067 1289.0538 3320.402 -2031.3482 1100.1024 2088.3804 -5219.831
Loop time of 4.1637 on 1 procs for 100 steps with 79596 atoms
Performance: 20.751 ns/day, 1.157 hours/ns, 24.017 timesteps/s
99.5% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 3.705 | 3.705 | 3.705 | 0.0 | 88.98
Bond | 0.29317 | 0.29317 | 0.29317 | 0.0 | 7.04
Neigh | 0.078491 | 0.078491 | 0.078491 | 0.0 | 1.89
Comm | 0.0019462 | 0.0019462 | 0.0019462 | 0.0 | 0.05
Output | 0.00099817 | 0.00099817 | 0.00099817 | 0.0 | 0.02
Modify | 0.079874 | 0.079874 | 0.079874 | 0.0 | 1.92
Other | | 0.004224 | | | 0.10
Nlocal: 79596 ave 79596 max 79596 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 2567 ave 2567 max 2567 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 1.1337e+06 ave 1.1337e+06 max 1.1337e+06 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 1133696
Ave neighs/atom = 14.243128
Ave special neighs/atom = 5.9402985
Neighbor list builds = 2
Dangerous builds = 0
Total wall time: 0:00:05

126
examples/PACKAGES/mesont/log.31Aug2022.film_mesocnt.g++.4 Normal file
View File

@ -0,0 +1,126 @@
LAMMPS (3 Aug 2022)
# initialisation
units metal
dimension 3
boundary p p s
atom_style full
special_bonds lj 1 1 1
neigh_modify every 5 delay 0 check yes
newton on
read_data data.film_mesocnt
Reading data file ...
orthogonal box = (-2500 -2500 -300) to (2500 2500 402.42)
2 by 2 by 1 MPI processor grid
reading atoms ...
79596 atoms
scanning bonds ...
1 = max bonds/atom
scanning angles ...
1 = max angles/atom
reading bonds ...
79200 bonds
reading angles ...
78804 angles
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 1 1 1
special bond factors coul: 0 0 0
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.005 seconds
read_data CPU = 0.156 seconds
# force field
bond_style mesocnt
bond_coeff 1 C 10 10 10
angle_style mesocnt
angle_coeff 1 buckling C 10 10 10
pair_style mesocnt 40 chain
pair_coeff * * C_10_10.mesocnt 1
Reading mesocnt potential file C_10_10.mesocnt with DATE: 2022-07-12
# output
compute epair all pe pair
compute ebond all pe bond
compute eangle all pe angle
compute epair_atom all pe/atom pair
compute ebond_atom all pe/atom bond
compute eangle_atom all pe/atom angle
fix angle_mesocnt_buckled all property/atom i_buckled ghost yes
thermo_style custom step temp etotal ke pe c_ebond c_eangle c_epair
thermo 10
#dump custom all custom 100 film_mesocnt.lmp id mol type x y z c_ebond_atom c_eangle_atom c_epair_atom i_buckled
# simulation setup
velocity all create 600.0 2022 loop geom
timestep 0.01
fix nvt all nvt temp 300.0 300.0 1
run 100
Neighbor list info ...
update: every = 5 steps, delay = 0 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 42
ghost atom cutoff = 42
binsize = 21, bins = 239 239 4
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair mesocnt, perpetual
attributes: full, newton on
pair build: full/bin
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 16.47 | 16.58 | 16.89 Mbytes
Step Temp TotEng KinEng PotEng c_ebond c_eangle c_epair
0 600 1355.8735 6173.0767 -4817.2032 28.668731 21.29199 -4867.1639
10 389.92732 1373.1839 4011.7521 -2638.5683 848.77485 1387.6166 -4874.9597
20 311.7468 1388.0161 3207.3949 -1819.3788 1211.534 1868.6817 -4899.5945
30 291.75586 1385.2114 3001.7189 -1616.5075 1184.1423 2133.2088 -4933.8586
40 320.42607 1364.2644 3296.6912 -1932.4267 951.62534 2088.33 -4972.382
50 341.37701 1346.0547 3512.2441 -2166.1894 956.59158 1891.9196 -5014.7005
60 333.15461 1337.5518 3427.6483 -2090.0965 1137.7331 1825.9946 -5053.8242
70 324.47061 1328.7153 3338.3033 -2009.588 1133.3639 1953.8288 -5096.7807
80 324.39487 1315.4948 3337.5241 -2022.0293 979.16213 2141.4697 -5142.6611
90 323.39973 1302.3471 3327.2856 -2024.9385 972.70286 2185.7686 -5183.41
100 322.73067 1289.0538 3320.402 -2031.3482 1100.1024 2088.3804 -5219.831
Loop time of 1.25052 on 4 procs for 100 steps with 79596 atoms
Performance: 69.091 ns/day, 0.347 hours/ns, 79.967 timesteps/s
99.6% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.87036 | 0.97558 | 1.1135 | 9.0 | 78.01
Bond | 0.071751 | 0.076357 | 0.084244 | 1.7 | 6.11
Neigh | 0.023232 | 0.023239 | 0.023244 | 0.0 | 1.86
Comm | 0.0046002 | 0.15227 | 0.26319 | 24.1 | 12.18
Output | 0.00032696 | 0.00037811 | 0.00045537 | 0.0 | 0.03
Modify | 0.019263 | 0.020646 | 0.023155 | 1.1 | 1.65
Other | | 0.00204 | | | 0.16
Nlocal: 19899 ave 21951 max 18670 min
Histogram: 1 1 0 1 0 0 0 0 0 1
Nghost: 1323.5 ave 1412 max 1255 min
Histogram: 1 0 1 0 1 0 0 0 0 1
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 283424 ave 325466 max 258171 min
Histogram: 1 0 2 0 0 0 0 0 0 1
Total # of neighbors = 1133696
Ave neighs/atom = 14.243128
Ave special neighs/atom = 5.9402985
Neighbor list builds = 2
Dangerous builds = 0
Total wall time: 0:00:02

126
examples/PACKAGES/mesont/log.9Jan20.cnt.g++.1
View File

@ -1,126 +0,0 @@
LAMMPS (09 Jan 2020)
#Initialisation
units nano
dimension 3
boundary p p p
atom_style full
comm_modify cutoff 11.0
neighbor 7.80 bin
newton on
#Read data
read_data cnt.data
orthogonal box = (0 0 0) to (600 600 60)
1 by 1 by 1 MPI processor grid
reading atoms ...
500 atoms
scanning bonds ...
1 = max bonds/atom
scanning angles ...
1 = max angles/atom
reading bonds ...
498 bonds
reading angles ...
496 angles
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.000180006 secs
read_data CPU = 0.00125766 secs
replicate 1 2 2
orthogonal box = (0 0 0) to (600 1200 120)
1 by 1 by 1 MPI processor grid
2000 atoms
1992 bonds
1984 angles
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.00054121 secs
replicate CPU = 0.000902414 secs
#Force field
bond_style harmonic
bond_coeff 1 268896.77 2.0
angle_style harmonic
angle_coeff 1 46562.17 180.0
pair_style mesocnt
pair_coeff * * 10_10.cnt
Reading potential file 10_10.cnt with DATE: 2020-01-13
#Output
thermo 1000
dump xyz all xyz 1000 cnt.xyz
#Simulation setup
timestep 1.0e-05
#Nose-Hoover thermostat
fix nvt all nvt temp 300 300 0.001
run 10000
WARNING: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions (src/pair.cpp:226)
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.177
ghost atom cutoff = 11
binsize = 5.0885, bins = 118 236 24
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair mesocnt, perpetual
attributes: full, newton on
pair build: full/bin
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 11.21 | 11.21 | 11.21 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -4632.0136 5.7939238e-17 -4632.0136 -0.00015032304
1000 298.82235 -19950.274 13208.089 5628.7024 -0.0056205182
2000 300.43933 -28320.212 11980.296 -3902.0877 -0.0045324757
3000 300.4263 -36049.855 11338.405 -12274.161 -0.0018833539
4000 299.13368 -43471.21 11926.882 -19160.553 -0.00043030866
5000 293.77858 -50083.893 12334.927 -25586.884 -0.0015653738
6000 296.4851 -56330.135 12325.63 -31730.376 -0.0012795986
7000 298.20879 -62120.359 12582.297 -37192.574 -0.0013845796
8000 299.45547 -67881.692 13058.926 -42425.669 -0.00021100885
9000 301.82622 -73333.698 13598.257 -47240.197 -0.0006009197
10000 307.16873 -78292.306 13818.929 -51756.96 -0.0005609903
Loop time of 4.0316 on 1 procs for 10000 steps with 2000 atoms
Performance: 2143.072 ns/day, 0.011 hours/ns, 2480.408 timesteps/s
99.9% CPU use with 1 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.5955 | 2.5955 | 2.5955 | 0.0 | 64.38
Bond | 1.1516 | 1.1516 | 1.1516 | 0.0 | 28.57
Neigh | 0.001163 | 0.001163 | 0.001163 | 0.0 | 0.03
Comm | 0.0019577 | 0.0019577 | 0.0019577 | 0.0 | 0.05
Output | 0.020854 | 0.020854 | 0.020854 | 0.0 | 0.52
Modify | 0.21637 | 0.21637 | 0.21637 | 0.0 | 5.37
Other | | 0.04409 | | | 1.09
Nlocal: 2000 ave 2000 max 2000 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 0 ave 0 max 0 min
Histogram: 1 0 0 0 0 0 0 0 0 0
FullNghs: 13320 ave 13320 max 13320 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 13320
Ave neighs/atom = 6.66
Ave special neighs/atom = 5.952
Neighbor list builds = 1
Dangerous builds = 0
Total wall time: 0:00:05

126
examples/PACKAGES/mesont/log.9Jan20.cnt.g++.4
View File

@ -1,126 +0,0 @@
LAMMPS (09 Jan 2020)
#Initialisation
units nano
dimension 3
boundary p p p
atom_style full
comm_modify cutoff 11.0
neighbor 7.80 bin
newton on
#Read data
read_data cnt.data
orthogonal box = (0 0 0) to (600 600 60)
2 by 2 by 1 MPI processor grid
reading atoms ...
500 atoms
scanning bonds ...
1 = max bonds/atom
scanning angles ...
1 = max angles/atom
reading bonds ...
498 bonds
reading angles ...
496 angles
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.000354767 secs
read_data CPU = 0.00286365 secs
replicate 1 2 2
orthogonal box = (0 0 0) to (600 1200 120)
1 by 4 by 1 MPI processor grid
2000 atoms
1992 bonds
1984 angles
2 = max # of 1-2 neighbors
2 = max # of 1-3 neighbors
4 = max # of 1-4 neighbors
6 = max # of special neighbors
special bonds CPU = 0.00019598 secs
replicate CPU = 0.00055337 secs
#Force field
bond_style harmonic
bond_coeff 1 268896.77 2.0
angle_style harmonic
angle_coeff 1 46562.17 180.0
pair_style mesocnt
pair_coeff * * 10_10.cnt
Reading potential file 10_10.cnt with DATE: 2020-01-13
#Output
thermo 1000
dump xyz all xyz 1000 cnt.xyz
#Simulation setup
timestep 1.0e-05
#Nose-Hoover thermostat
fix nvt all nvt temp 300 300 0.001
run 10000
WARNING: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions (src/pair.cpp:226)
Neighbor list info ...
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 10.177
ghost atom cutoff = 11
binsize = 5.0885, bins = 118 236 24
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair mesocnt, perpetual
attributes: full, newton on
pair build: full/bin
stencil: full/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.725 | 2.725 | 2.725 Mbytes
Step Temp E_pair E_mol TotEng Press
0 0 -4632.0136 5.7939238e-17 -4632.0136 -0.00015032304
1000 298.82235 -19950.274 13208.089 5628.7024 -0.0056205182
2000 300.43861 -28320.205 11980.287 -3902.1202 -0.0045324738
3000 300.41076 -36049.308 11339.149 -12273.513 -0.0018848513
4000 299.13326 -43471.424 11927.668 -19159.998 -0.00042845101
5000 293.78857 -50083.216 12333.969 -25586.752 -0.0015664633
6000 296.45482 -56329.621 12326.419 -31730.328 -0.0012773686
7000 298.19097 -62119.086 12581.4 -37192.937 -0.0013862831
8000 299.46424 -67880.989 13057.62 -42425.908 -0.00020874264
9000 301.80677 -73332.208 13597.237 -47240.532 -0.00060074773
10000 307.17104 -78292.912 13818.889 -51757.51 -0.00056148282
Loop time of 1.23665 on 4 procs for 10000 steps with 2000 atoms
Performance: 6986.607 ns/day, 0.003 hours/ns, 8086.351 timesteps/s
96.1% CPU use with 4 MPI tasks x no OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.66321 | 0.68439 | 0.71413 | 2.5 | 55.34
Bond | 0.28561 | 0.29434 | 0.30976 | 1.7 | 23.80
Neigh | 0.00043321 | 0.00043637 | 0.00043917 | 0.0 | 0.04
Comm | 0.026656 | 0.05346 | 0.097228 | 12.7 | 4.32
Output | 0.0070224 | 0.0073031 | 0.0081415 | 0.6 | 0.59
Modify | 0.12769 | 0.15394 | 0.18743 | 6.5 | 12.45
Other | | 0.04279 | | | 3.46
Nlocal: 500 ave 504 max 496 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost: 22 ave 24 max 20 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 0 ave 0 max 0 min
Histogram: 4 0 0 0 0 0 0 0 0 0
FullNghs: 3330 ave 3368 max 3292 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Total # of neighbors = 13320
Ave neighs/atom = 6.66
Ave special neighs/atom = 5.952
Neighbor list builds = 1
Dangerous builds = 0
Total wall time: 0:00:02

11
src/.gitignore vendored
View File

@ -268,6 +268,15 @@
/pair_mesont_tpm.cpp
/pair_mesont_tpm.h
/bond_mesocnt.cpp
/bond_mesocnt.h
/angle_mesocnt.cpp
/angle_mesocnt.h
/pair_mesocnt.cpp
/pair_mesocnt.h
/pair_mesocnt_viscous.cpp
/pair_mesocnt_viscous.h
/atom_vec_bpm_sphere.cpp
/atom_vec_bpm_sphere.h
/bond_bpm.cpp
@ -1240,8 +1249,6 @@
/pair_meam_spline.h
/pair_meam_sw_spline.cpp
/pair_meam_sw_spline.h
/pair_mesocnt.cpp
/pair_mesocnt.h
/pair_mm3_switch3_coulgauss_long.cpp
/pair_mm3_switch3_coulgauss_long.h
/pair_morse_soft.cpp

397
src/MESONT/angle_mesocnt.cpp Normal file
View File

@ -0,0 +1,397 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Philipp Kloza (University of Cambridge)
pak37@cam.ac.uk
------------------------------------------------------------------------- */
#include "angle_mesocnt.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
using namespace LAMMPS_NS;
using MathConst::DEG2RAD;
using MathConst::MY_2PI;
using MathConst::MY_PI;
using MathConst::RAD2DEG;
static constexpr double SMALL = 0.001;
static constexpr double A_CC = 1.421;
/* ---------------------------------------------------------------------- */
AngleMesoCNT::AngleMesoCNT(LAMMPS *_lmp) : Angle(_lmp)
{
buckling = nullptr;
kh = nullptr;
kb = nullptr;
thetab = nullptr;
}
/* ---------------------------------------------------------------------- */
AngleMesoCNT::~AngleMesoCNT()
{
if (allocated && !copymode) {
memory->destroy(setflag);
memory->destroy(buckling);
memory->destroy(kh);
memory->destroy(kb);
memory->destroy(thetab);
memory->destroy(theta0);
}
}
/* ---------------------------------------------------------------------- */
void AngleMesoCNT::compute(int eflag, int vflag)
{
int i1, i2, i3, n, type;
double delx1, dely1, delz1, delx2, dely2, delz2;
double eangle, f1[3], f3[3];
double dtheta, tk;
double rsq1, rsq2, r1, r2, c, s, a, a11, a12, a22;
eangle = 0.0;
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
int **anglelist = neighbor->anglelist;
int nanglelist = neighbor->nanglelist;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
int flag, cols;
int index = atom->find_custom("buckled", flag, cols);
int *buckled = atom->ivector[index];
for (n = 0; n < nanglelist; n++) {
i1 = anglelist[n][0];
i2 = anglelist[n][1];
i3 = anglelist[n][2];
type = anglelist[n][3];
// 1st bond
delx1 = x[i1][0] - x[i2][0];
dely1 = x[i1][1] - x[i2][1];
delz1 = x[i1][2] - x[i2][2];
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
r1 = sqrt(rsq1);
// 2nd bond
delx2 = x[i3][0] - x[i2][0];
dely2 = x[i3][1] - x[i2][1];
delz2 = x[i3][2] - x[i2][2];
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
r2 = sqrt(rsq2);
// angle (cos and sin)
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
c /= r1 * r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
s = sqrt(1.0 - c * c);
if (s < SMALL) s = SMALL;
s = 1.0 / s;
// force & energy
dtheta = acos(c) - theta0[type];
// harmonic bending
if (!buckling[type] || fabs(dtheta) < thetab[type]) {
tk = kh[type] * dtheta;
if (eflag) eangle = tk * dtheta;
a = -2.0 * tk * s;
buckled[i2] = 0;
}
// bending buckling
else {
if (eflag)
eangle = kb[type] * fabs(dtheta) + thetab[type] * (kh[type] * thetab[type] - kb[type]);
a = kb[type] * s;
buckled[i2] = 1;
}
a11 = a * c / rsq1;
a12 = -a / (r1 * r2);
a22 = a * c / rsq2;
f1[0] = a11 * delx1 + a12 * delx2;
f1[1] = a11 * dely1 + a12 * dely2;
f1[2] = a11 * delz1 + a12 * delz2;
f3[0] = a22 * delx2 + a12 * delx1;
f3[1] = a22 * dely2 + a12 * dely1;
f3[2] = a22 * delz2 + a12 * delz1;
// apply force to each of 3 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += f1[0];
f[i1][1] += f1[1];
f[i1][2] += f1[2];
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= f1[0] + f3[0];
f[i2][1] -= f1[1] + f3[1];
f[i2][2] -= f1[2] + f3[2];
}
if (newton_bond || i3 < nlocal) {
f[i3][0] += f3[0];
f[i3][1] += f3[1];
f[i3][2] += f3[2];
}
if (evflag)
ev_tally(i1, i2, i3, nlocal, newton_bond, eangle, f1, f3, delx1, dely1, delz1, delx2, dely2,
delz2);
}
}
/* ---------------------------------------------------------------------- */
void AngleMesoCNT::allocate()
{
allocated = 1;
const int np1 = atom->nangletypes + 1;
memory->create(buckling, np1, "angle:buckling");
memory->create(kh, np1, "angle:kh");
memory->create(kb, np1, "angle:kb");
memory->create(thetab, np1, "angle:thetab");
memory->create(theta0, np1, "angle:theta0");
memory->create(setflag, np1, "angle:setflag");
for (int i = 1; i < np1; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
void AngleMesoCNT::coeff(int narg, char **arg)
{
if (narg < 1) utils::missing_cmd_args(FLERR, "angle_coeff", error);
int buckling_one;
if (strcmp(arg[1], "buckling") == 0)
buckling_one = 1;
else if (strcmp(arg[1], "harmonic") == 0)
buckling_one = 0;
else
error->all(FLERR,
"Unknown first argument {} for angle coefficients, must be 'buckling' or 'harmonic'",
arg[1]);
// units, eV to energy unit conversion
double ang = force->angstrom;
double eunit;
if (strcmp(update->unit_style, "real") == 0)
eunit = 23.06054966;
else if (strcmp(update->unit_style, "metal") == 0)
eunit = 1.0;
else if (strcmp(update->unit_style, "si") == 0)
eunit = 1.6021765e-19;
else if (strcmp(update->unit_style, "cgs") == 0)
eunit = 1.6021765e-12;
else if (strcmp(update->unit_style, "electron") == 0)
eunit = 3.674932248e-2;
else if (strcmp(update->unit_style, "micro") == 0)
eunit = 1.6021765e-4;
else if (strcmp(update->unit_style, "nano") == 0)
eunit = 1.6021765e2;
else
error->all(FLERR, "Angle style mesocnt does not support {} units", update->unit_style);
// set parameters
double kh_one, kb_one, thetab_one;
if (strcmp(arg[2], "custom") == 0) {
if (buckling_one) {
if (narg != 6) error->all(FLERR, "Incorrect number of args for 'custom' angle coefficients");
kb_one = utils::numeric(FLERR, arg[4], false, lmp);
thetab_one = utils::numeric(FLERR, arg[5], false, lmp);
} else if (narg != 4)
error->all(FLERR, "Incorrect number of args for 'custom' angle coefficients");
kh_one = utils::numeric(FLERR, arg[3], false, lmp);
} else if (strcmp(arg[2], "C") == 0) {
if (narg != 6)
error->all(FLERR, "Incorrect number of args for 'C' preset in angle coefficients");
int n = utils::inumeric(FLERR, arg[3], false, lmp);
int m = utils::inumeric(FLERR, arg[4], false, lmp);
double l = utils::numeric(FLERR, arg[5], false, lmp);
double r_ang = sqrt(3.0 * (n * n + n * m + m * m)) * A_CC / MY_2PI;
// empirical parameters
double k = 63.80 * pow(r_ang, 2.93) * eunit * ang;
kh_one = 0.5 * k / l;
if (buckling_one) {
kb_one = 0.7 * k / (275.0 * ang);
thetab_one = 180.0 / MY_PI * atan(l / (275.0 * ang));
}
} else
error->all(FLERR, "Unknown {} preset in angle coefficients", arg[2]);
// set safe default values for buckling parameters if buckling is disabled
if (!buckling_one) {
kb_one = 0.0;
thetab_one = 180.0;
}
if (!allocated) allocate();
int ilo, ihi;
utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
// convert thetab from degrees to radians
int count = 0;
for (int i = ilo; i <= ihi; i++) {
buckling[i] = buckling_one;
kh[i] = kh_one;
kb[i] = kb_one;
thetab[i] = DEG2RAD * thetab_one;
theta0[i] = DEG2RAD * 180.0;
setflag[i] = 1;
count++;
}
if (count == 0) error->all(FLERR, "Invalid angle type {}", arg[0]);
}
/* ---------------------------------------------------------------------- */
void AngleMesoCNT::init_style()
{
std::string id_fix = "angle_mesocnt_buckled";
if (!modify->get_fix_by_id(id_fix))
modify->add_fix(id_fix + " all property/atom i_buckled ghost yes");
}
/* ---------------------------------------------------------------------- */
double AngleMesoCNT::equilibrium_angle(int /*i*/)
{
return 180.0;
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void AngleMesoCNT::write_restart(FILE *fp)
{
fwrite(&buckling[1], sizeof(int), atom->nangletypes, fp);
fwrite(&kh[1], sizeof(double), atom->nangletypes, fp);
fwrite(&kb[1], sizeof(double), atom->nangletypes, fp);
fwrite(&thetab[1], sizeof(double), atom->nangletypes, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
void AngleMesoCNT::read_restart(FILE *fp)
{
allocate();
if (comm->me == 0) {
utils::sfread(FLERR, &buckling[1], sizeof(int), atom->nangletypes, fp, nullptr, error);
utils::sfread(FLERR, &kh[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
utils::sfread(FLERR, &kb[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
utils::sfread(FLERR, &thetab[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
}
MPI_Bcast(&buckling[1], atom->nangletypes, MPI_INT, 0, world);
MPI_Bcast(&kh[1], atom->nangletypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&kb[1], atom->nangletypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&thetab[1], atom->nangletypes, MPI_DOUBLE, 0, world);
for (int i = 1; i <= atom->nangletypes; i++) {
theta0[i] = 180.0;
setflag[i] = 1;
}
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void AngleMesoCNT::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nangletypes; i++)
fprintf(fp, "%d %d %g %g %g\n", i, buckling[i], kh[i], kb[i], RAD2DEG * thetab[i]);
}
/* ---------------------------------------------------------------------- */
double AngleMesoCNT::single(int type, int i1, int i2, int i3)
{
double **x = atom->x;
double delx1 = x[i1][0] - x[i2][0];
double dely1 = x[i1][1] - x[i2][1];
double delz1 = x[i1][2] - x[i2][2];
domain->minimum_image(delx1, dely1, delz1);
double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1);
double delx2 = x[i3][0] - x[i2][0];
double dely2 = x[i3][1] - x[i2][1];
double delz2 = x[i3][2] - x[i2][2];
domain->minimum_image(delx2, dely2, delz2);
double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2);
double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
c /= r1 * r2;
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
double dtheta = acos(c) - theta0[type];
// harmonic bending
if (!buckling[type] || dtheta < thetab[type]) {
double tk = kh[type] * dtheta;
return tk * dtheta;
}
// bending buckling
else
return kb[type] * dtheta + thetab[type] * (kh[type] * thetab[type] - kb[type]);
}

55
src/MESONT/angle_mesocnt.h Normal file
View File

@ -0,0 +1,55 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Philipp Kloza (University of Cambridge)
pak37@cam.ac.uk
------------------------------------------------------------------------- */
#ifdef ANGLE_CLASS
// clang-format off
AngleStyle(mesocnt, AngleMesoCNT);
// clang-format on
#else
#ifndef LMP_ANGLE_MESOCNT_H
#define LMP_ANGLE_MESOCNT_H
#include "angle.h"
namespace LAMMPS_NS {
class AngleMesoCNT : public Angle {
public:
AngleMesoCNT(class LAMMPS *);
~AngleMesoCNT() override;
void compute(int, int) override;
void coeff(int, char **) override;
void init_style() override;
double equilibrium_angle(int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_data(FILE *) override;
double single(int, int, int, int) override;
protected:
int *buckling;
double *kh, *kb, *thetab, *theta0;
virtual void allocate();
};
} // namespace LAMMPS_NS
#endif
#endif

118
src/MESONT/bond_mesocnt.cpp Normal file
View File

@ -0,0 +1,118 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Philipp Kloza (University of Cambridge)
pak37@cam.ac.uk
------------------------------------------------------------------------- */
#include "bond_mesocnt.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using MathConst::MY_2PI;
static constexpr double A_CC = 1.421;
/* ---------------------------------------------------------------------- */
BondMesoCNT::BondMesoCNT(LAMMPS *_lmp) : BondHarmonic(_lmp)
{
born_matrix_enable = 1;
}
/* ---------------------------------------------------------------------- */
BondMesoCNT::~BondMesoCNT()
{
if (allocated && !copymode) {
memory->destroy(setflag);
memory->destroy(k);
memory->destroy(r0);
}
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
void BondMesoCNT::coeff(int narg, char **arg)
{
if (narg < 1) utils::missing_cmd_args(FLERR, "bond_coeff", error);
// units, eV to energy unit conversion
double ang = force->angstrom;
double eunit;
if (strcmp(update->unit_style, "real") == 0)
eunit = 23.06054966;
else if (strcmp(update->unit_style, "metal") == 0)
eunit = 1.0;
else if (strcmp(update->unit_style, "si") == 0)
eunit = 1.6021765e-19;
else if (strcmp(update->unit_style, "cgs") == 0)
eunit = 1.6021765e-12;
else if (strcmp(update->unit_style, "electron") == 0)
eunit = 3.674932248e-2;
else if (strcmp(update->unit_style, "micro") == 0)
eunit = 1.6021765e-4;
else if (strcmp(update->unit_style, "nano") == 0)
eunit = 1.6021765e2;
else
error->all(FLERR, "Bond style mesocnt does not support {} units", update->unit_style);
// set parameters
double k_one, r0_one;
if (strcmp(arg[1], "custom") == 0) {
if (narg != 4) error->all(FLERR, "Incorrect number of args for 'custom' bond coefficients");
k_one = utils::numeric(FLERR, arg[2], false, lmp);
r0_one = utils::numeric(FLERR, arg[3], false, lmp);
} else if (strcmp(arg[1], "C") == 0) {
if (narg != 5)
error->all(FLERR, "Incorrect number of args for 'C' preset in bond coefficients");
int n = utils::inumeric(FLERR, arg[2], false, lmp);
int m = utils::inumeric(FLERR, arg[3], false, lmp);
r0_one = utils::numeric(FLERR, arg[4], false, lmp);
double r_ang = sqrt(3.0 * (n * n + n * m + m * m)) * A_CC / MY_2PI;
k_one = 0.5 * (86.64 + 100.56 * r_ang) * eunit / (ang * r0_one);
} else
error->all(FLERR, "Unknown {} preset in bond coefficients", arg[1]);
if (!allocated) allocate();
int ilo, ihi;
utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
k[i] = k_one;
r0[i] = r0_one;
setflag[i] = 1;
count++;
}
if (count == 0) error->all(FLERR, "Invalid bond type {}", arg[0]);
}

42
src/MESONT/bond_mesocnt.h Normal file
View File

@ -0,0 +1,42 @@
/* -*- 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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Philipp Kloza (University of Cambridge)
pak37@cam.ac.uk
------------------------------------------------------------------------- */
#ifdef BOND_CLASS
// clang-format off
BondStyle(mesocnt, BondMesoCNT);
// clang-format on
#else
#ifndef LMP_BOND_MESOCNT_H
#define LMP_BOND_MESOCNT_H
#include "bond_harmonic.h"
namespace LAMMPS_NS {
class BondMesoCNT : public BondHarmonic {
public:
BondMesoCNT(class LAMMPS *);
~BondMesoCNT() override;
void coeff(int, char **) override;
};
} // namespace LAMMPS_NS
#endif
#endif

1410
src/MESONT/pair_mesocnt.cpp
View File

File diff suppressed because it is too large Load Diff

50
src/MESONT/pair_mesocnt.h
View File

@ -11,8 +11,15 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Philipp Kloza (University of Cambridge)
pak37@cam.ac.uk
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(mesocnt, PairMesoCNT);
// clang-format on
#else
#ifndef LMP_PAIR_MESOCNT_H
@ -32,15 +39,21 @@ class PairMesoCNT : public Pair {
void init_style() override;
double init_one(int, int) override;
int pack_forward_comm(int, int *, double *, int, int *) override;
void unpack_forward_comm(int, int, double *) override;
protected:
int nend_types;
int uinf_points, gamma_points, phi_points, usemi_points;
int *reduced_nlist, *numchainlist;
int **reduced_neighlist, **nchainlist, **endlist;
int *end_types, *reduced_nlist, *numchainlist, *selfid;
int **reduced_neighlist, **nchainlist, **endlist, **selfpos;
int ***chainlist;
bool segment_flag, neigh_flag;
double ang, ang_inv, eunit, funit;
double delta1, delta2;
double r, rsq, d, rc, rcsq, rc0, cutoff, cutoffsq;
double r, rsq, d, rc, rcsq, rc0, cutoff, cutoffsq, neigh_cutoff;
double r_ang, rsq_ang, d_ang, rc_ang, rcsq_ang, cutoff_ang, cutoffsq_ang;
double sig, sig_ang, comega, ctheta;
double hstart_uinf, hstart_gamma, hstart_phi, psistart_phi, hstart_usemi, xistart_usemi;
@ -50,18 +63,24 @@ class PairMesoCNT : public Pair {
double *param, *w, *wnode;
double **dq_w;
double ***q1_dq_w, ***q2_dq_w;
double *gl_nodes_finf, *gl_nodes_fsemi;
double *gl_weights_finf, *gl_weights_fsemi;
double *uinf_data, *gamma_data, **phi_data, **usemi_data;
double **uinf_coeff, **gamma_coeff, ****phi_coeff, ****usemi_coeff;
double **flocal, **fglobal, **basis;
bool match_end(int);
int count_chains(int *, int);
void allocate();
void bond_neigh();
void bond_neigh_id();
void bond_neigh_topo();
void neigh_common(int, int, int &, int *);
void chain_lengths(int *, int, int *);
void chain_split(int *, int, int *, int **, int *);
void sort(int *, int);
void read_file(const char *);
void read_data(PotentialFileReader &, double *, double &, double &, int);
void read_data(PotentialFileReader &, double **, double &, double &, double &, double &, int);
@ -69,22 +88,29 @@ class PairMesoCNT : public Pair {
void spline_coeff(double *, double **, double, int);
void spline_coeff(double **, double ****, double, double, int);
double spline(double, double, double, double **, int);
double dspline(double, double, double, double **, int);
double spline(double, double, double, double, double, double, double ****, int);
double dxspline(double, double, double, double, double, double, double ****, int);
double dyspline(double, double, double, double, double, double, double ****, int);
void geometry(const double *, const double *, const double *, const double *, const double *,
double *, double *, double *, double **);
void weight(const double *, const double *, const double *, const double *, double &, double *,
double *, double *, double *);
void finf(const double *, double &, double **);
void fsemi(const double *, double &, double &, double **);
// Legendre-Gauss integration
double legendre(int, double);
void gl_init_nodes(int, double *);
void gl_init_weights(int, double *, double *);
// inlined functions for efficiency
inline void weight(const double *, const double *, const double *, const double *, double &,
double *, double *, double *, double *);
inline double spline(double, double, double, double **, int);
inline double dspline(double, double, double, double **, int);
inline double spline(double, double, double, double, double, double, double ****, int);
inline double dxspline(double, double, double, double, double, double, double ****, int);
inline double dyspline(double, double, double, double, double, double, double ****, int);
inline double heaviside(double x)
{
if (x > 0)

871
src/MESONT/pair_mesocnt_viscous.cpp Normal file
View File

@ -0,0 +1,871 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Philipp Kloza (University of Cambridge)
pak37@cam.ac.uk
------------------------------------------------------------------------- */
#include "pair_mesocnt_viscous.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace MathExtra;
using MathConst::MY_PI;
#define SELF_CUTOFF 3
#define RHOMIN 10.0
#define QUAD_FINF 129
#define QUAD_FSEMI 10
/* ---------------------------------------------------------------------- */
void PairMesoCNTViscous::compute(int eflag, int vflag)
{
ev_init(eflag, vflag);
int i, j, k, i1, i2, j1, j2;
int endflag, endindex;
int clen, numchain;
int *end, *nchain;
int **chain;
double fend, lp, scale, sumw, sumw_inv;
double evdwl, evdwl_chain;
double vtot, fvisc_tot;
double *r1, *r2, *q1, *q2, *qe;
double *vr1, *vr2, *vq1, *vq2;
double vr[3], vp1[3], vp2[3], vp[3], vrel[3], fvisc[3];
double ftotal[3], ftorque[3], torque[3], delr1[3], delr2[3], delqe[3];
double t1[3], t2[3], t3[3];
double dr1_sumw[3], dr2_sumw[3];
double dr1_w[3], dr2_w[3], dq1_w[3], dq2_w[3];
double fgrad_r1_p1[3], fgrad_r1_p2[3], fgrad_r2_p1[3], fgrad_r2_p2[3];
double fgrad_q_p1[3], fgrad_q_p2[3];
double q1_dr1_w[3][3], q1_dr2_w[3][3], q2_dr1_w[3][3], q2_dr2_w[3][3];
double dr1_p1[3][3], dr1_p2[3][3], dr2_p1[3][3], dr2_p2[3][3];
double dq_p1[3][3], dq_p2[3][3];
double temp[3][3];
evdwl = 0.0;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
int flag, cols;
int buckled_index = atom->find_custom("buckled", flag, cols);
// update bond neighbor list when necessary
if (update->ntimestep == neighbor->lastcall) {
if (neigh_flag)
bond_neigh_topo();
else
bond_neigh_id();
}
// iterate over all bonds
for (i = 0; i < nbondlist; i++) {
i1 = bondlist[i][0];
i2 = bondlist[i][1];
r1 = x[i1];
r2 = x[i2];
vr1 = v[i1];
vr2 = v[i2];
add3(vr1, vr2, vr);
scale3(0.5, vr);
numchain = numchainlist[i];
end = endlist[i];
nchain = nchainlist[i];
chain = chainlist[i];
// iterate over all neighbouring chains
for (j = 0; j < numchain; j++) {
clen = nchain[j];
if (clen < 2) continue;
// check if segment-segment interactions are necessary
endflag = end[j];
int buckled = 0;
if (buckled_index != -1) {
for (k = 0; k < clen; k++) {
if (atom->ivector[buckled_index][chain[j][k]]) {
buckled = 1;
break;
}
}
}
if (segment_flag || buckled || j == selfid[i] || endflag == 3) {
zero3(vp1);
zero3(vp2);
sumw = 0.0;
for (k = 0; k < clen - 1; k++) {
// segment-segment interaction
// exclude SELF_CUTOFF neighbors in self-chain
if (j == selfid[i]) {
int min11 = abs(k - selfpos[i][0]);
int min12 = abs(k - selfpos[i][1]);
int min21 = abs(k + 1 - selfpos[i][0]);
int min22 = abs(k + 1 - selfpos[i][1]);
int min = min11;
if (min12 < min) min = min12;
if (min21 < min) min = min21;
if (min22 < min) min = min22;
if (min < SELF_CUTOFF) continue;
}
j1 = chain[j][k];
j2 = chain[j][k + 1];
j1 &= NEIGHMASK;
j2 &= NEIGHMASK;
q1 = x[j1];
q2 = x[j2];
vq1 = v[j1];
vq2 = v[j2];
// weighted velocity for friction
w[k] = weight(r1, r2, q1, q2);
if (w[k] == 0.0) continue;
sumw += w[k];
scaleadd3(w[k], vq1, vp1, vp1);
scaleadd3(w[k], vq2, vp2, vp2);
// check if orientation of segment needs to be flipped to prevent overlap
geometry(r1, r2, q1, q2, q1, p, m, param, basis);
if (param[0] > cutoff) continue;
double calpha = cos(param[1]);
double salpha = sin(param[1]);
double hsq = param[0] * param[0];
double ceta = calpha * param[4];
double seta = salpha * param[4];
double dsq1 = hsq + seta * seta;
if (ceta < param[2]) {
double dceta = ceta - param[2];
dsq1 += dceta * dceta;
} else if (ceta > param[3]) {
double dceta = ceta - param[3];
dsq1 += dceta * dceta;
}
ceta = calpha * param[5];
seta = salpha * param[5];
double dsq2 = hsq + seta * seta;
if (ceta < param[2]) {
double dceta = ceta - param[2];
dsq2 += dceta * dceta;
} else if (ceta > param[3]) {
double dceta = ceta - param[3];
dsq2 += dceta * dceta;
}
if (dsq1 > cutoffsq && dsq2 > cutoffsq) continue;
int jj1, jj2;
// do flip if necessary
if (dsq1 < dsq2) {
jj1 = j1;
jj2 = j2;
} else {
if (param[1] > MY_PI)
param[1] -= MY_PI;
else
param[1] += MY_PI;
double eta2 = -param[5];
param[5] = -param[4];
param[4] = eta2;
param[6] = eta2;
negate3(m);
jj1 = j2;
jj2 = j1;
}
// first force contribution
fsemi(param, evdwl, fend, flocal);
if (evdwl == 0.0) continue;
// transform to global coordinate system
matvec(basis[0], basis[1], basis[2], flocal[0], fglobal[0]);
matvec(basis[0], basis[1], basis[2], flocal[1], fglobal[1]);
// forces acting on segment 2
add3(fglobal[0], fglobal[1], ftotal);
scaleadd3(fend, m, ftotal, ftotal);
scale3(-0.5, ftotal);
sub3(r1, p, delr1);
sub3(r2, p, delr2);
cross3(delr1, fglobal[0], t1);
cross3(delr2, fglobal[1], t2);
add3(t1, t2, torque);
cross3(torque, m, ftorque);
lp = param[5] - param[4];
scale3(1.0 / lp, ftorque);
add3(ftotal, ftorque, fglobal[2]);
sub3(ftotal, ftorque, fglobal[3]);
// add forces to nodes
scaleadd3(0.5, fglobal[0], f[i1], f[i1]);
scaleadd3(0.5, fglobal[1], f[i2], f[i2]);
scaleadd3(0.5, fglobal[2], f[jj1], f[jj1]);
scaleadd3(0.5, fglobal[3], f[jj2], f[jj2]);
scaleadd3(0.5 * fend, m, f[jj1], f[jj1]);
// add energy
if (eflag_either) {
evdwl = 0.5 * evdwl;
double evdwl_atom = 0.25 * evdwl;
if (eflag_global) eng_vdwl += evdwl;
if (eflag_atom) {
eatom[i1] += evdwl_atom;
eatom[i2] += evdwl_atom;
eatom[jj1] += evdwl_atom;
eatom[jj2] += evdwl_atom;
}
}
// second force contribution
param[6] += lp;
fsemi(param, evdwl, fend, flocal);
if (evdwl == 0.0) continue;
// transform to global coordinate system
matvec(basis[0], basis[1], basis[2], flocal[0], fglobal[0]);
matvec(basis[0], basis[1], basis[2], flocal[1], fglobal[1]);
// forces acting on segment 2
add3(fglobal[0], fglobal[1], ftotal);
scaleadd3(fend, m, ftotal, ftotal);
scale3(-0.5, ftotal);
scaleadd3(lp, m, p, p);
sub3(r1, p, delr1);
sub3(r2, p, delr2);
cross3(delr1, fglobal[0], t1);
cross3(delr2, fglobal[1], t2);
add3(t1, t2, torque);
cross3(torque, m, ftorque);
scale3(1.0 / lp, ftorque);
add3(ftotal, ftorque, fglobal[2]);
sub3(ftotal, ftorque, fglobal[3]);
// add forces to nodes
scaleadd3(-0.5, fglobal[0], f[i1], f[i1]);
scaleadd3(-0.5, fglobal[1], f[i2], f[i2]);
scaleadd3(-0.5, fglobal[2], f[jj2], f[jj2]);
scaleadd3(0.5, fglobal[3], f[jj1], f[jj1]);
sub3(f[jj2], fglobal[3], f[jj2]);
scaleadd3(-0.5 * fend, m, f[jj2], f[jj2]);
// add energy
if (eflag_either) {
evdwl = -0.5 * evdwl;
double evdwl_atom = 0.25 * evdwl;
if (eflag_global) eng_vdwl += evdwl;
if (eflag_atom) {
eatom[i1] += evdwl_atom;
eatom[i2] += evdwl_atom;
eatom[jj1] += evdwl_atom;
eatom[jj2] += evdwl_atom;
}
}
}
if (sumw == 0.0) continue;
sumw_inv = 1.0 / sumw;
// mean chain velocity and relative velocity
add3(vp1, vp2, vp);
scale3(0.5 * sumw_inv, vp);
sub3(vp, vr, vrel);
vtot = dot3(vrel, basis[2]);
// friction forces
if (vtot == 0.0)
fvisc_tot = 0.0;
else {
fvisc_tot = fvisc_max / (1.0 + exp(-kvisc * (fabs(vtot) - vvisc))) - fvisc_shift;
fvisc_tot *= 0.25 * (1 - s(sumw)) * (param[3] - param[2]);
if (vtot < 0) fvisc_tot = -fvisc_tot;
}
scale3(fvisc_tot, basis[2], fvisc);
// add friction forces to current segment
add3(fvisc, f[i1], f[i1]);
add3(fvisc, f[i2], f[i2]);
// add friction forces to neighbor chain
for (k = 0; k < clen - 1; k++) {
if (w[k] == 0.0) continue;
j1 = chain[j][k];
j2 = chain[j][k + 1];
j1 &= NEIGHMASK;
j2 &= NEIGHMASK;
scale = w[k] * sumw_inv;
scaleadd3(-scale, fvisc, f[j1], f[j1]);
scaleadd3(-scale, fvisc, f[j2], f[j2]);
}
} else {
// segment-chain interaction
// assign end position
endflag = end[j];
if (endflag == 1) {
endindex = chain[j][0];
qe = x[endindex];
} else if (endflag == 2) {
endindex = chain[j][clen - 1];
qe = x[endindex];
}
// compute substitute straight (semi-)infinite CNT
zero3(p1);
zero3(p2);
zero3(dr1_sumw);
zero3(dr2_sumw);
zero3(vp1);
zero3(vp2);
zeromat3(q1_dr1_w);
zeromat3(q2_dr1_w);
zeromat3(q1_dr2_w);
zeromat3(q2_dr2_w);
for (k = 0; k < clen; k++) {
wnode[k] = 0.0;
zero3(dq_w[k]);
zeromat3(q1_dq_w[k]);
zeromat3(q2_dq_w[k]);
}
sumw = 0.0;
for (k = 0; k < clen - 1; k++) {
j1 = chain[j][k];
j2 = chain[j][k + 1];
j1 &= NEIGHMASK;
j2 &= NEIGHMASK;
q1 = x[j1];
q2 = x[j2];
vq1 = v[j1];
vq2 = v[j2];
weight(r1, r2, q1, q2, w[k], dr1_w, dr2_w, dq1_w, dq2_w);
if (w[k] == 0.0) {
if (endflag == 1 && k == 0)
endflag = 0;
else if (endflag == 2 && k == clen - 2)
endflag = 0;
continue;
}
sumw += w[k];
wnode[k] += w[k];
wnode[k + 1] += w[k];
scaleadd3(w[k], q1, p1, p1);
scaleadd3(w[k], q2, p2, p2);
// weighted velocity for friction
scaleadd3(w[k], vq1, vp1, vp1);
scaleadd3(w[k], vq2, vp2, vp2);
// weight gradient terms
add3(dr1_w, dr1_sumw, dr1_sumw);
add3(dr2_w, dr2_sumw, dr2_sumw);
outer3(q1, dr1_w, temp);
plus3(temp, q1_dr1_w, q1_dr1_w);
outer3(q2, dr1_w, temp);
plus3(temp, q2_dr1_w, q2_dr1_w);
outer3(q1, dr2_w, temp);
plus3(temp, q1_dr2_w, q1_dr2_w);
outer3(q2, dr2_w, temp);
plus3(temp, q2_dr2_w, q2_dr2_w);
add3(dq1_w, dq_w[k], dq_w[k]);
add3(dq2_w, dq_w[k + 1], dq_w[k + 1]);
outer3(q1, dq1_w, temp);
plus3(temp, q1_dq_w[k], q1_dq_w[k]);
outer3(q1, dq2_w, temp);
plus3(temp, q1_dq_w[k + 1], q1_dq_w[k + 1]);
outer3(q2, dq1_w, temp);
plus3(temp, q2_dq_w[k], q2_dq_w[k]);
outer3(q2, dq2_w, temp);
plus3(temp, q2_dq_w[k + 1], q2_dq_w[k + 1]);
}
if (sumw == 0.0) continue;
sumw_inv = 1.0 / sumw;
scale3(sumw_inv, p1);
scale3(sumw_inv, p2);
// compute geometry and forces
if (endflag == 0) {
// infinite CNT case
geometry(r1, r2, p1, p2, nullptr, p, m, param, basis);
if (param[0] > cutoff) continue;
if (!(param[2] < 0 && param[3] > 0)) {
double salpha = sin(param[1]);
double sxi1 = salpha * param[2];
double sxi2 = salpha * param[3];
double hsq = param[0] * param[0];
if (sxi1 * sxi1 + hsq > cutoffsq && sxi2 * sxi2 + hsq > cutoffsq) continue;
}
finf(param, evdwl, flocal);
} else {
// semi-infinite CNT case
// endflag == 1: CNT end at start of chain
// endflag == 2: CNT end at end of chain
if (endflag == 1)
geometry(r1, r2, p1, p2, qe, p, m, param, basis);
else
geometry(r1, r2, p2, p1, qe, p, m, param, basis);
if (param[0] > cutoff) continue;
if (!(param[2] < 0 && param[3] > 0)) {
double hsq = param[0] * param[0];
double calpha = cos(param[1]);
double etamin = calpha * param[2];
double dsq1;
if (etamin < param[6]) {
dsq1 = hsq + pow(sin(param[1]) * param[6], 2);
dsq1 += pow(param[2] - calpha * param[6], 2);
} else
dsq1 = hsq + pow(sin(param[1]) * param[2], 2);
etamin = calpha * param[3];
double dsq2;
if (etamin < param[6]) {
dsq2 = hsq + pow(sin(param[1]) * param[6], 2);
dsq2 += pow(param[3] - calpha * param[6], 2);
} else
dsq2 = hsq + pow(sin(param[1]) * param[3], 2);
if (dsq1 > cutoffsq && dsq2 > cutoffsq) continue;
}
fsemi(param, evdwl, fend, flocal);
}
if (evdwl == 0.0) continue;
evdwl *= 0.5;
// transform to global coordinate system
matvec(basis[0], basis[1], basis[2], flocal[0], fglobal[0]);
matvec(basis[0], basis[1], basis[2], flocal[1], fglobal[1]);
// mean chain velocity and relative velocity
add3(vp1, vp2, vp);
scale3(0.5 * sumw_inv, vp);
sub3(vp, vr, vrel);
vtot = dot3(vrel, basis[2]);
// friction forces
if (vtot == 0.0)
fvisc_tot = 0.0;
else {
fvisc_tot = fvisc_max / (1.0 + exp(-kvisc * (fabs(vtot) - vvisc))) - fvisc_shift;
fvisc_tot *= 0.25 * (1 - s(sumw)) * (param[3] - param[2]);
if (vtot < 0) fvisc_tot = -fvisc_tot;
}
scale3(fvisc_tot, basis[2], fvisc);
// add friction forces to current segment
add3(fvisc, f[i1], f[i1]);
add3(fvisc, f[i2], f[i2]);
// forces acting on approximate chain
add3(fglobal[0], fglobal[1], ftotal);
if (endflag) scaleadd3(fend, m, ftotal, ftotal);
scale3(-0.5, ftotal);
sub3(r1, p, delr1);
sub3(r2, p, delr2);
cross3(delr1, fglobal[0], t1);
cross3(delr2, fglobal[1], t2);
add3(t1, t2, torque);
// additional torque contribution from chain end
if (endflag) {
sub3(qe, p, delqe);
cross3(delqe, m, t3);
scale3(fend, t3);
add3(t3, torque, torque);
}
cross3(torque, m, ftorque);
lp = param[5] - param[4];
scale3(1.0 / lp, ftorque);
if (endflag == 2) {
add3(ftotal, ftorque, fglobal[3]);
sub3(ftotal, ftorque, fglobal[2]);
} else {
add3(ftotal, ftorque, fglobal[2]);
sub3(ftotal, ftorque, fglobal[3]);
}
scale3(0.5, fglobal[0]);
scale3(0.5, fglobal[1]);
scale3(0.5, fglobal[2]);
scale3(0.5, fglobal[3]);
// weight gradient terms acting on current segment
outer3(p1, dr1_sumw, temp);
minus3(q1_dr1_w, temp, dr1_p1);
outer3(p2, dr1_sumw, temp);
minus3(q2_dr1_w, temp, dr1_p2);
outer3(p1, dr2_sumw, temp);
minus3(q1_dr2_w, temp, dr2_p1);
outer3(p2, dr2_sumw, temp);
minus3(q2_dr2_w, temp, dr2_p2);
transpose_matvec(dr1_p1, fglobal[2], fgrad_r1_p1);
transpose_matvec(dr1_p2, fglobal[3], fgrad_r1_p2);
transpose_matvec(dr2_p1, fglobal[2], fgrad_r2_p1);
transpose_matvec(dr2_p2, fglobal[3], fgrad_r2_p2);
// add forces to nodes in current segment
add3(fglobal[0], f[i1], f[i1]);
add3(fglobal[1], f[i2], f[i2]);
scaleadd3(sumw_inv, fgrad_r1_p1, f[i1], f[i1]);
scaleadd3(sumw_inv, fgrad_r1_p2, f[i1], f[i1]);
scaleadd3(sumw_inv, fgrad_r2_p1, f[i2], f[i2]);
scaleadd3(sumw_inv, fgrad_r2_p2, f[i2], f[i2]);
// add forces in approximate chain
for (k = 0; k < clen - 1; k++) {
if (w[k] == 0.0) continue;
j1 = chain[j][k];
j2 = chain[j][k + 1];
j1 &= NEIGHMASK;
j2 &= NEIGHMASK;
scale = w[k] * sumw_inv;
scaleadd3(scale, fglobal[2], f[j1], f[j1]);
scaleadd3(scale, fglobal[3], f[j2], f[j2]);
// friction forces
scaleadd3(-scale, fvisc, f[j1], f[j1]);
scaleadd3(-scale, fvisc, f[j2], f[j2]);
}
// weight gradient terms acting on approximate chain
// iterate over nodes instead of segments
for (k = 0; k < clen; k++) {
if (wnode[k] == 0.0) continue;
j1 = chain[j][k];
j1 &= NEIGHMASK;
outer3(p1, dq_w[k], temp);
minus3(q1_dq_w[k], temp, dq_p1);
outer3(p2, dq_w[k], temp);
minus3(q2_dq_w[k], temp, dq_p2);
transpose_matvec(dq_p1, fglobal[2], fgrad_q_p1);
transpose_matvec(dq_p2, fglobal[3], fgrad_q_p2);
scaleadd3(sumw_inv, fgrad_q_p1, f[j1], f[j1]);
scaleadd3(sumw_inv, fgrad_q_p2, f[j1], f[j1]);
}
// force on node at CNT end
if (endflag) scaleadd3(0.5 * fend, m, f[endindex], f[endindex]);
// compute energy
if (eflag_either) {
if (eflag_global) eng_vdwl += evdwl;
if (eflag_atom) {
eatom[i1] += 0.25 * evdwl;
eatom[i2] += 0.25 * evdwl;
for (k = 0; k < clen - 1; k++) {
if (w[k] == 0.0) continue;
j1 = chain[j][k];
j2 = chain[j][k + 1];
j1 &= NEIGHMASK;
j2 &= NEIGHMASK;
evdwl_chain = 0.5 * w[k] * sumw_inv * evdwl;
eatom[j1] += evdwl_chain;
eatom[j2] += evdwl_chain;
}
}
}
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairMesoCNTViscous::coeff(int narg, char **arg)
{
if (narg < 6) utils::missing_cmd_args(FLERR, "pair_coeff", error);
read_file(arg[2]);
fvisc_max = utils::numeric(FLERR, arg[3], false, lmp);
kvisc = utils::numeric(FLERR, arg[4], false, lmp);
vvisc = utils::numeric(FLERR, arg[5], false, lmp);
fvisc_shift = fvisc_max / (1.0 + exp(kvisc * vvisc));
nend_types = narg - 6;
if (!allocated) allocate();
// end atom types
for (int i = 6; i < narg; i++) end_types[i - 6] = utils::inumeric(FLERR, arg[i], false, lmp);
// units, eV to energy unit conversion
ang = force->angstrom;
ang_inv = 1.0 / ang;
if (strcmp(update->unit_style, "real") == 0)
eunit = 23.06054966;
else if (strcmp(update->unit_style, "metal") == 0)
eunit = 1.0;
else if (strcmp(update->unit_style, "si") == 0)
eunit = 1.6021765e-19;
else if (strcmp(update->unit_style, "cgs") == 0)
eunit = 1.6021765e-12;
else if (strcmp(update->unit_style, "electron") == 0)
eunit = 3.674932248e-2;
else if (strcmp(update->unit_style, "micro") == 0)
eunit = 1.6021765e-4;
else if (strcmp(update->unit_style, "nano") == 0)
eunit = 1.6021765e2;
else
error->all(FLERR, "Pair style mesocnt/viscous does not support {} units", update->unit_style);
funit = eunit * ang_inv;
// potential variables
sig = sig_ang * ang;
r = r_ang * ang;
rsq = r * r;
d = 2.0 * r;
d_ang = 2.0 * r_ang;
rc = 3.0 * sig;
cutoff = rc + d;
cutoffsq = cutoff * cutoff;
cutoff_ang = cutoff * ang_inv;
cutoffsq_ang = cutoff_ang * cutoff_ang;
comega = 0.275 * (1.0 - 1.0 / (1.0 + 0.59 * r_ang));
ctheta = 0.35 + 0.0226 * (r_ang - 6.785);
// compute spline coefficients
spline_coeff(uinf_data, uinf_coeff, delh_uinf, uinf_points);
spline_coeff(gamma_data, gamma_coeff, delh_gamma, gamma_points);
spline_coeff(phi_data, phi_coeff, delh_phi, delpsi_phi, phi_points);
spline_coeff(usemi_data, usemi_coeff, delh_usemi, delxi_usemi, usemi_points);
memory->destroy(uinf_data);
memory->destroy(gamma_data);
memory->destroy(phi_data);
memory->destroy(usemi_data);
// compute Gauss-Legendre quadrature nodes and weights
gl_init_nodes(QUAD_FINF, gl_nodes_finf);
gl_init_nodes(QUAD_FSEMI, gl_nodes_fsemi);
gl_init_weights(QUAD_FINF, gl_nodes_finf, gl_weights_finf);
gl_init_weights(QUAD_FSEMI, gl_nodes_fsemi, gl_weights_fsemi);
int ntypes = atom->ntypes;
for (int i = 1; i <= ntypes; i++)
for (int j = i; j <= ntypes; j++) setflag[i][j] = 1;
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairMesoCNTViscous::init_style()
{
if (atom->tag_enable == 0) error->all(FLERR, "Pair style mesocnt/viscous requires atom IDs");
if (force->newton_pair == 0)
error->all(FLERR, "Pair style mesocnt/viscous requires newton pair on");
if (force->special_lj[1] == 0.0 || force->special_lj[2] == 0.0 || force->special_lj[3] == 0.0)
error->all(FLERR, "Pair mesocnt/viscous requires all special_bond lj values to be non-zero");
if (comm->ghost_velocity == 0)
error->all(FLERR, "Pair mesocnt/viscous requires ghost atoms store velocity");
// need a full neighbor list
neighbor->add_request(this, NeighConst::REQ_FULL);
}
/* ----------------------------------------------------------------------
weight for averaged friction from CNT chain
------------------------------------------------------------------------- */
inline double PairMesoCNTViscous::weight(const double *r1, const double *r2, const double *p1,
const double *p2)
{
double dr, dp, rhoc, rhomin, rho;
double r[3], p[3];
add3(r1, r2, r);
add3(p1, p2, p);
dr = sqrt(0.25 * distsq3(r1, r2) + rsq);
dp = sqrt(0.25 * distsq3(p1, p2) + rsq);
rhoc = dr + dp + rc;
rhomin = RHOMIN * ang;
rho = 0.5 * sqrt(distsq3(r, p));
return s((rho - rhomin) / (rhoc - rhomin));
}
/* ----------------------------------------------------------------------
weight for substitute CNT chain
computes gradients with respect to positions
------------------------------------------------------------------------- */
inline void PairMesoCNTViscous::weight(const double *r1, const double *r2, const double *p1,
const double *p2, double &w, double *dr1_w, double *dr2_w,
double *dp1_w, double *dp2_w)
{
double dr, dp, rhoc, rhomin, rho, frac, arg, factor;
double r[3], p[3];
double dr_rho[3], dr_rhoc[3], dp_rhoc[3];
add3(r1, r2, r);
add3(p1, p2, p);
scale3(0.5, r);
scale3(0.5, p);
dr = sqrt(0.25 * distsq3(r1, r2) + rsq);
dp = sqrt(0.25 * distsq3(p1, p2) + rsq);
rhoc = dr + dp + rc;
rhomin = RHOMIN * ang;
rho = sqrt(distsq3(r, p));
frac = 1.0 / (rhoc - rhomin);
arg = frac * (rho - rhomin);
w = s(arg);
if (w == 0.0 || w == 1.0) {
zero3(dr1_w);
zero3(dr2_w);
zero3(dp1_w);
zero3(dp2_w);
} else {
factor = ds(arg) * frac;
sub3(r, p, dr_rho);
sub3(r1, r2, dr_rhoc);
sub3(p1, p2, dp_rhoc);
scale3(0.5 / rho, dr_rho);
scale3(0.25 / dr, dr_rhoc);
scale3(0.25 / dp, dp_rhoc);
scaleadd3(-arg, dr_rhoc, dr_rho, dr1_w);
scaleadd3(arg, dr_rhoc, dr_rho, dr2_w);
negate3(dr_rho);
scaleadd3(-arg, dp_rhoc, dr_rho, dp1_w);
scaleadd3(arg, dp_rhoc, dr_rho, dp2_w);
scale3(factor, dr1_w);
scale3(factor, dr2_w);
scale3(factor, dp1_w);
scale3(factor, dp2_w);
}
}

50
src/MESONT/pair_mesocnt_viscous.h Normal file
View File

@ -0,0 +1,50 @@
/* -*- 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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Philipp Kloza (University of Cambridge)
pak37@cam.ac.uk
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(mesocnt/viscous, PairMesoCNTViscous);
// clang-format on
#else
#ifndef LMP_PAIR_MESOCNT_VISCOUS_H
#define LMP_PAIR_MESOCNT_VISCOUS_H
#include "pair_mesocnt.h"
namespace LAMMPS_NS {
class PairMesoCNTViscous : public PairMesoCNT {
public:
using PairMesoCNT::PairMesoCNT;
void compute(int, int) override;
void coeff(int, char **) override;
void init_style() override;
protected:
double fvisc_max, kvisc, vvisc, fvisc_shift;
inline double weight(const double *, const double *, const double *, const double *);
inline void weight(const double *, const double *, const double *, const double *, double &,
double *, double *, double *, double *);
};
} // namespace LAMMPS_NS
#endif
#endif

2
src/MESONT/potentials.txt
View File

@ -1,4 +1,4 @@
# list of potential files to be fetched when this package is installed
# potential file md5sum
C_10_10.mesocnt 028de73ec828b7830d762702eda571c1
C_10_10.mesocnt 68b5ca26283968fd9889aa0a37f7b7fb
TABTP_10_10.mesont 744a739da49ad5e78492c1fc9fd9f8c1