Merge branch 'develop' into gcmc-eatom
This commit is contained in:
Aidan Thompson
1
README
1
README
@ -34,6 +34,7 @@ lib additional provided or external libraries
|
||||
potentials interatomic potential files
|
||||
python Python module for LAMMPS
|
||||
src source files
|
||||
third_party Copies of thirdparty software bundled with LAMMPS
|
||||
tools pre- and post-processing tools
|
||||
unittest test programs for use with CTest
|
||||
.github Git and GitHub related files and tools
|
||||
|
||||
@ -44,6 +44,7 @@ set(LAMMPS_DOC_DIR ${LAMMPS_DIR}/doc)
|
||||
set(LAMMPS_TOOLS_DIR ${LAMMPS_DIR}/tools)
|
||||
set(LAMMPS_PYTHON_DIR ${LAMMPS_DIR}/python)
|
||||
set(LAMMPS_POTENTIALS_DIR ${LAMMPS_DIR}/potentials)
|
||||
set(LAMMPS_THIRDPARTY_DIR ${LAMMPS_DIR}/third_party)
|
||||
|
||||
set(LAMMPS_DOWNLOADS_URL "https://download.lammps.org" CACHE STRING "Base URL for LAMMPS downloads")
|
||||
set(LAMMPS_POTENTIALS_URL "${LAMMPS_DOWNLOADS_URL}/potentials")
|
||||
@ -131,7 +132,7 @@ endif()
|
||||
|
||||
# silence nvcc warnings
|
||||
if((PKG_KOKKOS) AND (Kokkos_ENABLE_CUDA) AND NOT (CMAKE_CXX_COMPILER_ID STREQUAL "Clang"))
|
||||
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}" "-Xcudafe --diag_suppress=unrecognized_pragma,--diag_suppress=128")
|
||||
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Xcudafe --diag_suppress=unrecognized_pragma,--diag_suppress=128")
|
||||
endif()
|
||||
|
||||
# we *require* C++11 without extensions but prefer C++17.
|
||||
@ -370,6 +371,7 @@ endforeach()
|
||||
# packages with special compiler needs or external libs
|
||||
######################################################
|
||||
target_include_directories(lammps PUBLIC $<BUILD_INTERFACE:${LAMMPS_SOURCE_DIR}>)
|
||||
target_include_directories(lammps PUBLIC $<BUILD_INTERFACE:${LAMMPS_THIRDPARTY_DIR}>)
|
||||
|
||||
if(PKG_ADIOS)
|
||||
# The search for ADIOS2 must come before MPI because
|
||||
|
||||
@ -30,7 +30,7 @@ function(check_omp_h_include)
|
||||
if(OpenMP_CXX_FOUND)
|
||||
set(CMAKE_REQUIRED_FLAGS ${OpenMP_CXX_FLAGS})
|
||||
set(CMAKE_REQUIRED_INCLUDES ${OpenMP_CXX_INCLUDE_DIRS})
|
||||
set(CMAKE_REQUIRED_LINK_OPTIONS ${OpenMP_CXX_FLAGS})
|
||||
separate_arguments(CMAKE_REQUIRED_LINK_OPTIONS NATIVE_COMMAND ${OpenMP_CXX_FLAGS}) # needs to be a list
|
||||
set(CMAKE_REQUIRED_LIBRARIES ${OpenMP_CXX_LIBRARIES})
|
||||
# there are all kinds of problems with finding omp.h
|
||||
# for Clang and derived compilers so we pretend it is there.
|
||||
|
||||
@ -23,7 +23,7 @@ set(MPI_CXX_COMPILER "mpicxx" CACHE STRING "" FORCE)
|
||||
set(CMAKE_CXX_STANDARD 17 CACHE STRING "" FORCE)
|
||||
|
||||
# set(_intel_sycl_flags " -w -fsycl -flink-huge-device-code -fsycl-targets=spir64_gen "
|
||||
set(_intel_sycl_flags " -w -fsycl -fsycl-device-code-split=per_kernel -fsycl-targets=spir64_gen "
|
||||
set(_intel_sycl_flags " -w -fsycl -fsycl-device-code-split=per_kernel -fsycl-targets=spir64_gen ")
|
||||
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${_intel_sycl_flags}" CACHE STRING "" FORCE)
|
||||
|
||||
#set(CMAKE_EXE_LINKER_FLAGS "-fsycl -flink-huge-device-code -fsycl-targets=spir64_gen " CACHE STRING "" FORCE)
|
||||
|
||||
@ -179,6 +179,7 @@ OPT.
|
||||
* :doc:`lj/long/dipole/long <pair_dipole>`
|
||||
* :doc:`lj/long/tip4p/long (o) <pair_lj_long>`
|
||||
* :doc:`lj/mdf <pair_mdf>`
|
||||
* :doc:`lj/pirani (o) <pair_lj_pirani>`
|
||||
* :doc:`lj/relres (o) <pair_lj_relres>`
|
||||
* :doc:`lj/spica (gko) <pair_spica>`
|
||||
* :doc:`lj/spica/coul/long (gko) <pair_spica>`
|
||||
|
||||
@ -12,19 +12,10 @@ several advantages:
|
||||
LAMMPS. For that, you should first create your own :doc:`fork on
|
||||
GitHub <Howto_github>`, though.
|
||||
|
||||
You must have `git <git_>`_ installed on your system to use the
|
||||
commands explained below to communicate with the git servers on
|
||||
GitHub. For people still using subversion (svn), GitHub also
|
||||
provides `limited support for subversion clients <svn_>`_.
|
||||
|
||||
.. note::
|
||||
|
||||
As of October 2016, the official home of public LAMMPS development is
|
||||
on GitHub. The previously advertised LAMMPS git repositories on
|
||||
git.lammps.org and bitbucket.org are now offline or deprecated.
|
||||
You must have `git <git_>`_ installed on your system to use the commands
|
||||
explained below to communicate with the git servers on GitHub.
|
||||
|
||||
.. _git: https://git-scm.com
|
||||
.. _svn: https://help.github.com/en/github/importing-your-projects-to-github/working-with-subversion-on-github
|
||||
|
||||
You can follow the LAMMPS development on 4 different git branches:
|
||||
|
||||
|
||||
@ -1250,10 +1250,10 @@ tabulate tool
|
||||
|
||||
.. versionadded:: 22Dec2022
|
||||
|
||||
The ``tabulate`` folder contains Python scripts scripts to generate tabulated
|
||||
potential files for LAMMPS. The bulk of the code is in the ``tabulate`` module
|
||||
in the ``tabulate.py`` file. Some example files demonstrating its use are
|
||||
included. See the README file for more information.
|
||||
The ``tabulate`` folder contains Python scripts scripts to generate and
|
||||
visualize tabulated potential files for LAMMPS. The bulk of the code is in the
|
||||
``tabulate`` module in the ``tabulate.py`` file. Some example files
|
||||
demonstrating its use are included. See the README file for more information.
|
||||
|
||||
----------
|
||||
|
||||
@ -1276,7 +1276,7 @@ Those scripts were written by Steve Plimpton sjplimp at gmail.com
|
||||
valgrind tool
|
||||
-------------
|
||||
|
||||
The ``valgrind`` folder contains additional suppressions fur LAMMPS when
|
||||
The ``valgrind`` folder contains additional suppressions for LAMMPS when
|
||||
using `valgrind's <https://valgrind.org/>`_ ` `memcheck tool
|
||||
<https://valgrind.org/info/tools.html#memcheck>`_ to search for memory
|
||||
access violation and memory leaks. These suppressions are automatically
|
||||
|
||||
@ -215,6 +215,9 @@ for an overview of LAMMPS output options.
|
||||
The vector or array will be floating point values that correspond to
|
||||
the specified attribute.
|
||||
|
||||
Any settings with the *store/local* option are not saved to a restart
|
||||
file and must be redefined.
|
||||
|
||||
The single() function of this bond style returns 0.0 for the energy
|
||||
of a bonded interaction, since energy is not conserved in these
|
||||
dissipative potentials. It also returns only the normal component of
|
||||
|
||||
@ -215,6 +215,9 @@ for an overview of LAMMPS output options.
|
||||
The vector or array will be floating point values that correspond to
|
||||
the specified attribute.
|
||||
|
||||
Any settings with the *store/local* option are not saved to a restart
|
||||
file and must be redefined.
|
||||
|
||||
The potential energy and the single() function of this bond style return
|
||||
:math:`k (r - r_0)^2 / 2` as a proxy of the energy of a bonded interaction,
|
||||
ignoring any volumetric/smoothing factors or dissipative forces. The single()
|
||||
|
||||
@ -88,7 +88,7 @@ The *phase* property indicates whether the particle is in a fluid state,
|
||||
a value of 0, or a solid state, a value of 1.
|
||||
|
||||
The *surface* property indicates the surface designation produced by
|
||||
the *interface/reconstruct* option of :doc:`fix rheo <fix_rheo>`. Bulk
|
||||
the *surface/detection* option of :doc:`fix rheo <fix_rheo>`. Bulk
|
||||
particles have a value of 0, surface particles have a value of 1, and
|
||||
splash particles have a value of 2. The *surface/r* property is the
|
||||
distance from the surface, up to the kernel cutoff length. Surface particles
|
||||
|
||||
@ -214,6 +214,8 @@ formulas for the meaning of these parameters:
|
||||
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
|
||||
| :doc:`lj/mdf <pair_mdf>` | epsilon,sigma | type pairs |
|
||||
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
|
||||
| :doc:`lj/pirani <pair_lj_pirani>` | alpha, beta, gamma, rm, epsilon | type pairs |
|
||||
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
|
||||
| :doc:`lj/sf/dipole/sf <pair_dipole>` | epsilon,sigma,scale | type pairs |
|
||||
+------------------------------------------------------------------------------+--------------------------------------------------+-------------+
|
||||
| :doc:`lubricate <pair_lubricate>` | mu | global |
|
||||
|
||||
163
doc/src/pair_lj_pirani.rst
Normal file
163
doc/src/pair_lj_pirani.rst
Normal file
@ -0,0 +1,163 @@
|
||||
.. index:: pair_style lj/pirani
|
||||
.. index:: pair_style lj/pirani/omp
|
||||
|
||||
pair_style lj/pirani command
|
||||
============================
|
||||
|
||||
Accelerator Variants: *lj/pirani/omp*
|
||||
|
||||
Syntax
|
||||
""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_style lj/pirani cutoff
|
||||
|
||||
* lj/pirani = name of the pair style
|
||||
* cutoff = global cutoff (distance units)
|
||||
|
||||
Examples
|
||||
""""""""
|
||||
|
||||
.. code-block:: LAMMPS
|
||||
|
||||
pair_style lj/pirani 10.0
|
||||
pair_coeff 1 1 4.0 7.0 6.0 3.5 0.0045
|
||||
|
||||
Description
|
||||
"""""""""""
|
||||
|
||||
.. versionadded:: TBD
|
||||
|
||||
Pair style *lj/pirani* computes pairwise interactions from an Improved
|
||||
Lennard-Jones (ILJ) potential according to :ref:`(Pirani) <Pirani>`.
|
||||
The ILJ force field is adequate to model both equilibrium and
|
||||
non-equilibrium properties of matter, in gaseous and condensed phases,
|
||||
and at gas-surface interfaces. In particular, its use improves the
|
||||
description of elementary process dynamics where the traditional
|
||||
Lennard-Jones (LJ) formulation is usually applied.
|
||||
|
||||
|
||||
.. math::
|
||||
|
||||
x = r/R_m \\
|
||||
n_x = \alpha*x^2 + \beta \\
|
||||
\gamma \equiv m \\
|
||||
|
||||
V(x) = \varepsilon \cdot \left( \frac{\gamma}{ n_x - \gamma} \left(\frac{1}{x} \right)^{n_x}
|
||||
- \frac{n_x}{n_x - \gamma} \left(\frac{1}{x} \right)^{\gamma} \right) \qquad r < r_c
|
||||
|
||||
:math:`r_c` is the cutoff.
|
||||
|
||||
|
||||
An additional parameter, :math:`\alpha`, has been introduced in order to
|
||||
be able to recover the traditional Lennard-Jones 12-6 with a specific
|
||||
choice of parameters. With :math:`R_m \equiv r_0 = \sigma \cdot 2^{1 /
|
||||
6}`, :math:`\alpha = 0`, :math:`\beta = 12` and :math:`\gamma = 6` it is
|
||||
straightforward to prove that LJ 12-6 is obtained. Also, it can be
|
||||
verified that using :math:`\alpha= 4`, :math:`\beta= 8` and
|
||||
:math:`\gamma = 6`, at the equilibrium distance, the first and second
|
||||
derivatives of ILJ match those of LJ 12-6. The parameter :math:`R_m`
|
||||
corresponds to the equilibrium distance and :math:`\epsilon` to the well
|
||||
depth.
|
||||
|
||||
|
||||
This potential provides some advantages with respect to the standard LJ
|
||||
potential, as explained in :ref:`(Pirani) <Pirani>`: it provides a more
|
||||
realistic description of the long range behavior and an attenuation of
|
||||
the hardness of the repulsive wall.
|
||||
|
||||
This force field can be used for neutral-neutral (:math:`\gamma = 6`),
|
||||
ion-neutral (:math:`\gamma = 4`) or ion-ion systems (:math:`\gamma =
|
||||
1`). Notice that this implementation does not include explicit
|
||||
electrostatic interactions. If these are desired, this pair style
|
||||
should be used along with a Coulomb pair style like
|
||||
:doc:`pair styles coul/cut or coul/long <pair_coul>` by using
|
||||
:doc:`pair style hybrid/overlay <pair_hybrid>` and a suitable
|
||||
kspace style :doc:`<kspace_style>`, if needed.
|
||||
|
||||
As discussed in :ref:`(Pirani) <Pirani>`, analysis of a variety of
|
||||
systems showed that :math:`\alpha= 4` generally works very well. In
|
||||
some special cases (e.g. those involving very small multiple charged
|
||||
ions) this factor may take a slightly different value. The parameter
|
||||
:math:`\beta` codifies the hardness (polarizability) of the interacting
|
||||
partners, and for neutral-neutral systems it usually ranges from 6
|
||||
to 11. Moreover, the modulation of :math:`\beta` can model additional
|
||||
interaction effects, such as charge transfer in the perturbative limit,
|
||||
and can mitigate the effect of some uncertainty in the data used to
|
||||
build up the potential function.
|
||||
|
||||
The following coefficients must be defined for each pair of atoms
|
||||
types via the :doc:`pair_coeff <pair_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:
|
||||
|
||||
* :math:`\alpha` (dimensionless)
|
||||
* :math:`\beta` (dimensionless)
|
||||
* :math:`\gamma` (dimensionless)
|
||||
* :math:`R_m` (distance units)
|
||||
* :math:`\epsilon` (energy units)
|
||||
* cutoff (distance units)
|
||||
|
||||
The last coefficient is optional. If not specified, the global cutoff is used.
|
||||
|
||||
----------
|
||||
|
||||
.. include:: accel_styles.rst
|
||||
|
||||
----------
|
||||
|
||||
Mixing, shift, table, tail correction, restart, rRESPA info
|
||||
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
|
||||
|
||||
This pair style does not support mixing. Thus, coefficients for all I,J
|
||||
pairs must be specified explicitly.
|
||||
|
||||
This pair style supports the :doc:`pair_modify <pair_modify>` shift
|
||||
option for the energy of the pair interaction.
|
||||
|
||||
The :doc:`pair_modify <pair_modify>` table options are not relevant for
|
||||
this pair style.
|
||||
|
||||
This pair style does not support the :doc:`pair_modify <pair_modify>`
|
||||
tail option for adding long-range tail corrections to energy and
|
||||
pressure.
|
||||
|
||||
This pair style writes its information to :doc:`binary restart files
|
||||
<restart>`, so pair_style and pair_coeff commands do not need to be
|
||||
specified in an input script that reads a restart file.
|
||||
|
||||
This pair style supports the use of the *inner*, *middle*, and
|
||||
*outer* keywords of the :doc:`run_style respa <run_style>` command,
|
||||
meaning the pairwise forces can be partitioned by distance at different
|
||||
levels of the rRESPA hierarchy. See the :doc:`run_style <run_style>`
|
||||
command for details.
|
||||
|
||||
|
||||
----------
|
||||
|
||||
Restrictions
|
||||
""""""""""""
|
||||
|
||||
This pair style is only enabled if LAMMPS was built with the EXTRA-PAIR
|
||||
package. See the :doc:`Build package <Build_package>` page for more
|
||||
info.
|
||||
|
||||
Related commands
|
||||
""""""""""""""""
|
||||
|
||||
* :doc:`pair_coeff <pair_coeff>`
|
||||
* :doc:`pair_style lj/cut <pair_lj>`
|
||||
|
||||
Default
|
||||
"""""""
|
||||
|
||||
none
|
||||
|
||||
--------------
|
||||
|
||||
.. _Pirani:
|
||||
|
||||
**(Pirani)** F. Pirani, S. Brizi, L. Roncaratti, P. Casavecchia, D. Cappelletti and F. Vecchiocattivi,
|
||||
Phys. Chem. Chem. Phys., 2008, 10, 5489-5503.
|
||||
@ -272,6 +272,7 @@ accelerated styles exist.
|
||||
* :doc:`lj/long/dipole/long <pair_dipole>` - long-range LJ and long-range point dipoles
|
||||
* :doc:`lj/long/tip4p/long <pair_lj_long>` - long-range LJ and long-range Coulomb for TIP4P water
|
||||
* :doc:`lj/mdf <pair_mdf>` - LJ potential with a taper function
|
||||
* :doc:`lj/pirani <pair_lj_pirani>` - Improved LJ potential
|
||||
* :doc:`lj/relres <pair_lj_relres>` - LJ using multiscale Relative Resolution (RelRes) methodology :ref:`(Chaimovich) <Chaimovich2>`.
|
||||
* :doc:`lj/spica <pair_spica>` - LJ for SPICA coarse-graining
|
||||
* :doc:`lj/spica/coul/long <pair_spica>` - LJ for SPICA coarse-graining with long-range Coulomb
|
||||
|
||||
@ -389,6 +389,7 @@ Bretonnet
|
||||
Briels
|
||||
Brien
|
||||
Brilliantov
|
||||
Brizi
|
||||
Broadwell
|
||||
Broglie
|
||||
brownian
|
||||
@ -431,10 +432,12 @@ Camiloni
|
||||
Campana
|
||||
Cangi
|
||||
Cao
|
||||
Cappelletti
|
||||
Capolungo
|
||||
Caro
|
||||
cartesian
|
||||
Cas
|
||||
Casavecchia
|
||||
CasP
|
||||
Caswell
|
||||
Cates
|
||||
@ -2922,6 +2925,7 @@ perp
|
||||
Perram
|
||||
persp
|
||||
Persp
|
||||
perturbative
|
||||
peru
|
||||
Peskin
|
||||
Pettifor
|
||||
@ -2960,6 +2964,8 @@ pimdb
|
||||
Piola
|
||||
pIp
|
||||
pipelining
|
||||
Pirani
|
||||
pirani
|
||||
Pisarev
|
||||
Pishevar
|
||||
Pitera
|
||||
@ -3346,6 +3352,7 @@ Rockett
|
||||
rocksalt
|
||||
Rodrigues
|
||||
Rohart
|
||||
Roncaratti
|
||||
Ronchetti
|
||||
Ronevich
|
||||
Rosati
|
||||
@ -4058,6 +4065,7 @@ vdW
|
||||
vdwl
|
||||
vec
|
||||
Vecchio
|
||||
Vecchiocattivi
|
||||
vectorial
|
||||
vectorization
|
||||
Vectorization
|
||||
|
||||
@ -1,18 +0,0 @@
|
||||
Python generated LAMMPS data file
|
||||
|
||||
2 atoms
|
||||
1 atom types
|
||||
|
||||
0 0.08 xlo xhi
|
||||
0 0.04 ylo yhi
|
||||
0 0.08 zlo zhi
|
||||
|
||||
Atoms # sphere
|
||||
|
||||
1 1 0.004 2500 0.04 0.02 0.04 0 0 0
|
||||
2 1 0.004 2500 0.04 0.02 0.04416 0 0 0
|
||||
|
||||
Velocities
|
||||
|
||||
1 0.0 0.0 1 0 0 0
|
||||
2 0.0 0.0 -1 0 0 0
|
||||
@ -1,24 +1,43 @@
|
||||
units si
|
||||
atom_style sphere
|
||||
comm_modify vel yes
|
||||
boundary p p p
|
||||
|
||||
boundary p p f
|
||||
region box block 0 80e-3 0 40e-3 0 80e-3 open 3 open 4
|
||||
region box block -0.01 0.01 -0.01 0.01 0 0.08
|
||||
create_box 2 box
|
||||
|
||||
read_data data.particles add append
|
||||
group mb type 1
|
||||
create_atoms 1 single 0.0 0.0 0.04
|
||||
create_atoms 1 single 0.0 0.0 0.04416
|
||||
set group all diameter 0.004 density 2500
|
||||
|
||||
pair_style granular
|
||||
pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping coeff_restitution
|
||||
# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping coeff_restitution
|
||||
comm_modify vel yes
|
||||
pair_style granular
|
||||
pair_coeff * * hertz/material 1e6 0.8 0.4 &
|
||||
tangential mindlin NULL 0.0 0.5 &
|
||||
damping coeff_restitution
|
||||
#pair_coeff * * hooke 1e6 0.5 &
|
||||
# tangential mindlin 1 1.0 0.0 &
|
||||
# damping coeff_restitution
|
||||
|
||||
timestep 1e-9
|
||||
fix 1 all nve/sphere
|
||||
compute s all stress/atom NULL pair
|
||||
|
||||
#dump 1 all custom 2000000 op.dump id x y z vx vy vz
|
||||
#dump_modify 1 pad 8
|
||||
thermo_style custom step ke
|
||||
run_style verlet
|
||||
run 10000000
|
||||
group a1 id 1
|
||||
group a2 id 2
|
||||
|
||||
velocity a1 set 0 0 1.0
|
||||
velocity a2 set 0 0 -1.0
|
||||
|
||||
compute z1 a1 reduce sum z
|
||||
compute z2 a2 reduce sum z
|
||||
compute v1 a1 reduce sum vz
|
||||
compute v2 a2 reduce sum vz
|
||||
compute f1 a1 reduce sum fz
|
||||
compute f2 a2 reduce sum fz
|
||||
variable dz equal c_z1 - c_z2
|
||||
|
||||
# dump 1 all custom 2000000 op.dump id x y z vx vy vz
|
||||
|
||||
thermo 10000
|
||||
thermo_style custom step ke v_dz c_v1 c_v2 c_f1 c_f2
|
||||
|
||||
run 1000000
|
||||
|
||||
@ -1,80 +0,0 @@
|
||||
LAMMPS (17 Apr 2024 - Development - patch_17Apr2024-93-g4e7bddaa0b)
|
||||
using 1 OpenMP thread(s) per MPI task
|
||||
units si
|
||||
atom_style sphere
|
||||
|
||||
boundary p p f
|
||||
region box block 0 80e-3 0 40e-3 0 80e-3 open 3 open 4
|
||||
create_box 2 box
|
||||
Created orthogonal box = (0 0 0) to (0.08 0.04 0.08)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
|
||||
read_data data.particles add append
|
||||
Reading data file ...
|
||||
orthogonal box = (0 0 0) to (0.08 0.04 0.08)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
reading atoms ...
|
||||
2 atoms
|
||||
reading velocities ...
|
||||
2 velocities
|
||||
read_data CPU = 0.002 seconds
|
||||
group mb type 1
|
||||
2 atoms in group mb
|
||||
|
||||
pair_style granular
|
||||
pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping coeff_restitution
|
||||
# pair_coeff * * hooke 1e6 0.5 tangential mindlin 1 1.0 0.0 damping coeff_restitution
|
||||
comm_modify vel yes
|
||||
|
||||
timestep 1e-9
|
||||
fix 1 all nve/sphere
|
||||
compute s all stress/atom NULL pair
|
||||
|
||||
#dump 1 all custom 2000000 op.dump id x y z vx vy vz
|
||||
#dump_modify 1 pad 8
|
||||
thermo_style custom step ke
|
||||
run_style verlet
|
||||
run 10000000
|
||||
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
|
||||
Neighbor list info ...
|
||||
update: every = 1 steps, delay = 0 steps, check = yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 0.005
|
||||
ghost atom cutoff = 0.005
|
||||
binsize = 0.0025, bins = 32 16 32
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair granular, perpetual
|
||||
attributes: half, newton on, size, history
|
||||
pair build: half/size/bin/atomonly/newton
|
||||
stencil: half/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 10.1 | 10.1 | 10.1 Mbytes
|
||||
Step KinEng
|
||||
0 8.3775804e-05
|
||||
10000000 5.3616513e-05
|
||||
Loop time of 5.99782 on 1 procs for 10000000 steps with 2 atoms
|
||||
|
||||
77.9% CPU use with 1 MPI tasks x 1 OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.60235 | 0.60235 | 0.60235 | 0.0 | 10.04
|
||||
Neigh | 0.00021965 | 0.00021965 | 0.00021965 | 0.0 | 0.00
|
||||
Comm | 1.7939 | 1.7939 | 1.7939 | 0.0 | 29.91
|
||||
Output | 2.5955e-05 | 2.5955e-05 | 2.5955e-05 | 0.0 | 0.00
|
||||
Modify | 1.7622 | 1.7622 | 1.7622 | 0.0 | 29.38
|
||||
Other | | 1.839 | | | 30.66
|
||||
|
||||
Nlocal: 2 ave 2 max 2 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
|
||||
|
||||
Total # of neighbors = 0
|
||||
Ave neighs/atom = 0
|
||||
Neighbor list builds = 14
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:06
|
||||
195
examples/granular/log.4Feb25.restitution.g++.1
Normal file
195
examples/granular/log.4Feb25.restitution.g++.1
Normal file
@ -0,0 +1,195 @@
|
||||
LAMMPS (4 Feb 2025)
|
||||
units si
|
||||
atom_style sphere
|
||||
comm_modify vel yes
|
||||
boundary p p p
|
||||
|
||||
region box block -0.01 0.01 -0.01 0.01 0 0.08
|
||||
create_box 2 box
|
||||
Created orthogonal box = (-0.01 -0.01 0) to (0.01 0.01 0.08)
|
||||
1 by 1 by 1 MPI processor grid
|
||||
|
||||
create_atoms 1 single 0.0 0.0 0.04
|
||||
Created 1 atoms
|
||||
using lattice units in orthogonal box = (-0.01 -0.01 0) to (0.01 0.01 0.08)
|
||||
create_atoms CPU = 0.000 seconds
|
||||
create_atoms 1 single 0.0 0.0 0.04416
|
||||
Created 1 atoms
|
||||
using lattice units in orthogonal box = (-0.01 -0.01 0) to (0.01 0.01 0.08)
|
||||
create_atoms CPU = 0.000 seconds
|
||||
set group all diameter 0.004 density 2500
|
||||
Setting atom values ...
|
||||
2 settings made for diameter
|
||||
2 settings made for density
|
||||
|
||||
pair_style granular
|
||||
pair_coeff * * hertz/material 1e6 0.8 0.4 tangential mindlin NULL 0.0 0.5 damping coeff_restitution
|
||||
#pair_coeff * * hooke 1e6 0.5 # tangential mindlin 1 1.0 0.0 # damping coeff_restitution
|
||||
|
||||
timestep 1e-9
|
||||
fix 1 all nve/sphere
|
||||
|
||||
group a1 id 1
|
||||
1 atoms in group a1
|
||||
group a2 id 2
|
||||
1 atoms in group a2
|
||||
|
||||
velocity a1 set 0 0 1.0
|
||||
velocity a2 set 0 0 -1.0
|
||||
|
||||
compute z1 a1 reduce sum z
|
||||
compute z2 a2 reduce sum z
|
||||
compute v1 a1 reduce sum vz
|
||||
compute v2 a2 reduce sum vz
|
||||
compute f1 a1 reduce sum fz
|
||||
compute f2 a2 reduce sum fz
|
||||
variable dz equal c_z1 - c_z2
|
||||
|
||||
# dump 1 all custom 2000000 op.dump id x y z vx vy vz
|
||||
|
||||
thermo 10000
|
||||
thermo_style custom step ke v_dz c_v1 c_v2 c_f1 c_f2
|
||||
|
||||
run 1000000
|
||||
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
|
||||
Neighbor list info ...
|
||||
update: every = 1 steps, delay = 0 steps, check = yes
|
||||
max neighbors/atom: 2000, page size: 100000
|
||||
master list distance cutoff = 0.005
|
||||
ghost atom cutoff = 0.005
|
||||
binsize = 0.0025, bins = 8 8 32
|
||||
1 neighbor lists, perpetual/occasional/extra = 1 0 0
|
||||
(1) pair granular, perpetual
|
||||
attributes: half, newton on, size, history
|
||||
pair build: half/size/bin/atomonly/newton
|
||||
stencil: half/bin/3d
|
||||
bin: standard
|
||||
Per MPI rank memory allocation (min/avg/max) = 9.998 | 9.998 | 9.998 Mbytes
|
||||
Step KinEng v_dz c_v1 c_v2 c_f1 c_f2
|
||||
0 8.3775804e-05 -0.00416 1 -1 0 0
|
||||
10000 8.3775804e-05 -0.00414 1 -1 0 0
|
||||
20000 8.3775804e-05 -0.00412 1 -1 0 0
|
||||
30000 8.3775804e-05 -0.0041 1 -1 0 0
|
||||
40000 8.3775804e-05 -0.00408 1 -1 0 0
|
||||
50000 8.3775804e-05 -0.00406 1 -1 0 0
|
||||
60000 8.3775804e-05 -0.00404 1 -1 0 0
|
||||
70000 8.3775804e-05 -0.00402 1 -1 0 0
|
||||
80000 8.3775804e-05 -0.004 1 -1 0 0
|
||||
90000 8.3411065e-05 -0.003980019 0.99782075 -0.99782075 -0.023914747 0.023914747
|
||||
100000 8.2852688e-05 -0.0039600945 0.9944753 -0.9944753 -0.032005131 0.032005131
|
||||
110000 8.2139641e-05 -0.0039402463 0.99018672 -0.99018672 -0.039875404 0.039875404
|
||||
120000 8.1272296e-05 -0.0039204934 0.98494496 -0.98494496 -0.048007824 0.048007824
|
||||
130000 8.0246788e-05 -0.0039008551 0.97871113 -0.97871113 -0.056503872 0.056503872
|
||||
140000 7.9058986e-05 -0.0038813518 0.97144075 -0.97144075 -0.065373554 0.065373554
|
||||
150000 7.7705654e-05 -0.0038620047 0.9630903 -0.9630903 -0.074595879 0.074595879
|
||||
160000 7.6184906e-05 -0.0038428357 0.95361959 -0.95361959 -0.084137355 0.084137355
|
||||
170000 7.4496418e-05 -0.0038238676 0.94299284 -0.94299284 -0.093958893 0.093958893
|
||||
180000 7.2641536e-05 -0.0038051239 0.93117907 -0.93117907 -0.10401872 0.10401872
|
||||
190000 7.0623328e-05 -0.0037866285 0.91815243 -0.91815243 -0.11427372 0.11427372
|
||||
200000 6.8446602e-05 -0.003768406 0.90389221 -0.90389221 -0.12468011 0.12468011
|
||||
210000 6.6117901e-05 -0.0037504812 0.88838298 -0.88838298 -0.13519381 0.13519381
|
||||
220000 6.3645478e-05 -0.0037328791 0.87161455 -0.87161455 -0.14577066 0.14577066
|
||||
230000 6.1039243e-05 -0.003715625 0.85358204 -0.85358204 -0.1563666 0.1563666
|
||||
240000 5.8310702e-05 -0.0036987442 0.83428576 -0.83428576 -0.16693776 0.16693776
|
||||
250000 5.5472871e-05 -0.003682262 0.8137313 -0.8137313 -0.17744066 0.17744066
|
||||
260000 5.2540172e-05 -0.0036662033 0.79192936 -0.79192936 -0.18783225 0.18783225
|
||||
270000 4.9528314e-05 -0.003650593 0.76889577 -0.76889577 -0.19807012 0.19807012
|
||||
280000 4.6454158e-05 -0.0036354555 0.74465137 -0.74465137 -0.20811258 0.20811258
|
||||
290000 4.3335566e-05 -0.0036208149 0.71922195 -0.71922195 -0.21791884 0.21791884
|
||||
300000 4.0191232e-05 -0.0036066944 0.69263807 -0.69263807 -0.22744912 0.22744912
|
||||
310000 3.7040511e-05 -0.0035931168 0.66493499 -0.66493499 -0.23666485 0.23666485
|
||||
320000 3.390323e-05 -0.0035801042 0.63615249 -0.63615249 -0.24552878 0.24552878
|
||||
330000 3.0799488e-05 -0.0035676776 0.60633473 -0.60633473 -0.25400513 0.25400513
|
||||
340000 2.7749462e-05 -0.0035558573 0.57553002 -0.57553002 -0.26205979 0.26205979
|
||||
350000 2.4773197e-05 -0.0035446626 0.54379063 -0.54379063 -0.2696604 0.2696604
|
||||
360000 2.1890403e-05 -0.0035341116 0.51117262 -0.51117262 -0.27677654 0.27677654
|
||||
370000 1.9120254e-05 -0.0035242212 0.47773551 -0.47773551 -0.28337983 0.28337983
|
||||
380000 1.6481181e-05 -0.0035150072 0.44354212 -0.44354212 -0.28944409 0.28944409
|
||||
390000 1.3990689e-05 -0.0035064841 0.40865822 -0.40865822 -0.29494544 0.29494544
|
||||
400000 1.1665166e-05 -0.003498665 0.37315233 -0.37315233 -0.2998624 0.2998624
|
||||
410000 9.5197195e-06 -0.0034915617 0.33709536 -0.33709536 -0.304176 0.304176
|
||||
420000 7.5680136e-06 -0.0034851844 0.30056032 -0.30056032 -0.30786987 0.30786987
|
||||
430000 5.8221324e-06 -0.003479542 0.26362205 -0.26362205 -0.31093026 0.31093026
|
||||
440000 4.2924559e-06 -0.0034746417 0.22635684 -0.22635684 -0.31334618 0.31334618
|
||||
450000 2.9875585e-06 -0.0034704894 0.18884214 -0.18884214 -0.31510936 0.31510936
|
||||
460000 1.9141264e-06 -0.0034670891 0.15115621 -0.15115621 -0.31621432 0.31621432
|
||||
470000 1.0768988e-06 -0.0034644437 0.11337783 -0.11337783 -0.31665837 0.31665837
|
||||
480000 4.786302e-07 -0.0034625541 0.075585893 -0.075585893 -0.31644161 0.31644161
|
||||
490000 1.2007709e-07 -0.0034614198 0.037859142 -0.037859142 -0.31556689 0.31556689
|
||||
500000 6.3727744e-12 -0.0034610388 0.0002758068 -0.0002758068 -0.31403982 0.31403982
|
||||
510000 1.1522726e-07 -0.0034614073 -0.03708671 0.03708671 -0.31186864 0.31186864
|
||||
520000 4.6064472e-07 -0.0034625203 -0.074152149 0.074152149 -0.30906426 0.30906426
|
||||
530000 1.029334e-06 -0.0034643709 -0.1108457 0.1108457 -0.30564009 0.30564009
|
||||
540000 1.812635e-06 -0.0034669512 -0.14709431 0.14709431 -0.30161199 0.30161199
|
||||
550000 2.8002645e-06 -0.0034702513 -0.18282695 0.18282695 -0.29699817 0.29699817
|
||||
560000 3.9804448e-06 -0.0034742603 -0.21797491 0.21797491 -0.29181905 0.29181905
|
||||
570000 5.3400475e-06 -0.0034789659 -0.25247202 0.25247202 -0.28609717 0.28609717
|
||||
580000 6.8647484e-06 -0.0034843545 -0.28625495 0.28625495 -0.27985701 0.27985701
|
||||
590000 8.5391931e-06 -0.003490411 -0.31926339 0.31926339 -0.2731249 0.2731249
|
||||
600000 1.0347171e-05 -0.0034971194 -0.35144026 0.35144026 -0.2659288 0.2659288
|
||||
610000 1.2271792e-05 -0.0035044627 -0.38273192 0.38273192 -0.25829822 0.25829822
|
||||
620000 1.429567e-05 -0.0035124225 -0.41308835 0.41308835 -0.25026402 0.25026402
|
||||
630000 1.6401103e-05 -0.0035209797 -0.44246327 0.44246327 -0.24185823 0.24185823
|
||||
640000 1.8570255e-05 -0.0035301142 -0.47081428 0.47081428 -0.23311393 0.23311393
|
||||
650000 2.078533e-05 -0.0035398052 -0.498103 0.498103 -0.22406507 0.22406507
|
||||
660000 2.3028745e-05 -0.003550031 -0.52429514 0.52429514 -0.21474628 0.21474628
|
||||
670000 2.5283288e-05 -0.0035607695 -0.54936056 0.54936056 -0.20519274 0.20519274
|
||||
680000 2.7532278e-05 -0.0035719977 -0.57327337 0.57327337 -0.19544002 0.19544002
|
||||
690000 2.9759697e-05 -0.0035836926 -0.59601193 0.59601193 -0.18552392 0.18552392
|
||||
700000 3.1950329e-05 -0.0035958303 -0.61755887 0.61755887 -0.17548034 0.17548034
|
||||
710000 3.408987e-05 -0.0036083869 -0.63790112 0.63790112 -0.16534511 0.16534511
|
||||
720000 3.6165032e-05 -0.0036213382 -0.65702988 0.65702988 -0.15515392 0.15515392
|
||||
730000 3.8163631e-05 -0.00363466 -0.67494058 0.67494058 -0.14494218 0.14494218
|
||||
740000 4.0074659e-05 -0.0036483277 -0.69163285 0.69163285 -0.13474489 0.13474489
|
||||
750000 4.1888343e-05 -0.0036623172 -0.7071105 0.7071105 -0.1245966 0.1245966
|
||||
760000 4.3596185e-05 -0.0036766041 -0.7213814 0.7213814 -0.11453133 0.11453133
|
||||
770000 4.5190996e-05 -0.0036911645 -0.73445747 0.73445747 -0.1045825 0.1045825
|
||||
780000 4.6666906e-05 -0.0037059745 -0.74635458 0.74635458 -0.094782939 0.094782939
|
||||
790000 4.8019368e-05 -0.0037210109 -0.75709246 0.75709246 -0.085164857 0.085164857
|
||||
800000 4.9245153e-05 -0.0037362507 -0.76669467 0.76669467 -0.075759891 0.075759891
|
||||
810000 5.0342325e-05 -0.0037516713 -0.77518852 0.77518852 -0.066599178 0.066599178
|
||||
820000 5.1310215e-05 -0.003767251 -0.78260499 0.78260499 -0.057713474 0.057713474
|
||||
830000 5.2149388e-05 -0.0037829686 -0.78897875 0.78897875 -0.049133335 0.049133335
|
||||
840000 5.2861601e-05 -0.0037988035 -0.79434809 0.79434809 -0.040889384 0.040889384
|
||||
850000 5.3449761e-05 -0.0038147361 -0.79875498 0.79875498 -0.03301269 0.03301269
|
||||
860000 5.3917883e-05 -0.0038307476 -0.80224517 0.80224517 -0.025535302 0.025535302
|
||||
870000 5.4271056e-05 -0.0038468201 -0.80486832 0.80486832 -0.018491033 0.018491033
|
||||
880000 5.4515415e-05 -0.0038629369 -0.80667827 0.80667827 -0.011916591 0.011916591
|
||||
890000 5.4658137e-05 -0.0038790822 -0.80773352 0.80773352 -0.00585333 0.00585333
|
||||
900000 5.4707463e-05 -0.0038952416 -0.80809791 0.80809791 -0.00035001143 0.00035001143
|
||||
910000 5.4672786e-05 -0.003911402 -0.80784176 0.80784176 0.0045325009 -0.0045325009
|
||||
920000 5.4564826e-05 -0.0039275517 -0.80704376 0.80704376 0.0087126344 -0.0087126344
|
||||
930000 5.4396015e-05 -0.0039436807 -0.80579439 0.80579439 0.012070422 -0.012070422
|
||||
940000 5.4181322e-05 -0.0039597811 -0.80420264 0.80420264 0.014404618 -0.014404618
|
||||
950000 5.3940325e-05 -0.0039758475 -0.80241211 0.80241211 0.015295834 -0.015295834
|
||||
960000 5.3704817e-05 -0.0039918778 -0.8006585 0.8006585 0.013305187 -0.013305187
|
||||
970000 5.3616513e-05 -0.0040078808 -0.79999999 0.79999999 0 0
|
||||
980000 5.3616513e-05 -0.0040238808 -0.79999999 0.79999999 0 0
|
||||
990000 5.3616513e-05 -0.0040398808 -0.79999999 0.79999999 0 0
|
||||
1000000 5.3616513e-05 -0.0040558808 -0.79999999 0.79999999 0 0
|
||||
Loop time of 0.572025 on 1 procs for 1000000 steps with 2 atoms
|
||||
|
||||
99.3% CPU use with 1 MPI tasks x no OpenMP threads
|
||||
|
||||
MPI task timing breakdown:
|
||||
Section | min time | avg time | max time |%varavg| %total
|
||||
---------------------------------------------------------------
|
||||
Pair | 0.1248 | 0.1248 | 0.1248 | 0.0 | 21.82
|
||||
Neigh | 0 | 0 | 0 | 0.0 | 0.00
|
||||
Comm | 0.21166 | 0.21166 | 0.21166 | 0.0 | 37.00
|
||||
Output | 0.00057419 | 0.00057419 | 0.00057419 | 0.0 | 0.10
|
||||
Modify | 0.12695 | 0.12695 | 0.12695 | 0.0 | 22.19
|
||||
Other | | 0.108 | | | 18.89
|
||||
|
||||
Nlocal: 2 ave 2 max 2 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: 1 ave 1 max 1 min
|
||||
Histogram: 1 0 0 0 0 0 0 0 0 0
|
||||
|
||||
Total # of neighbors = 1
|
||||
Ave neighs/atom = 0.5
|
||||
Neighbor list builds = 0
|
||||
Dangerous builds = 0
|
||||
Total wall time: 0:00:01
|
||||
2
src/.gitignore
vendored
2
src/.gitignore
vendored
@ -1801,4 +1801,6 @@
|
||||
/pair_smtbq.h
|
||||
/pair_vashishta*.cpp
|
||||
/pair_vashishta*.h
|
||||
/pair_lj_pirani.cpp
|
||||
/pair_lj_pirani.h
|
||||
|
||||
|
||||
@ -397,6 +397,7 @@ double BondBPM::equilibrium_distance(int /*i*/)
|
||||
void BondBPM::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&overlay_flag, sizeof(int), 1, fp);
|
||||
fwrite(&break_flag, sizeof(int), 1, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -405,8 +406,12 @@ void BondBPM::write_restart(FILE *fp)
|
||||
|
||||
void BondBPM::read_restart(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) utils::sfread(FLERR, &overlay_flag, sizeof(int), 1, fp, nullptr, error);
|
||||
if (comm->me == 0) {
|
||||
utils::sfread(FLERR, &overlay_flag, sizeof(int), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &break_flag, sizeof(int), 1, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&overlay_flag, 1, MPI_INT, 0, world);
|
||||
MPI_Bcast(&break_flag, 1, MPI_INT, 0, world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -28,6 +28,7 @@
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "neighbor.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
@ -483,6 +484,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
|
||||
int newton_bond = force->newton_bond;
|
||||
|
||||
double **bondstore = fix_bond_history->bondstore;
|
||||
const bool allow_breaks = (update->setupflag == 0) && break_flag;
|
||||
|
||||
for (n = 0; n < nbondlist; n++) {
|
||||
|
||||
@ -527,7 +529,7 @@ void BondBPMRotational::compute(int eflag, int vflag)
|
||||
breaking = elastic_forces(i1, i2, type, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
|
||||
torque1on2, torque2on1);
|
||||
|
||||
if ((breaking >= 1.0) && break_flag) {
|
||||
if ((breaking >= 1.0) && allow_breaks) {
|
||||
bondlist[n][2] = 0;
|
||||
process_broken(i1, i2);
|
||||
continue;
|
||||
@ -763,6 +765,7 @@ void BondBPMRotational::read_restart(FILE *fp)
|
||||
void BondBPMRotational::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&smooth_flag, sizeof(int), 1, fp);
|
||||
fwrite(&normalize_flag, sizeof(int), 1, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -771,8 +774,12 @@ void BondBPMRotational::write_restart_settings(FILE *fp)
|
||||
|
||||
void BondBPMRotational::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) utils::sfread(FLERR, &smooth_flag, sizeof(int), 1, fp, nullptr, error);
|
||||
if (comm->me == 0) {
|
||||
utils::sfread(FLERR, &smooth_flag, sizeof(int), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &normalize_flag, sizeof(int), 1, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&smooth_flag, 1, MPI_INT, 0, world);
|
||||
MPI_Bcast(&normalize_flag, 1, MPI_INT, 0, world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -26,6 +26,7 @@
|
||||
#include "memory.h"
|
||||
#include "modify.h"
|
||||
#include "neighbor.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
@ -218,6 +219,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
|
||||
double invdim = 1.0 / dim;
|
||||
|
||||
double **bondstore = fix_bond_history->bondstore;
|
||||
const bool allow_breaks = (update->setupflag == 0) && break_flag;
|
||||
|
||||
for (n = 0; n < nbondlist; n++) {
|
||||
|
||||
@ -249,7 +251,7 @@ void BondBPMSpring::compute(int eflag, int vflag)
|
||||
r = sqrt(rsq);
|
||||
e = (r - r0) / r0;
|
||||
|
||||
if ((fabs(e) > ecrit[type]) && break_flag) {
|
||||
if ((fabs(e) > ecrit[type]) && allow_breaks) {
|
||||
bondlist[n][2] = 0;
|
||||
process_broken(i1, i2);
|
||||
|
||||
|
||||
@ -26,6 +26,7 @@
|
||||
#include "force.h"
|
||||
#include "memory.h"
|
||||
#include "neighbor.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
@ -185,6 +186,7 @@ void BondBPMSpringPlastic::compute(int eflag, int vflag)
|
||||
int newton_bond = force->newton_bond;
|
||||
|
||||
double **bondstore = fix_bond_history->bondstore;
|
||||
const bool allow_breaks = (update->setupflag == 0) && break_flag;
|
||||
|
||||
for (n = 0; n < nbondlist; n++) {
|
||||
|
||||
@ -217,7 +219,7 @@ void BondBPMSpringPlastic::compute(int eflag, int vflag)
|
||||
r = sqrt(rsq);
|
||||
e = (r - r0) / r0;
|
||||
|
||||
if ((fabs(e) > ecrit[type]) && break_flag) {
|
||||
if ((fabs(e) > ecrit[type]) && allow_breaks) {
|
||||
bondlist[n][2] = 0;
|
||||
process_broken(i1, i2);
|
||||
continue;
|
||||
|
||||
@ -26,6 +26,7 @@
|
||||
#include "neigh_list.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <set>
|
||||
#include <utility>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
@ -267,18 +268,40 @@ void FixUpdateSpecialBonds::post_run()
|
||||
|
||||
void FixUpdateSpecialBonds::add_broken_bond(int i, int j)
|
||||
{
|
||||
auto tag_pair = std::make_pair(atom->tag[i], atom->tag[j]);
|
||||
tagint *tag = atom->tag;
|
||||
int mintag = MIN(tag[i], tag[j]);
|
||||
int maxtag = MAX(tag[i], tag[j]);
|
||||
auto tag_pair = std::make_pair(mintag, maxtag);
|
||||
new_broken_pairs.push_back(tag_pair);
|
||||
broken_pairs.push_back(tag_pair);
|
||||
|
||||
// cancel out if created->destroyed before nlist rebuild
|
||||
// however, still cannot break + create in the same timestep
|
||||
auto const &it = created_pairs.find(tag_pair);
|
||||
if (it != created_pairs.end()) {
|
||||
created_pairs.erase(it);
|
||||
} else {
|
||||
broken_pairs.insert(tag_pair);
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void FixUpdateSpecialBonds::add_created_bond(int i, int j)
|
||||
{
|
||||
auto tag_pair = std::make_pair(atom->tag[i], atom->tag[j]);
|
||||
tagint *tag = atom->tag;
|
||||
int mintag = MIN(tag[i], tag[j]);
|
||||
int maxtag = MAX(tag[i], tag[j]);
|
||||
auto tag_pair = std::make_pair(mintag, maxtag);
|
||||
new_created_pairs.push_back(tag_pair);
|
||||
created_pairs.push_back(tag_pair);
|
||||
|
||||
// cancel out if destroyed->created before nlist rebuild
|
||||
// however, still cannot break + create in the same timestep
|
||||
auto const &it = broken_pairs.find(tag_pair);
|
||||
if (it != broken_pairs.end()) {
|
||||
broken_pairs.erase(it);
|
||||
} else {
|
||||
created_pairs.insert(tag_pair);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
||||
@ -22,6 +22,7 @@ FixStyle(UPDATE_SPECIAL_BONDS,FixUpdateSpecialBonds);
|
||||
|
||||
#include "fix.h"
|
||||
|
||||
#include <set>
|
||||
#include <utility>
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
@ -39,15 +40,15 @@ class FixUpdateSpecialBonds : public Fix {
|
||||
void write_restart(FILE *) override;
|
||||
|
||||
protected:
|
||||
// Create two arrays to store bonds broken this timestep (new)
|
||||
// and since the last neighbor list build
|
||||
// Create array to store bonds broken this timestep (new)
|
||||
// and a set for those broken since the last neighbor list build
|
||||
std::vector<std::pair<tagint, tagint>> new_broken_pairs;
|
||||
std::vector<std::pair<tagint, tagint>> broken_pairs;
|
||||
std::set<std::pair<tagint, tagint>> broken_pairs;
|
||||
|
||||
// Create two arrays to store newly created this timestep (new)
|
||||
// and since the last neighbor list build
|
||||
// Create arrays to store newly created this timestep (new)
|
||||
// and a set for those created since the last neighbor list build
|
||||
std::vector<std::pair<tagint, tagint>> new_created_pairs;
|
||||
std::vector<std::pair<tagint, tagint>> created_pairs;
|
||||
std::set<std::pair<tagint, tagint>> created_pairs;
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
888
src/EXTRA-PAIR/pair_lj_pirani.cpp
Normal file
888
src/EXTRA-PAIR/pair_lj_pirani.cpp
Normal file
@ -0,0 +1,888 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
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: Mateo Rodríguez (mateorsuarez@gmail.com) (IFF-CSIC)
|
||||
Work done at the Molecular Interactions Group (INTERMOL) of the
|
||||
Fundamental Physics Institute (http://intermol.iff.csic.es/).
|
||||
Optimization of the code: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "pair_lj_pirani.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "memory.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neighbor.h"
|
||||
#include "respa.h"
|
||||
#include "update.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using MathSpecial::square;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJPirani::PairLJPirani(LAMMPS *lmp) : Pair(lmp)
|
||||
{
|
||||
respa_enable = 1;
|
||||
born_matrix_enable = 0;
|
||||
writedata = 1;
|
||||
}
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJPirani::~PairLJPirani()
|
||||
{
|
||||
if (copymode) return;
|
||||
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(cutsq);
|
||||
memory->destroy(cut);
|
||||
memory->destroy(alpha);
|
||||
memory->destroy(beta);
|
||||
memory->destroy(gamma);
|
||||
memory->destroy(rm);
|
||||
memory->destroy(epsilon);
|
||||
memory->destroy(offset);
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
allocate all arrays
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->ntypes + 1;
|
||||
|
||||
memory->create(setflag, n, n, "pair:setflag");
|
||||
for (int i = 1; i < n; i++)
|
||||
for (int j = i; j < n; j++) setflag[i][j] = 0;
|
||||
|
||||
memory->create(cutsq, n, n, "pair:cutsq");
|
||||
|
||||
memory->create(cut, n, n, "pair:cut");
|
||||
memory->create(alpha, n, n, "pair:alpha");
|
||||
memory->create(beta, n, n, "pair:beta");
|
||||
memory->create(gamma, n, n, "pair:gamma");
|
||||
memory->create(rm, n, n, "pair:rm");
|
||||
memory->create(epsilon, n, n, "pair:epsilon");
|
||||
memory->create(offset, n, n, "pair:offset");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::compute(int eflag, int vflag)
|
||||
|
||||
{
|
||||
int i, j, ii, jj, inum, jnum, itype, jtype;
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
|
||||
double rsq, forcelj, factor_lj;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
double r, rx, n_x;
|
||||
double pow_rx_n_x, pow_rx_gamma;
|
||||
double filj1, filj2, filj3, filj4, filj5, filj6, forceilj;
|
||||
double ilj1, ilj2;
|
||||
double fxtmp, fytmp, fztmp;
|
||||
|
||||
evdwl = 0.0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
|
||||
r = sqrt(rsq);
|
||||
|
||||
rx = r / rm[itype][jtype];
|
||||
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
|
||||
pow_rx_n_x = pow(1.0 / rx, n_x);
|
||||
pow_rx_gamma = pow(1.0 / rx, gamma[itype][jtype]);
|
||||
|
||||
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow_rx_n_x /
|
||||
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
|
||||
|
||||
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow_rx_gamma /
|
||||
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
|
||||
|
||||
filj3 = -2.0 * alpha[itype][jtype] * rx * pow_rx_gamma /
|
||||
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) *
|
||||
log(1 / rx) * pow_rx_n_x / (n_x - gamma[itype][jtype]);
|
||||
|
||||
filj5 = -1.0 * gamma[itype][jtype] * n_x * pow_rx_n_x / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj6 = +1.0 * gamma[itype][jtype] * n_x * pow_rx_gamma / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
// F = -dV/dr
|
||||
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
|
||||
fpair = factor_lj * forceilj / r; // F_x = -x/r * dV/dr (chain rule)
|
||||
|
||||
fxtmp += delx * fpair;
|
||||
fytmp += dely * fpair;
|
||||
fztmp += delz * fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx * fpair;
|
||||
f[j][1] -= dely * fpair;
|
||||
f[j][2] -= delz * fpair;
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
ilj1 = epsilon[itype][jtype] * gamma[itype][jtype] * pow(1 / rx, n_x) /
|
||||
(n_x - gamma[itype][jtype]);
|
||||
ilj2 = -epsilon[itype][jtype] * n_x * pow(1 / rx, gamma[itype][jtype]) /
|
||||
(n_x - gamma[itype][jtype]);
|
||||
|
||||
evdwl = ilj1 + ilj2 - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
|
||||
}
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
if (vflag_fdotr) virial_fdotr_compute();
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::compute_inner()
|
||||
{
|
||||
int i, j, ii, jj, inum, jnum, itype, jtype;
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
|
||||
double rsq, r2inv, r6inv, forcelj, factor_lj, rsw;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
double r, rx, n_x;
|
||||
double pow_rx_n_x, pow_rx_gamma;
|
||||
double filj1, filj2, filj3, filj4, filj5, filj6, forceilj;
|
||||
double ilj1, ilj2;
|
||||
double fxtmp, fytmp, fztmp;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
inum = list->inum_inner;
|
||||
ilist = list->ilist_inner;
|
||||
numneigh = list->numneigh_inner;
|
||||
firstneigh = list->firstneigh_inner;
|
||||
|
||||
double cut_out_on = cut_respa[0];
|
||||
double cut_out_off = cut_respa[1];
|
||||
|
||||
double cut_out_diff = cut_out_off - cut_out_on;
|
||||
double cut_out_on_sq = cut_out_on * cut_out_on;
|
||||
double cut_out_off_sq = cut_out_off * cut_out_off;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
|
||||
if (rsq < cut_out_off_sq) {
|
||||
jtype = type[j];
|
||||
r = sqrt(rsq);
|
||||
|
||||
rx = r / rm[itype][jtype];
|
||||
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
|
||||
pow_rx_n_x = pow(1.0 / rx, n_x);
|
||||
pow_rx_gamma = pow(1.0 / rx, gamma[itype][jtype]);
|
||||
|
||||
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow_rx_n_x /
|
||||
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
|
||||
|
||||
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow_rx_gamma /
|
||||
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
|
||||
|
||||
filj3 = -2.0 * alpha[itype][jtype] * rx * pow_rx_gamma /
|
||||
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) *
|
||||
log(1 / rx) * pow_rx_n_x / (n_x - gamma[itype][jtype]);
|
||||
|
||||
filj5 = -1.0 * gamma[itype][jtype] * n_x * pow_rx_n_x / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj6 = +1.0 * gamma[itype][jtype] * n_x * pow_rx_gamma / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
// F = -dV/dr
|
||||
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
|
||||
fpair = factor_lj * forceilj / r; // F_x = -x/r * dV/dr (chain rule)
|
||||
|
||||
if (rsq > cut_out_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff;
|
||||
fpair *= 1.0 - rsw * rsw * (3.0 - 2.0 * rsw);
|
||||
}
|
||||
|
||||
fxtmp += delx * fpair;
|
||||
fytmp += dely * fpair;
|
||||
fztmp += delz * fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx * fpair;
|
||||
f[j][1] -= dely * fpair;
|
||||
f[j][2] -= delz * fpair;
|
||||
}
|
||||
}
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::compute_middle()
|
||||
|
||||
{
|
||||
int i, j, ii, jj, inum, jnum, itype, jtype;
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
|
||||
double rsq, r2inv, r6inv, forcelj, factor_lj, rsw;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
double r, rx, n_x;
|
||||
double pow_rx_n_x, pow_rx_gamma;
|
||||
double filj1, filj2, filj3, filj4, filj5, filj6, forceilj;
|
||||
double ilj1, ilj2;
|
||||
double fxtmp, fytmp, fztmp;
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
inum = list->inum_middle;
|
||||
ilist = list->ilist_middle;
|
||||
numneigh = list->numneigh_middle;
|
||||
firstneigh = list->firstneigh_middle;
|
||||
|
||||
double cut_in_off = cut_respa[0];
|
||||
double cut_in_on = cut_respa[1];
|
||||
double cut_out_on = cut_respa[2];
|
||||
double cut_out_off = cut_respa[3];
|
||||
|
||||
double cut_in_diff = cut_in_on - cut_in_off;
|
||||
double cut_out_diff = cut_out_off - cut_out_on;
|
||||
double cut_in_off_sq = cut_in_off * cut_in_off;
|
||||
double cut_in_on_sq = cut_in_on * cut_in_on;
|
||||
double cut_out_on_sq = cut_out_on * cut_out_on;
|
||||
double cut_out_off_sq = cut_out_off * cut_out_off;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
|
||||
if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
|
||||
jtype = type[j];
|
||||
r = sqrt(rsq);
|
||||
|
||||
rx = r / rm[itype][jtype];
|
||||
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
|
||||
pow_rx_n_x = pow(1.0 / rx, n_x);
|
||||
pow_rx_gamma = pow(1.0 / rx, gamma[itype][jtype]);
|
||||
|
||||
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow_rx_n_x /
|
||||
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
|
||||
|
||||
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow_rx_gamma /
|
||||
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
|
||||
|
||||
filj3 = -2.0 * alpha[itype][jtype] * rx * pow_rx_gamma /
|
||||
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) *
|
||||
log(1 / rx) * pow_rx_n_x / (n_x - gamma[itype][jtype]);
|
||||
|
||||
filj5 = -1.0 * gamma[itype][jtype] * n_x * pow_rx_n_x / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj6 = +1.0 * gamma[itype][jtype] * n_x * pow_rx_gamma / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
// F = -dV/dr
|
||||
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
|
||||
fpair = factor_lj * forceilj / r;
|
||||
if (rsq < cut_in_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff;
|
||||
fpair *= rsw * rsw * (3.0 - 2.0 * rsw);
|
||||
}
|
||||
if (rsq > cut_out_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff;
|
||||
fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0);
|
||||
}
|
||||
|
||||
fxtmp += delx * fpair;
|
||||
fytmp += dely * fpair;
|
||||
fztmp += delz * fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx * fpair;
|
||||
f[j][1] -= dely * fpair;
|
||||
f[j][2] -= delz * fpair;
|
||||
}
|
||||
}
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::compute_outer(int eflag, int vflag)
|
||||
|
||||
{
|
||||
int i, j, ii, jj, inum, jnum, itype, jtype;
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
|
||||
double rsq, r2inv, r6inv, forcelj, factor_lj, rsw;
|
||||
int *ilist, *jlist, *numneigh, **firstneigh;
|
||||
|
||||
double r, rx, n_x;
|
||||
double pow_rx_n_x, pow_rx_gamma;
|
||||
double filj1, filj2, filj3, filj4, filj5, filj6, forceilj;
|
||||
double ilj1, ilj2;
|
||||
double fxtmp, fytmp, fztmp;
|
||||
|
||||
evdwl = 0.0;
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
int *type = atom->type;
|
||||
int nlocal = atom->nlocal;
|
||||
double *special_lj = force->special_lj;
|
||||
int newton_pair = force->newton_pair;
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
numneigh = list->numneigh;
|
||||
firstneigh = list->firstneigh;
|
||||
|
||||
double cut_in_off = cut_respa[2];
|
||||
double cut_in_on = cut_respa[3];
|
||||
|
||||
double cut_in_diff = cut_in_on - cut_in_off;
|
||||
double cut_in_off_sq = cut_in_off * cut_in_off;
|
||||
double cut_in_on_sq = cut_in_on * cut_in_on;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (ii = 0; ii < inum; ii++) {
|
||||
i = ilist[ii];
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j][0];
|
||||
dely = ytmp - x[j][1];
|
||||
delz = ztmp - x[j][2];
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsq[itype][jtype]) {
|
||||
if (rsq > cut_in_off_sq) {
|
||||
r = sqrt(rsq);
|
||||
|
||||
rx = r / rm[itype][jtype];
|
||||
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
|
||||
pow_rx_n_x = pow(1.0 / rx, n_x);
|
||||
pow_rx_gamma = pow(1.0 / rx, gamma[itype][jtype]);
|
||||
|
||||
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow_rx_n_x /
|
||||
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
|
||||
|
||||
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow_rx_gamma /
|
||||
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
|
||||
|
||||
filj3 = -2.0 * alpha[itype][jtype] * rx * pow_rx_gamma /
|
||||
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) *
|
||||
log(1 / rx) * pow_rx_n_x / (n_x - gamma[itype][jtype]);
|
||||
|
||||
filj5 = -1.0 * gamma[itype][jtype] * n_x * pow_rx_n_x / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj6 =
|
||||
+1.0 * gamma[itype][jtype] * n_x * pow_rx_gamma / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
// F = -dV/dr
|
||||
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
|
||||
fpair = factor_lj * forceilj / r;
|
||||
if (rsq < cut_in_on_sq) {
|
||||
rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff;
|
||||
fpair *= rsw * rsw * (3.0 - 2.0 * rsw);
|
||||
}
|
||||
|
||||
fxtmp += delx * fpair;
|
||||
fytmp += dely * fpair;
|
||||
fztmp += delz * fpair;
|
||||
if (newton_pair || j < nlocal) {
|
||||
f[j][0] -= delx * fpair;
|
||||
f[j][1] -= dely * fpair;
|
||||
f[j][2] -= delz * fpair;
|
||||
}
|
||||
}
|
||||
|
||||
if (eflag) {
|
||||
|
||||
r = sqrt(rsq);
|
||||
|
||||
rx = r / rm[itype][jtype];
|
||||
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
|
||||
|
||||
ilj1 = epsilon[itype][jtype] * gamma[itype][jtype] * pow(1 / rx, n_x) /
|
||||
(n_x - gamma[itype][jtype]);
|
||||
ilj2 = -epsilon[itype][jtype] * n_x * pow(1 / rx, gamma[itype][jtype]) /
|
||||
(n_x - gamma[itype][jtype]);
|
||||
|
||||
evdwl = ilj1 + ilj2 - offset[itype][jtype];
|
||||
evdwl *= factor_lj;
|
||||
}
|
||||
|
||||
if (vflag) {
|
||||
if (rsq <= cut_in_off_sq) {
|
||||
|
||||
r = sqrt(rsq);
|
||||
|
||||
rx = r / rm[itype][jtype];
|
||||
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
|
||||
pow_rx_n_x = pow(1.0 / rx, n_x);
|
||||
pow_rx_gamma = pow(1.0 / rx, gamma[itype][jtype]);
|
||||
|
||||
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow_rx_n_x /
|
||||
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
|
||||
|
||||
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow_rx_gamma /
|
||||
(square(n_x - gamma[itype][jtype]) * rm[itype][jtype]);
|
||||
|
||||
filj3 = -2.0 * alpha[itype][jtype] * rx * pow_rx_gamma /
|
||||
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) *
|
||||
log(1 / rx) * pow_rx_n_x / (n_x - gamma[itype][jtype]);
|
||||
|
||||
filj5 =
|
||||
-1.0 * gamma[itype][jtype] * n_x * pow_rx_n_x / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj6 =
|
||||
+1.0 * gamma[itype][jtype] * n_x * pow_rx_gamma / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
// F = -dV/dr
|
||||
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
|
||||
fpair = factor_lj * forceilj / r;
|
||||
|
||||
} else if (rsq < cut_in_on_sq)
|
||||
fpair = factor_lj * forceilj / r;
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
|
||||
}
|
||||
}
|
||||
f[i][0] += fxtmp;
|
||||
f[i][1] += fytmp;
|
||||
f[i][2] += fztmp;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
global settings
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 1)
|
||||
error->all(FLERR, "Pair style ilj/cut must have exactly one argument: cutoff distance");
|
||||
|
||||
cut_global = utils::numeric(FLERR, arg[0], false, lmp);
|
||||
|
||||
// reset cutoffs that have been explicitly set
|
||||
|
||||
if (allocated) {
|
||||
int i, j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++)
|
||||
if (setflag[i][j]) cut[i][j] = cut_global;
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
set coeffs for one or more type pairs
|
||||
------------------------------------------------------------------------- */
|
||||
/*
|
||||
|
||||
7 or 8 coefficients: 5 for the ILJ, 2 for the pair, 1 for the cutoff (optional)
|
||||
|
||||
*/
|
||||
void PairLJPirani::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg < 7 || narg > 8) error->all(FLERR, "Incorrect args for pair coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo, ihi, jlo, jhi;
|
||||
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
|
||||
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
|
||||
|
||||
double alpha_one = utils::numeric(FLERR, arg[2], false, lmp);
|
||||
double beta_one = utils::numeric(FLERR, arg[3], false, lmp);
|
||||
double gamma_one = utils::numeric(FLERR, arg[4], false, lmp);
|
||||
double rm_one = utils::numeric(FLERR, arg[5], false, lmp);
|
||||
double epsilon_one = utils::numeric(FLERR, arg[6], false, lmp);
|
||||
|
||||
double cut_one = cut_global;
|
||||
if (narg == 8) cut_one = utils::numeric(FLERR, arg[7], false, lmp);
|
||||
|
||||
if (rm_one <= 0.0 || epsilon_one < 0.0 || gamma_one <= 0.0)
|
||||
error->all(FLERR, "Illegal ILJ coefficients");
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo, i); j <= jhi; j++) {
|
||||
alpha[i][j] = alpha_one;
|
||||
beta[i][j] = beta_one;
|
||||
gamma[i][j] = gamma_one;
|
||||
rm[i][j] = rm_one;
|
||||
epsilon[i][j] = epsilon_one;
|
||||
cut[i][j] = cut_one;
|
||||
setflag[i][j] = 1;
|
||||
count++;
|
||||
}
|
||||
}
|
||||
|
||||
// Initialize symmetric entries
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
for (int j = MAX(jlo, i); j <= jhi; j++) {
|
||||
alpha[j][i] = alpha[i][j];
|
||||
beta[j][i] = beta[i][j];
|
||||
gamma[j][i] = gamma[i][j];
|
||||
rm[j][i] = rm[i][j];
|
||||
epsilon[j][i] = epsilon[i][j];
|
||||
cut[j][i] = cut[i][j];
|
||||
setflag[j][i] = setflag[i][j];
|
||||
}
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init specific to this pair style
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::init_style()
|
||||
{
|
||||
// request regular or rRESPA neighbor list
|
||||
|
||||
int list_style = NeighConst::REQ_DEFAULT;
|
||||
|
||||
if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) {
|
||||
auto respa = dynamic_cast<Respa *>(update->integrate);
|
||||
if (respa->level_inner >= 0) list_style = NeighConst::REQ_RESPA_INOUT;
|
||||
if (respa->level_middle >= 0) list_style = NeighConst::REQ_RESPA_ALL;
|
||||
}
|
||||
neighbor->add_request(this, list_style);
|
||||
|
||||
// set rRESPA cutoffs
|
||||
|
||||
if (utils::strmatch(update->integrate_style, "^respa") &&
|
||||
(dynamic_cast<Respa *>(update->integrate))->level_inner >= 0)
|
||||
cut_respa = (dynamic_cast<Respa *>(update->integrate))->cutoff;
|
||||
else
|
||||
cut_respa = nullptr;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
init for one type pair i,j and corresponding j,i
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
double PairLJPirani::init_one(int i, int j)
|
||||
{
|
||||
if (setflag[i][j] == 0) error->all(FLERR, "All pair coeffs are not set");
|
||||
|
||||
if (offset_flag && (cut[i][j] > 0.0)) {
|
||||
double r = cut[i][j] / rm[i][j];
|
||||
double nx = alpha[i][j] * r * r + beta[i][j];
|
||||
offset[i][j] = epsilon[i][j] *
|
||||
((gamma[i][j] / (nx - gamma[i][j])) * pow(1 / r, nx) -
|
||||
(nx / (nx - gamma[i][j])) * pow(1 / r, gamma[i][j]));
|
||||
} else
|
||||
offset[i][j] = 0.0;
|
||||
|
||||
alpha[j][i] = alpha[i][j];
|
||||
beta[j][i] = beta[i][j];
|
||||
gamma[j][i] = gamma[i][j];
|
||||
rm[j][i] = rm[i][j];
|
||||
epsilon[j][i] = epsilon[i][j];
|
||||
offset[j][i] = offset[i][j];
|
||||
|
||||
// check interior rRESPA cutoff
|
||||
|
||||
if (cut_respa && cut[i][j] < cut_respa[3])
|
||||
error->all(FLERR, "Pair cutoff < Respa interior cutoff");
|
||||
|
||||
return cut[i][j];
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::write_restart(FILE *fp)
|
||||
{
|
||||
write_restart_settings(fp);
|
||||
|
||||
int i, j;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
fwrite(&setflag[i][j], sizeof(int), 1, fp);
|
||||
if (setflag[i][j]) {
|
||||
fwrite(&alpha[i][j], sizeof(double), 1, fp);
|
||||
fwrite(&beta[i][j], sizeof(double), 1, fp);
|
||||
fwrite(&gamma[i][j], sizeof(double), 1, fp);
|
||||
fwrite(&rm[i][j], sizeof(double), 1, fp);
|
||||
fwrite(&epsilon[i][j], sizeof(double), 1, fp);
|
||||
fwrite(&cut[i][j], sizeof(double), 1, fp);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to restart file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&cut_global, sizeof(double), 1, fp);
|
||||
fwrite(&offset_flag, sizeof(int), 1, fp);
|
||||
fwrite(&mix_flag, sizeof(int), 1, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::read_restart(FILE *fp)
|
||||
{
|
||||
read_restart_settings(fp);
|
||||
allocate();
|
||||
|
||||
int i, j;
|
||||
int me = comm->me;
|
||||
for (i = 1; i <= atom->ntypes; i++)
|
||||
for (j = i; j <= atom->ntypes; j++) {
|
||||
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
|
||||
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
|
||||
if (setflag[i][j]) {
|
||||
if (me == 0) {
|
||||
utils::sfread(FLERR, &alpha[i][j], sizeof(double), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &beta[i][j], sizeof(double), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &gamma[i][j], sizeof(double), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &rm[i][j], sizeof(double), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&alpha[i][j], 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&beta[i][j], 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&gamma[i][j], 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&rm[i][j], 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 reads from restart file, bcasts
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::read_restart_settings(FILE *fp)
|
||||
{
|
||||
int me = comm->me;
|
||||
if (me == 0) {
|
||||
utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
|
||||
MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
fprintf(fp, "%d %g %g %g %g %g\n", i, alpha[i][i], beta[i][i], gamma[i][i], rm[i][i],
|
||||
epsilon[i][i]);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
proc 0 writes all pairs to data file
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPirani::write_data_all(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->ntypes; i++)
|
||||
for (int j = i; j <= atom->ntypes; j++)
|
||||
fprintf(fp, "%d %d %g %g %g %g %g %g\n", i, j, alpha[i][j], beta[i][j], gamma[i][j], rm[i][j],
|
||||
epsilon[i][j], cut[i][j]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJPirani::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
|
||||
double /*factor_coul*/, double factor_lj, double &fforce)
|
||||
{
|
||||
double r, rx, n_x, filj1, filj2, filj3, filj4, filj5, filj6, forceilj;
|
||||
double ilj1, ilj2;
|
||||
|
||||
r = sqrt(rsq);
|
||||
rx = r / rm[itype][jtype];
|
||||
n_x = alpha[itype][jtype] * rx * rx + beta[itype][jtype];
|
||||
filj1 = -2.0 * alpha[itype][jtype] * gamma[itype][jtype] * rx * pow(1 / rx, n_x) /
|
||||
(pow(n_x - gamma[itype][jtype], 2.0) * rm[itype][jtype]);
|
||||
|
||||
filj2 = +2.0 * alpha[itype][jtype] * rx * n_x * pow(1 / rx, gamma[itype][jtype]) /
|
||||
(pow(n_x - gamma[itype][jtype], 2.0) * rm[itype][jtype]);
|
||||
|
||||
filj3 = -2.0 * alpha[itype][jtype] * rx * pow(1 / rx, gamma[itype][jtype]) /
|
||||
(rm[itype][jtype] * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj4 = +2.0 * alpha[itype][jtype] * gamma[itype][jtype] * (rx / rm[itype][jtype]) * log(1 / rx) *
|
||||
pow(1 / rx, n_x) / (n_x - gamma[itype][jtype]);
|
||||
|
||||
filj5 = -1.0 * gamma[itype][jtype] * n_x * pow(1 / rx, n_x) / (r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
filj6 = +1.0 * gamma[itype][jtype] * n_x * pow(1 / rx, gamma[itype][jtype]) /
|
||||
(r * (n_x - gamma[itype][jtype]));
|
||||
|
||||
forceilj = -epsilon[itype][jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
|
||||
|
||||
fforce = factor_lj * forceilj / r;
|
||||
|
||||
ilj1 =
|
||||
epsilon[itype][jtype] * gamma[itype][jtype] * pow(1 / rx, n_x) / (n_x - gamma[itype][jtype]);
|
||||
ilj2 =
|
||||
-epsilon[itype][jtype] * n_x * pow(1 / rx, gamma[itype][jtype]) / (n_x - gamma[itype][jtype]);
|
||||
return factor_lj * (ilj1 + ilj2 - offset[itype][jtype]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void *PairLJPirani::extract(const char *str, int &dim)
|
||||
{
|
||||
dim = 2;
|
||||
if (strcmp(str, "alpha") == 0) return (void *) alpha;
|
||||
if (strcmp(str, "beta") == 0) return (void *) beta;
|
||||
if (strcmp(str, "gamma") == 0) return (void *) gamma;
|
||||
if (strcmp(str, "rm") == 0) return (void *) rm;
|
||||
if (strcmp(str, "epsilon") == 0) return (void *) epsilon;
|
||||
return nullptr;
|
||||
}
|
||||
60
src/EXTRA-PAIR/pair_lj_pirani.h
Normal file
60
src/EXTRA-PAIR/pair_lj_pirani.h
Normal file
@ -0,0 +1,60 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
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 PAIR_CLASS
|
||||
// clang-format off
|
||||
PairStyle(lj/pirani,PairLJPirani);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_LJ_PIRANI_H
|
||||
#define LMP_PAIR_LJ_PIRANI_H
|
||||
|
||||
#include "pair.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairLJPirani : public Pair {
|
||||
public:
|
||||
PairLJPirani(class LAMMPS *);
|
||||
virtual ~PairLJPirani() override;
|
||||
void compute(int, int) override;
|
||||
void settings(int, char **) override;
|
||||
void coeff(int, char **) override;
|
||||
void init_style() override;
|
||||
double init_one(int, int) override;
|
||||
void write_restart(FILE *) override;
|
||||
void read_restart(FILE *) override;
|
||||
void write_restart_settings(FILE *) override;
|
||||
void read_restart_settings(FILE *) override;
|
||||
void write_data(FILE *) override;
|
||||
void write_data_all(FILE *) override;
|
||||
double single(int, int, int, int, double, double, double, double &) override;
|
||||
void *extract(const char *, int &) override;
|
||||
void compute_inner() override;
|
||||
void compute_middle() override;
|
||||
void compute_outer(int, int) override;
|
||||
|
||||
protected:
|
||||
double cut_global;
|
||||
double **cut;
|
||||
double **alpha, **beta, **gamma, **rm, **epsilon;
|
||||
double **offset;
|
||||
double *cut_respa;
|
||||
virtual void allocate();
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
194
src/OPENMP/pair_lj_pirani_omp.cpp
Normal file
194
src/OPENMP/pair_lj_pirani_omp.cpp
Normal file
@ -0,0 +1,194 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
This software is distributed under the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
Contributing author: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "pair_lj_pirani_omp.h"
|
||||
|
||||
#include "atom.h"
|
||||
#include "comm.h"
|
||||
#include "force.h"
|
||||
#include "math_special.h"
|
||||
#include "neigh_list.h"
|
||||
#include "suffix.h"
|
||||
|
||||
#include "omp_compat.h"
|
||||
using namespace LAMMPS_NS;
|
||||
using MathSpecial::square;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
PairLJPiraniOMP::PairLJPiraniOMP(LAMMPS *lmp) : PairLJPirani(lmp), ThrOMP(lmp, THR_PAIR)
|
||||
{
|
||||
suffix_flag |= Suffix::OMP;
|
||||
respa_enable = 0;
|
||||
cut_respa = nullptr;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairLJPiraniOMP::compute(int eflag, int vflag)
|
||||
{
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
const int nall = atom->nlocal + atom->nghost;
|
||||
const int nthreads = comm->nthreads;
|
||||
const int inum = list->inum;
|
||||
|
||||
#if defined(_OPENMP)
|
||||
#pragma omp parallel LMP_DEFAULT_NONE LMP_SHARED(eflag, vflag)
|
||||
#endif
|
||||
{
|
||||
int ifrom, ito, tid;
|
||||
|
||||
loop_setup_thr(ifrom, ito, tid, inum, nthreads);
|
||||
ThrData *thr = fix->get_thr(tid);
|
||||
thr->timer(Timer::START);
|
||||
ev_setup_thr(eflag, vflag, nall, eatom, vatom, nullptr, thr);
|
||||
|
||||
if (evflag) {
|
||||
if (eflag) {
|
||||
if (force->newton_pair)
|
||||
eval<1, 1, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<1, 1, 0>(ifrom, ito, thr);
|
||||
} else {
|
||||
if (force->newton_pair)
|
||||
eval<1, 0, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<1, 0, 0>(ifrom, ito, thr);
|
||||
}
|
||||
} else {
|
||||
if (force->newton_pair)
|
||||
eval<0, 0, 1>(ifrom, ito, thr);
|
||||
else
|
||||
eval<0, 0, 0>(ifrom, ito, thr);
|
||||
}
|
||||
thr->timer(Timer::PAIR);
|
||||
reduce_thr(this, eflag, vflag, thr);
|
||||
} // end of omp parallel region
|
||||
}
|
||||
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
void PairLJPiraniOMP::eval(int iifrom, int iito, ThrData *const thr)
|
||||
{
|
||||
const auto *_noalias const x = (dbl3_t *) atom->x[0];
|
||||
auto *_noalias const f = (dbl3_t *) thr->get_f()[0];
|
||||
const int *_noalias const type = atom->type;
|
||||
const double *_noalias const special_lj = force->special_lj;
|
||||
const int *_noalias const ilist = list->ilist;
|
||||
const int *_noalias const numneigh = list->numneigh;
|
||||
const int *const *const firstneigh = list->firstneigh;
|
||||
|
||||
double xtmp, ytmp, ztmp, delx, dely, delz, fxtmp, fytmp, fztmp;
|
||||
double rsq, forceilj, factor_lj, evdwl, fpair;
|
||||
|
||||
const int nlocal = atom->nlocal;
|
||||
int j, jj, jnum, jtype;
|
||||
|
||||
evdwl = 0.0;
|
||||
|
||||
// loop over neighbors of my atoms
|
||||
|
||||
for (int ii = iifrom; ii < iito; ++ii) {
|
||||
const int i = ilist[ii];
|
||||
const int itype = type[i];
|
||||
const int *_noalias const jlist = firstneigh[i];
|
||||
const double *_noalias const cutsqi = cutsq[itype];
|
||||
const double *_noalias const offseti = offset[itype];
|
||||
const double *_noalias const alphai = alpha[itype];
|
||||
const double *_noalias const betai = beta[itype];
|
||||
const double *_noalias const gammai = gamma[itype];
|
||||
const double *_noalias const rmi = rm[itype];
|
||||
const double *_noalias const epsiloni = epsilon[itype];
|
||||
|
||||
xtmp = x[i].x;
|
||||
ytmp = x[i].y;
|
||||
ztmp = x[i].z;
|
||||
jnum = numneigh[i];
|
||||
fxtmp = fytmp = fztmp = 0.0;
|
||||
|
||||
for (jj = 0; jj < jnum; jj++) {
|
||||
j = jlist[jj];
|
||||
factor_lj = special_lj[sbmask(j)];
|
||||
j &= NEIGHMASK;
|
||||
|
||||
delx = xtmp - x[j].x;
|
||||
dely = ytmp - x[j].y;
|
||||
delz = ztmp - x[j].z;
|
||||
rsq = delx * delx + dely * dely + delz * delz;
|
||||
jtype = type[j];
|
||||
|
||||
if (rsq < cutsqi[jtype]) {
|
||||
const double r = sqrt(rsq);
|
||||
|
||||
const double rx = r / rmi[jtype];
|
||||
const double n_x = alphai[jtype] * rx * rx + betai[jtype];
|
||||
const double pow_rx_n_x = pow(1.0 / rx, n_x);
|
||||
const double pow_rx_gamma = pow(1.0 / rx, gammai[jtype]);
|
||||
|
||||
double filj1 = -2.0 * alphai[jtype] * gammai[jtype] * rx * pow_rx_n_x /
|
||||
(square(n_x - gammai[jtype]) * rmi[jtype]);
|
||||
|
||||
double filj2 = 2.0 * alphai[jtype] * rx * n_x * pow_rx_gamma /
|
||||
(square(n_x - gammai[jtype]) * rmi[jtype]);
|
||||
|
||||
double filj3 =
|
||||
-2.0 * alphai[jtype] * rx * pow_rx_gamma / (rmi[jtype] * (n_x - gammai[jtype]));
|
||||
|
||||
double filj4 = 2.0 * alphai[jtype] * gammai[jtype] * (rx / rmi[jtype]) * log(1 / rx) *
|
||||
pow_rx_n_x / (n_x - gammai[jtype]);
|
||||
|
||||
double filj5 = -1.0 * gammai[jtype] * n_x * pow_rx_n_x / (r * (n_x - gammai[jtype]));
|
||||
|
||||
double filj6 = 1.0 * gammai[jtype] * n_x * pow_rx_gamma / (r * (n_x - gammai[jtype]));
|
||||
|
||||
// F = -dV/dr
|
||||
forceilj = -epsiloni[jtype] * (filj1 + filj2 + filj3 + filj4 + filj5 + filj6);
|
||||
fpair = factor_lj * forceilj / r; // F_x = -x/r * dV/dr (chain rule)
|
||||
|
||||
fxtmp += delx * fpair;
|
||||
fytmp += dely * fpair;
|
||||
fztmp += delz * fpair;
|
||||
if (NEWTON_PAIR || j < nlocal) {
|
||||
f[j].x -= delx * fpair;
|
||||
f[j].y -= dely * fpair;
|
||||
f[j].z -= delz * fpair;
|
||||
}
|
||||
|
||||
if (EFLAG) {
|
||||
double ilj1 = epsiloni[jtype] * gammai[jtype] * pow(1 / rx, n_x) / (n_x - gammai[jtype]);
|
||||
double ilj2 = -epsiloni[jtype] * n_x * pow(1 / rx, gammai[jtype]) / (n_x - gammai[jtype]);
|
||||
|
||||
evdwl = ilj1 + ilj2 - offseti[jtype];
|
||||
evdwl *= factor_lj;
|
||||
}
|
||||
|
||||
if (EVFLAG)
|
||||
ev_tally_thr(this, i, j, nlocal, NEWTON_PAIR, evdwl, 0.0, fpair, delx, dely, delz, thr);
|
||||
}
|
||||
}
|
||||
f[i].x += fxtmp;
|
||||
f[i].y += fytmp;
|
||||
f[i].z += fztmp;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairLJPiraniOMP::memory_usage()
|
||||
{
|
||||
double bytes = memory_usage_thr();
|
||||
bytes += PairLJPirani::memory_usage();
|
||||
|
||||
return bytes;
|
||||
}
|
||||
48
src/OPENMP/pair_lj_pirani_omp.h
Normal file
48
src/OPENMP/pair_lj_pirani_omp.h
Normal file
@ -0,0 +1,48 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
LAMMPS development team: developers@lammps.org
|
||||
|
||||
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: Axel Kohlmeyer (Temple U)
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifdef PAIR_CLASS
|
||||
// clang-format off
|
||||
PairStyle(lj/pirani,PairLJPirani);
|
||||
// clang-format on
|
||||
#else
|
||||
|
||||
#ifndef LMP_PAIR_LJ_PIRANI_OMP_H
|
||||
#define LMP_PAIR_LJ_PIRANI_OMP_H
|
||||
|
||||
#include "pair_lj_pirani.h"
|
||||
#include "thr_omp.h"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
|
||||
class PairLJPiraniOMP : public PairLJPirani, public ThrOMP {
|
||||
|
||||
public:
|
||||
PairLJPiraniOMP(class LAMMPS *);
|
||||
|
||||
void compute(int, int) override;
|
||||
double memory_usage() override;
|
||||
|
||||
private:
|
||||
template <int EVFLAG, int EFLAG, int NEWTON_PAIR>
|
||||
void eval(int ifrom, int ito, ThrData *const thr);
|
||||
};
|
||||
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
#endif
|
||||
@ -101,6 +101,6 @@ void sfree(void *ptr)
|
||||
if (ptr == nullptr) return;
|
||||
|
||||
free(ptr);
|
||||
ptr = nullptr;
|
||||
ptr = nullptr; // NOTE: this has no effect
|
||||
}
|
||||
} // namespace ReaxFF
|
||||
|
||||
@ -44,7 +44,7 @@ using namespace RHEO_NS;
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) :
|
||||
BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr), nbond(nullptr),
|
||||
BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr), dbond(nullptr),
|
||||
id_fix(nullptr), compute_surface(nullptr)
|
||||
{
|
||||
partial_flag = 1;
|
||||
@ -71,7 +71,6 @@ BondRHEOShell::BondRHEOShell(LAMMPS *_lmp) :
|
||||
modify->add_fix(fmt::format("{} all property/atom i_shell_nbond", id_fix));
|
||||
index_nb = atom->find_custom("shell_nbond", tmp1, tmp2);
|
||||
}
|
||||
nbond = atom->ivector[index_nb];
|
||||
|
||||
//Store non-persistent per atom quantities, intermediate
|
||||
|
||||
@ -181,6 +180,7 @@ void BondRHEOShell::compute(int eflag, int vflag)
|
||||
double **v = atom->v;
|
||||
double **f = atom->f;
|
||||
tagint *tag = atom->tag;
|
||||
int *nbond = atom->ivector[index_nb];
|
||||
int *status = atom->rheo_status;
|
||||
int **bondlist = neighbor->bondlist;
|
||||
int nbondlist = neighbor->nbondlist;
|
||||
|
||||
@ -45,7 +45,7 @@ class BondRHEOShell : public BondBPM {
|
||||
double *k, *ecrit, *gamma;
|
||||
double tform, rmax, rsurf;
|
||||
|
||||
int *dbond, *nbond;
|
||||
int *dbond;
|
||||
int index_nb, nmax_store;
|
||||
char *id_fix;
|
||||
|
||||
|
||||
@ -39,12 +39,10 @@ using namespace LAMMPS_NS;
|
||||
using namespace RHEO_NS;
|
||||
using namespace MathExtra;
|
||||
|
||||
static constexpr double EPSILON = 1e-1;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
|
||||
Compute(lmp, narg, arg), chi(nullptr), fp_store(nullptr), fix_rheo(nullptr), rho0(nullptr),
|
||||
Compute(lmp, narg, arg), chi(nullptr), fix_rheo(nullptr), rho0(nullptr),
|
||||
norm(nullptr), normwf(nullptr), id_fix_pa(nullptr), list(nullptr), compute_kernel(nullptr),
|
||||
fix_pressure(nullptr)
|
||||
|
||||
@ -64,13 +62,12 @@ ComputeRHEOInterface::ComputeRHEOInterface(LAMMPS *lmp, int narg, char **arg) :
|
||||
// between timesteps (fix property atom will handle callbacks)
|
||||
|
||||
int tmp1, tmp2;
|
||||
int index = atom->find_custom("fp_store", tmp1, tmp2);
|
||||
if (index == -1) {
|
||||
index_fp_store = atom->find_custom("fp_store", tmp1, tmp2);
|
||||
if (index_fp_store == -1) {
|
||||
id_fix_pa = utils::strdup(id + std::string("_fix_property_atom"));
|
||||
modify->add_fix(fmt::format("{} all property/atom d2_fp_store 3", id_fix_pa));
|
||||
index = atom->find_custom("fp_store", tmp1, tmp2);
|
||||
index_fp_store = atom->find_custom("fp_store", tmp1, tmp2);
|
||||
}
|
||||
fp_store = atom->darray[index];
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -114,7 +111,7 @@ void ComputeRHEOInterface::init_list(int /*id*/, NeighList *ptr)
|
||||
|
||||
void ComputeRHEOInterface::compute_peratom()
|
||||
{
|
||||
int a, i, j, ii, jj, jnum, itype, fluidi, fluidj, status_match;
|
||||
int a, i, j, ii, jj, jnum, fluidi, fluidj, status_match;
|
||||
double xtmp, ytmp, ztmp, rsq, w, dot, dx[3];
|
||||
|
||||
int inum, *ilist, *jlist, *numneigh, **firstneigh;
|
||||
@ -126,6 +123,7 @@ void ComputeRHEOInterface::compute_peratom()
|
||||
int newton = force->newton;
|
||||
int *status = atom->rheo_status;
|
||||
double *rho = atom->rho;
|
||||
double **fp_store = atom->darray[index_fp_store];
|
||||
|
||||
inum = list->inum;
|
||||
ilist = list->ilist;
|
||||
@ -151,7 +149,6 @@ void ComputeRHEOInterface::compute_peratom()
|
||||
xtmp = x[i][0];
|
||||
ytmp = x[i][1];
|
||||
ztmp = x[i][2];
|
||||
itype = type[i];
|
||||
fluidi = !(status[i] & PHASECHECK);
|
||||
jlist = firstneigh[i];
|
||||
jnum = numneigh[i];
|
||||
@ -212,9 +209,10 @@ void ComputeRHEOInterface::compute_peratom()
|
||||
if (status[i] & PHASECHECK) {
|
||||
if (normwf[i] != 0.0) {
|
||||
// Stores rho for solid particles 1+Pw in Adami Adams 2012
|
||||
rho[i] = MAX(EPSILON, fix_pressure->calc_rho(rho[i] / normwf[i], i));
|
||||
// cap out at a tenth of equilibrium
|
||||
rho[i] = MAX(0.1 * rho0[type[i]], fix_pressure->calc_rho(rho[i] / normwf[i], i));
|
||||
} else {
|
||||
rho[i] = rho0[itype];
|
||||
rho[i] = rho0[type[i]];
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -230,6 +228,7 @@ int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /
|
||||
{
|
||||
int m = 0;
|
||||
double *rho = atom->rho;
|
||||
double **fp_store = atom->darray[index_fp_store];
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
int j = list[i];
|
||||
@ -250,6 +249,8 @@ int ComputeRHEOInterface::pack_forward_comm(int n, int *list, double *buf, int /
|
||||
void ComputeRHEOInterface::unpack_forward_comm(int n, int first, double *buf)
|
||||
{
|
||||
double *rho = atom->rho;
|
||||
double **fp_store = atom->darray[index_fp_store];
|
||||
|
||||
int m = 0;
|
||||
int last = first + n;
|
||||
for (int i = first; i < last; i++) {
|
||||
@ -335,6 +336,7 @@ void ComputeRHEOInterface::store_forces()
|
||||
double *mass = atom->mass;
|
||||
double *rmass = atom->rmass;
|
||||
double **f = atom->f;
|
||||
double **fp_store = atom->darray[index_fp_store];
|
||||
|
||||
// When this is called, fp_store stores the pressure force
|
||||
// After this method, fp_store instead stores non-pressure forces
|
||||
|
||||
@ -40,7 +40,8 @@ class ComputeRHEOInterface : public Compute {
|
||||
double correct_rho(int);
|
||||
void store_forces();
|
||||
|
||||
double *chi, **fp_store;
|
||||
double *chi;
|
||||
int index_fp_store;
|
||||
class FixRHEO *fix_rheo;
|
||||
|
||||
private:
|
||||
|
||||
@ -220,6 +220,8 @@ double ComputeRHEOKernel::calc_dw(int i, int j, double delx, double dely, double
|
||||
return calc_dw_wendlandc4(delx, dely, delz, r, dWij, dWji);
|
||||
if (kernel_style == QUINTIC)
|
||||
return calc_dw_quintic(delx, dely, delz, r, dWij, dWji);
|
||||
if (kernel_style == RK0)
|
||||
return calc_dw_quintic(delx, dely, delz, r, dWij, dWji);
|
||||
|
||||
double wp;
|
||||
int corrections_i = check_corrections(i);
|
||||
|
||||
@ -502,7 +502,7 @@ void ComputeRHEOPropertyAtom::pack_total_stress(int n)
|
||||
|
||||
void ComputeRHEOPropertyAtom::pack_nbond_shell(int n)
|
||||
{
|
||||
int *nbond = fix_oxidation->nbond;
|
||||
int *nbond = atom->ivector[fix_oxidation->index_nb];
|
||||
int *mask = atom->mask;
|
||||
int nlocal = atom->nlocal;
|
||||
|
||||
|
||||
@ -39,6 +39,7 @@ using namespace RHEO_NS;
|
||||
using namespace FixConst;
|
||||
|
||||
static const char cite_rheo[] =
|
||||
"RHEO package: doi:10.1063/5.0228823\n\n"
|
||||
"@article{Palermo2024,\n"
|
||||
" journal = {Physics of Fluids},\n"
|
||||
" title = {Reproducing hydrodynamics and elastic objects: A hybrid mesh-free model framework for dynamic multi-phase flows},\n"
|
||||
@ -117,17 +118,17 @@ FixRHEO::FixRHEO(LAMMPS *lmp, int narg, char **arg) :
|
||||
shift_flag = 1;
|
||||
memory->create(shift_type, n + 1, "rheo:shift_type");
|
||||
for (i = 1; i <= n; i++) shift_type[i] = 1;
|
||||
while (iarg < narg) { // optional sub-arguments
|
||||
if (strcmp(arg[iarg], "scale/cross/type") == 0) {
|
||||
if (iarg + 3 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error);
|
||||
while (iarg + 1 < narg) { // optional sub-arguments
|
||||
if (strcmp(arg[iarg + 1], "scale/cross/type") == 0) {
|
||||
if (iarg + 4 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift scale/cross/type", error);
|
||||
shift_cross_type_flag = 1;
|
||||
shift_scale = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
|
||||
shift_cmin = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
|
||||
shift_wmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
|
||||
shift_scale = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
|
||||
shift_cmin = utils::numeric(FLERR, arg[iarg + 3], false, lmp);
|
||||
shift_wmin = utils::numeric(FLERR, arg[iarg + 4], false, lmp);
|
||||
iarg += 3;
|
||||
} else if (strcmp(arg[iarg], "exclude/type") == 0) {
|
||||
if (iarg + 1 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error);
|
||||
utils::bounds(FLERR, arg[iarg + 1], 1, n, nlo, nhi, error);
|
||||
} else if (strcmp(arg[iarg + 1], "exclude/type") == 0) {
|
||||
if (iarg + 2 >= narg) utils::missing_cmd_args(FLERR, "fix rheo shift exclude/type", error);
|
||||
utils::bounds(FLERR, arg[iarg + 2], 1, n, nlo, nhi, error);
|
||||
for (i = nlo; i <= nhi; i++) shift_type[i] = 0;
|
||||
iarg += 1;
|
||||
} else {
|
||||
|
||||
@ -38,6 +38,7 @@ using namespace FixConst;
|
||||
enum { NONE, CONSTANT };
|
||||
|
||||
static const char cite_rheo_oxide[] =
|
||||
"RHEO oxidation: doi:10.1016/j.apm.2024.02.027\n\n"
|
||||
"@article{ApplMathModel.130.310,\n"
|
||||
" title = {A hybrid smoothed-particle hydrodynamics model of oxide skins on molten aluminum},\n"
|
||||
" journal = {Applied Mathematical Modelling},\n"
|
||||
@ -109,7 +110,6 @@ void FixRHEOOxidation::init()
|
||||
int tmp1, tmp2;
|
||||
index_nb = atom->find_custom("shell_nbond", tmp1, tmp2);
|
||||
if (index_nb == -1) error->all(FLERR, "Must use bond style rheo/shell to use fix rheo/oxidation");
|
||||
nbond = atom->ivector[index_nb];
|
||||
|
||||
// need a half neighbor list
|
||||
auto req = neighbor->add_request(this, NeighConst::REQ_FULL);
|
||||
|
||||
@ -39,11 +39,11 @@ class FixRHEOOxidation : public Fix {
|
||||
void post_force(int) override;
|
||||
int pack_forward_comm(int, int *, double *, int, int *) override;
|
||||
void unpack_forward_comm(int, int, double *) override;
|
||||
int *nbond;
|
||||
double rsurf, cut;
|
||||
int index_nb;
|
||||
|
||||
private:
|
||||
int btype, index_nb;
|
||||
int btype;
|
||||
double cutsq;
|
||||
|
||||
class NeighList *list;
|
||||
|
||||
@ -221,11 +221,6 @@ FixRHEOThermal::FixRHEOThermal(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
FixRHEOThermal::~FixRHEOThermal()
|
||||
{
|
||||
// Remove custom property if it exists
|
||||
int tmp1, tmp2, index;
|
||||
index = atom->find_custom("rheo_conductivity", tmp1, tmp2);
|
||||
if (index != -1) atom->remove_custom(index, 1, 0);
|
||||
|
||||
memory->destroy(cv_style);
|
||||
memory->destroy(Tc_style);
|
||||
memory->destroy(kappa_style);
|
||||
|
||||
@ -115,7 +115,7 @@ void PairRHEO::compute(int eflag, int vflag)
|
||||
|
||||
double **fp_store, *chi;
|
||||
if (compute_interface) {
|
||||
fp_store = compute_interface->fp_store;
|
||||
fp_store = atom->darray[compute_interface->index_fp_store];
|
||||
chi = compute_interface->chi;
|
||||
|
||||
for (i = 0; i < atom->nmax; i++) {
|
||||
@ -539,7 +539,7 @@ double PairRHEO::init_one(int i, int j)
|
||||
|
||||
int PairRHEO::pack_reverse_comm(int n, int first, double *buf)
|
||||
{
|
||||
double **fp_store = compute_interface->fp_store;
|
||||
double **fp_store = atom->darray[compute_interface->index_fp_store];
|
||||
int m = 0;
|
||||
int last = first + n;
|
||||
for (int i = first; i < last; i++) {
|
||||
@ -555,7 +555,7 @@ int PairRHEO::pack_reverse_comm(int n, int first, double *buf)
|
||||
|
||||
void PairRHEO::unpack_reverse_comm(int n, int *list, double *buf)
|
||||
{
|
||||
double **fp_store = compute_interface->fp_store;
|
||||
double **fp_store = atom->darray[compute_interface->index_fp_store];
|
||||
int m = 0;
|
||||
for (int i = 0; i < n; i++) {
|
||||
int j = list[i];
|
||||
|
||||
@ -125,7 +125,7 @@ void DisplaceAtoms::command(int narg, char **arg)
|
||||
else if (strcmp(arg[2],"y") == 0) d_dim = 1;
|
||||
else if (strcmp(arg[2],"z") == 0) d_dim = 2;
|
||||
else error->all(FLERR, 2, "Unknown displace_atoms ramp dimension {}", arg[2]);
|
||||
if ((domain->dimension == 2) && (d_dim = 2))
|
||||
if ((domain->dimension == 2) && (d_dim == 2))
|
||||
error->all(FLERR, 2, "Must not displace atoms in z-direction with 2d system");
|
||||
|
||||
double d_lo,d_hi;
|
||||
|
||||
@ -78,8 +78,8 @@ DumpCustom::DumpCustom(LAMMPS *lmp, int narg, char **arg) :
|
||||
|
||||
int ioffset = 5;
|
||||
expand = 0;
|
||||
nfield = nargnew = utils::expand_args(FLERR,narg-5,&arg[5],1,earg,lmp);
|
||||
if (earg != &arg[5]) expand = 1;
|
||||
nfield = nargnew = utils::expand_args(FLERR,narg-ioffset,&arg[ioffset],1,earg,lmp);
|
||||
if (earg != &arg[ioffset]) expand = 1;
|
||||
|
||||
// allocate field vectors
|
||||
|
||||
|
||||
10
src/info.cpp
10
src/info.cpp
@ -44,7 +44,9 @@
|
||||
#include "region.h"
|
||||
#include "update.h"
|
||||
#include "variable.h"
|
||||
#ifndef FMT_STATIC_THOUSANDS_SEPARATOR
|
||||
#include "fmt/chrono.h"
|
||||
#endif
|
||||
|
||||
#include <cctype>
|
||||
#include <cmath>
|
||||
@ -270,8 +272,16 @@ void Info::command(int narg, char **arg)
|
||||
if (out == nullptr) return;
|
||||
|
||||
fputs("\nInfo-Info-Info-Info-Info-Info-Info-Info-Info-Info-Info\n",out);
|
||||
#if defined(FMT_STATIC_THOUSANDS_SEPARATOR)
|
||||
{
|
||||
time_t tv = time(nullptr);
|
||||
struct tm *now = localtime(&tv);
|
||||
utils::print(out, "Printed on {}", asctime(now));
|
||||
}
|
||||
#else
|
||||
std::tm now = fmt::localtime(std::time(nullptr));
|
||||
utils::print(out,"Printed on {}", std::asctime(&now));
|
||||
#endif
|
||||
|
||||
if (flags & CONFIG) {
|
||||
utils::print(out,"\nLAMMPS version: {} / {}\n", lmp->version, lmp->num_ver);
|
||||
|
||||
@ -17,7 +17,7 @@
|
||||
// wrapper around including the JSON parsing and writing class
|
||||
// Do NOT include in any header file
|
||||
|
||||
#include "../third_party/nlohmann/json.hpp"
|
||||
#include "nlohmann/json.hpp"
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
using json = ::nlohmann_lmp::json;
|
||||
|
||||
@ -24,6 +24,7 @@
|
||||
#include "atom.h"
|
||||
#include "atom_vec.h"
|
||||
#include "comm.h"
|
||||
#include "command.h"
|
||||
#include "compute.h"
|
||||
#include "domain.h"
|
||||
#include "dump.h"
|
||||
@ -42,6 +43,7 @@
|
||||
#include "molecule.h"
|
||||
#include "neigh_list.h"
|
||||
#include "neighbor.h"
|
||||
#include "neigh_request.h"
|
||||
#include "output.h"
|
||||
#include "pair.h"
|
||||
#if defined(LMP_PLUGIN)
|
||||
@ -6119,6 +6121,96 @@ int lammps_find_compute_neighlist(void *handle, const char *id, int reqid) {
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
// helper Command class for a single neighbor list build
|
||||
|
||||
namespace LAMMPS_NS {
|
||||
class NeighProxy : protected Command
|
||||
{
|
||||
public:
|
||||
NeighProxy(class LAMMPS *lmp) : Command(lmp), neigh_idx(-1) {};
|
||||
|
||||
void command(int, char **) override;
|
||||
int get_index() const { return neigh_idx; }
|
||||
protected:
|
||||
int neigh_idx;
|
||||
};
|
||||
}
|
||||
|
||||
void NeighProxy::command(int narg, char **arg)
|
||||
{
|
||||
neigh_idx = -1;
|
||||
if (narg != 3) return;
|
||||
auto *req = neighbor->add_request(this, arg[0]);
|
||||
int flags = atoi(arg[1]);
|
||||
double cutoff = atof(arg[2]);
|
||||
req->apply_flags(flags);
|
||||
if (cutoff > 0.0) req->set_cutoff(cutoff);
|
||||
lmp->init();
|
||||
|
||||
// setup domain, communication and neighboring
|
||||
// acquire ghosts and build standard neighbor lists
|
||||
|
||||
if (domain->triclinic) domain->x2lamda(atom->nlocal);
|
||||
domain->pbc();
|
||||
domain->reset_box();
|
||||
comm->setup();
|
||||
if (neighbor->style) neighbor->setup_bins();
|
||||
comm->exchange();
|
||||
comm->borders();
|
||||
if (domain->triclinic) domain->lamda2x(atom->nlocal + atom->nghost);
|
||||
neighbor->build(1);
|
||||
|
||||
// build neighbor list this command needs based on earlier request
|
||||
|
||||
auto list = neighbor->find_list(this);
|
||||
neighbor->build_one(list);
|
||||
|
||||
// find neigh list
|
||||
for (int i = 0; i < neighbor->nlist; i++) {
|
||||
NeighList *list = neighbor->lists[i];
|
||||
if (this == list->requestor) {
|
||||
neigh_idx = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/** Build a single neighbor list in between runs and return its index
|
||||
*
|
||||
* A neighbor list request is created and the neighbor list built from a
|
||||
* proxy command. The request ID is typically 0, but will be
|
||||
* > 0 in case a compute has multiple neighbor list requests.
|
||||
*
|
||||
* \param handle pointer to a previously created LAMMPS instance cast to ``void *``.
|
||||
* \param id Identifier of neighbor list request
|
||||
* \param flags Neighbor list flags
|
||||
* \param cutoff Custom neighbor list cutoff if > 0, use default cutoff if < 0
|
||||
* \return return neighbor list index if valid, otherwise -1 */
|
||||
|
||||
int lammps_request_single_neighlist(void *handle, const char *id, int flags, double cutoff) {
|
||||
auto lmp = (LAMMPS *)handle;
|
||||
int idx = -1;
|
||||
if (!lmp || !lmp->error || !lmp->neighbor) {
|
||||
lammps_last_global_errormessage = fmt::format("ERROR: {}(): Invalid LAMMPS handle\n", FNERR);
|
||||
return -1;
|
||||
}
|
||||
|
||||
BEGIN_CAPTURE
|
||||
{
|
||||
NeighProxy proxy(lmp);
|
||||
std::vector<std::string> args = {id, std::to_string(flags), std::to_string(cutoff)};
|
||||
std::vector<char *> c_args;
|
||||
std::transform(args.begin(), args.end(), std::back_inserter(c_args),
|
||||
[](const std::string& s) { return (char *)s.c_str(); });
|
||||
proxy.command(static_cast<int>(c_args.size()), c_args.data());
|
||||
idx = proxy.get_index();
|
||||
}
|
||||
END_CAPTURE
|
||||
return idx;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
/** Return the number of entries in the neighbor list with given index
|
||||
*
|
||||
* \param handle pointer to a previously created LAMMPS instance cast to ``void *``.
|
||||
|
||||
@ -109,6 +109,17 @@ enum _LMP_VAR_CONST {
|
||||
LMP_VAR_STRING = 3 /*!< return value will be a string (catch-all) */
|
||||
};
|
||||
|
||||
/** Neighbor list settings constants
|
||||
*
|
||||
* Must be kept in sync with the equivalent constants in ``python/lammps/constants.py``,
|
||||
* ``fortran/lammps.f90``, ``tools/swig/lammps.i``, and
|
||||
* ``examples/COUPLE/plugin/liblammpsplugin.h`` */
|
||||
|
||||
enum _LMP_NEIGH_CONST {
|
||||
LMP_NEIGH_HALF = 0, /*!< request (default) half neighbor list */
|
||||
LMP_NEIGH_FULL = 1, /*!< request full neighbor list */
|
||||
};
|
||||
|
||||
/* Ifdefs to allow this file to be included in C and C++ programs */
|
||||
|
||||
#ifdef __cplusplus
|
||||
@ -235,6 +246,7 @@ int lammps_create_atoms(void *handle, int n, const int64_t *id, const int *type,
|
||||
int lammps_find_pair_neighlist(void *handle, const char *style, int exact, int nsub, int request);
|
||||
int lammps_find_fix_neighlist(void *handle, const char *id, int request);
|
||||
int lammps_find_compute_neighlist(void *handle, const char *id, int request);
|
||||
int lammps_request_single_neighlist(void *handle, const char *id, int flags, double cutoff);
|
||||
int lammps_neighlist_num_elements(void *handle, int idx);
|
||||
void lammps_neighlist_element_neighbors(void *handle, int idx, int element, int *iatom,
|
||||
int *numneigh, int **neighbors);
|
||||
|
||||
@ -113,7 +113,10 @@ void Min::init()
|
||||
// create fix needed for storing atom-based quantities
|
||||
// will delete it at end of run
|
||||
|
||||
fix_minimize = dynamic_cast<FixMinimize *>(modify->add_fix("MINIMIZE all MINIMIZE"));
|
||||
if (lmp->kokkos)
|
||||
fix_minimize = dynamic_cast<FixMinimize *>(modify->add_fix("MINIMIZE all MINIMIZE/kk"));
|
||||
else
|
||||
fix_minimize = dynamic_cast<FixMinimize *>(modify->add_fix("MINIMIZE all MINIMIZE"));
|
||||
|
||||
// clear out extra global and per-atom dof
|
||||
// will receive requests for new per-atom dof during pair init()
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -15,7 +15,9 @@
|
||||
|
||||
#include "comm.h"
|
||||
#include "error.h"
|
||||
#ifndef FMT_STATIC_THOUSANDS_SEPARATOR
|
||||
#include "fmt/chrono.h"
|
||||
#endif
|
||||
#include "tokenizer.h"
|
||||
|
||||
#include <cstring>
|
||||
@ -264,8 +266,15 @@ void Timer::modify_params(int narg, char **arg)
|
||||
// format timeout setting
|
||||
std::string timeout = "off";
|
||||
if (_timeout >= 0.0) {
|
||||
#if defined(FMT_STATIC_THOUSANDS_SEPARATOR)
|
||||
char outstr[200];
|
||||
struct tm *tv = gmtime(&((time_t) _timeout));
|
||||
strftime(outstr, 200, "%02d:%M:%S", tv);
|
||||
timeout = outstr;
|
||||
#else
|
||||
std::tm tv = fmt::gmtime((std::time_t) _timeout);
|
||||
timeout = fmt::format("{:02d}:{:%M:%S}", tv.tm_yday * 24 + tv.tm_hour, tv);
|
||||
#endif
|
||||
}
|
||||
|
||||
utils::logmesg(lmp, "New timer settings: style={} mode={} timeout={}\n", timer_style[_level],
|
||||
|
||||
@ -20,7 +20,9 @@
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "fix.h"
|
||||
#ifndef FMT_STATIC_THOUSANDS_SEPARATOR
|
||||
#include "fmt/chrono.h"
|
||||
#endif
|
||||
#include "input.h"
|
||||
#include "label_map.h"
|
||||
#include "memory.h"
|
||||
@ -1953,8 +1955,15 @@ int utils::date2num(const std::string &date)
|
||||
std::string utils::current_date()
|
||||
{
|
||||
time_t tv = time(nullptr);
|
||||
#if defined(FMT_STATIC_THOUSANDS_SEPARATOR)
|
||||
char outstr[200];
|
||||
struct tm *today = localtime(&tv);
|
||||
strftime(outstr, 200, "%Y-%m-%d", today);
|
||||
return std::string(outstr);
|
||||
#else
|
||||
std::tm today = fmt::localtime(tv);
|
||||
return fmt::format("{:%Y-%m-%d}", today);
|
||||
#endif
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
|
||||
14
third_party/README
vendored
14
third_party/README
vendored
@ -1 +1,13 @@
|
||||
This folder contains copies of third-party software
|
||||
This folder contains copies of third-party software.
|
||||
These packages may be either included in the source
|
||||
code tarball or downloaded on-the-fly during compilation
|
||||
as needed. In some cases, some small customizations
|
||||
are made to avoid namespace conflicts where the
|
||||
included third party software needs workarounds for
|
||||
issues on platforms (OS or compiler related) to
|
||||
function consistently with LAMMPS.
|
||||
|
||||
The following software source packages are included
|
||||
|
||||
Folder: Package: Version: Original URL
|
||||
nlohmann/json JSON for Modern C++ 3.12.0 https://github.com/nlohmann/json/releases/tag/v3.12.0
|
||||
|
||||
@ -27,6 +27,7 @@ Please see the individual tabulation scripts in this folder for examples:
|
||||
| wall_harmonic_tabulate.py | creates a table for fix wall/table with a simple repulsive harmonic potential |
|
||||
| wall_multi_tabulate.py | creates a table for fix wall/table with multiple tables |
|
||||
| pair_bi_tabulate.py | creates a table from radial distribution file using Boltzmann Inversion |
|
||||
| plot_forces.py | plots and extracts tabulated forces from table files |
|
||||
|
||||
Common command line flags:
|
||||
|
||||
|
||||
278
tools/tabulate/plot_forces.py
Executable file
278
tools/tabulate/plot_forces.py
Executable file
@ -0,0 +1,278 @@
|
||||
#!/usr/bin/env python3
|
||||
# -*- coding: utf-8 -*-
|
||||
# Author: Germain Clavier (Unicaen), germain.clavier at unicaen.fr
|
||||
|
||||
"""
|
||||
Plot LAMMPS tabulated forces.
|
||||
"""
|
||||
|
||||
import argparse
|
||||
import numpy as np
|
||||
import os
|
||||
import logging
|
||||
import sys
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
|
||||
units = {
|
||||
"lj": {
|
||||
"Distance": "Reduced units",
|
||||
"Energy": "Reduced units",
|
||||
"Force": "Reduced units",
|
||||
"kb": 1,
|
||||
},
|
||||
"real": {
|
||||
"Distance": "[A]",
|
||||
"Energy": "[kcal/mol]",
|
||||
"Force": "[kcal/mol/A]",
|
||||
"kb": 0.001985875,
|
||||
},
|
||||
"metal": {
|
||||
"Distance": "[A]",
|
||||
"Energy": "[eV]",
|
||||
"Force": "[eV/A]",
|
||||
"kb": 8.6173332e-5,
|
||||
},
|
||||
"si": {"Distance": "[m]", "Energy": "[J]", "Force": "[N]", "kb": 1.380649e-23},
|
||||
"cgs": {
|
||||
"Distance": "[cm]",
|
||||
"Energy": "[ergs]",
|
||||
"Force": "[dynes]",
|
||||
"kb": 1.3806504e-16,
|
||||
},
|
||||
"electron": {
|
||||
"Distance": "[Bohr]",
|
||||
"Energy": "[Hartrees]",
|
||||
"Force": "[Hartree/Bohr]",
|
||||
"kb": 3.16681534e-6,
|
||||
},
|
||||
"micro": {
|
||||
"Distance": "[um]",
|
||||
"Energy": "[pg·um^2/us^2]",
|
||||
"Force": "[pg·um/us^2]",
|
||||
"kb": 1.3806504e-8,
|
||||
},
|
||||
"nano": {
|
||||
"Distance": "[nm]",
|
||||
"Energy": "[ag·nm^2/ns^2]",
|
||||
"Force": "[ag·nm/ns^2]",
|
||||
"kb": 0.013806504,
|
||||
},
|
||||
}
|
||||
|
||||
|
||||
def compute_energy(tp):
|
||||
r = tp[0]
|
||||
fo = tp[2]
|
||||
e = np.zeros(r.shape)
|
||||
for i, (ri, fi) in enumerate(zip(r, fo)):
|
||||
if i == 0:
|
||||
continue
|
||||
dr = ri - r[i - 1]
|
||||
e[i] = e[i - 1] - dr * fo[i - 1]
|
||||
e -= e[-1]
|
||||
return e
|
||||
|
||||
|
||||
def main():
|
||||
|
||||
parser = argparse.ArgumentParser(
|
||||
description="""
|
||||
Plots LAMMPS tabulated forces. This script takes a table
|
||||
file as an input and plots all the tabulated forces inside with their
|
||||
corresponding energy. The forces label is the token used to name the
|
||||
force in the file. It can be used to output all the forces in separate
|
||||
files and/or recompute the energy from forces through finite difference
|
||||
(assuming e(rc)=0). This script requires the matplotlib and numpy
|
||||
Python libraries. Bitmap format is not supported.
|
||||
"""
|
||||
)
|
||||
parser.add_argument(
|
||||
"-u",
|
||||
"--units",
|
||||
dest="units",
|
||||
default="real",
|
||||
help="Units of the file (LAMMPS units system)",
|
||||
)
|
||||
parser.add_argument(
|
||||
"-f",
|
||||
"--file",
|
||||
dest="infile",
|
||||
default="",
|
||||
help="File to read",
|
||||
)
|
||||
parser.add_argument(
|
||||
"-x",
|
||||
dest="xrange",
|
||||
default="",
|
||||
help="xrange separated by : (for negative values use the '=' sign: -x=-3:10)",
|
||||
)
|
||||
parser.add_argument(
|
||||
"-y",
|
||||
dest="yrange",
|
||||
default="",
|
||||
help="yrange separated by :",
|
||||
)
|
||||
parser.add_argument(
|
||||
"-t",
|
||||
dest="temp",
|
||||
default=None,
|
||||
type=float,
|
||||
help="temperature for KbT plot [default none]",
|
||||
)
|
||||
parser.add_argument(
|
||||
"-d",
|
||||
"--diff-num",
|
||||
dest="recompute",
|
||||
action="store_true",
|
||||
help="Recompute the energies from forces and distances through finite differences",
|
||||
)
|
||||
parser.add_argument(
|
||||
"-e",
|
||||
dest="extract",
|
||||
action="store_true",
|
||||
help="Extract the forces in separate files",
|
||||
)
|
||||
args = parser.parse_args()
|
||||
logging.basicConfig(level=logging.INFO)
|
||||
|
||||
##########
|
||||
# Manage arguments
|
||||
|
||||
udistance = units[args.units]["Distance"]
|
||||
uenergy = units[args.units]["Energy"]
|
||||
uforce = units[args.units]["Force"]
|
||||
kb = units[args.units]["kb"]
|
||||
rlabel = " ".join(["Rij", udistance])
|
||||
elabel = " ".join(["E", uenergy])
|
||||
flabel = " ".join(["F", uforce])
|
||||
etitle = "Energy"
|
||||
ftitle = "Force"
|
||||
font = "DejaVu Sans"
|
||||
fontsize = 30
|
||||
|
||||
infile = args.infile
|
||||
if not os.path.isfile(infile):
|
||||
logger.error("Input file not found")
|
||||
sys.exit(1)
|
||||
|
||||
toplot = []
|
||||
with open(infile, "r") as f:
|
||||
lines = iter(f.readlines())
|
||||
while True:
|
||||
try:
|
||||
r = []
|
||||
force = []
|
||||
ener = []
|
||||
tok = []
|
||||
while not tok:
|
||||
tok = next(lines).partition("#")[0].rstrip()
|
||||
logger.info("Found {} token".format(tok))
|
||||
infos = next(lines).split()
|
||||
npoints = int(infos[1])
|
||||
next(lines)
|
||||
if "bitmap" in infos:
|
||||
logger.info("Unsupported bitmap format for token {:s}".format(tok))
|
||||
for _ in range(npoints):
|
||||
continue
|
||||
else:
|
||||
for i in range(npoints):
|
||||
line = next(lines).split()
|
||||
r.append(float(line[1]))
|
||||
ener.append(float(line[2]))
|
||||
force.append(float(line[3]))
|
||||
r = np.array(r)
|
||||
ener = np.array(ener)
|
||||
force = np.array(force)
|
||||
toplot.append([r, ener, force, tok])
|
||||
tok = []
|
||||
next(lines)
|
||||
except StopIteration:
|
||||
break
|
||||
if args.recompute:
|
||||
etitle = "Estimated energy"
|
||||
for tp in toplot:
|
||||
tp[1] = compute_energy(tp)
|
||||
|
||||
fig, axes = plt.subplots(1, 2)
|
||||
|
||||
for tp in toplot:
|
||||
axes[0].plot(tp[0], tp[1], label=tp[3], linewidth=3)
|
||||
axes[1].plot(tp[0], tp[2], label=tp[3], linewidth=3)
|
||||
hmin, hmax = axes[1].get_xlim()
|
||||
axes[1].hlines(0, hmin, hmax, color="black", linewidth=3, linestyles="dashdot")
|
||||
|
||||
if args.temp:
|
||||
if args.temp > 0:
|
||||
hmin, hmax = axes[0].get_xlim()
|
||||
axes[0].hlines(
|
||||
kb * args.temp,
|
||||
hmin,
|
||||
hmax,
|
||||
color="orange",
|
||||
label=r"$k_BT$",
|
||||
linewidth=3,
|
||||
linestyles="dashdot",
|
||||
)
|
||||
axes[0].text(hmax / 2.0, kb * args.temp, "KbT", fontsize=0.7 * fontsize)
|
||||
logger.info("KbT value= {:e} {:s}".format(kb * args.temp, uenergy))
|
||||
else:
|
||||
logger.info("Invalid temperature value: {:e}".format(args.temp))
|
||||
|
||||
if args.xrange:
|
||||
xmin, xmax = list(map(float, args.xrange.split(":")))
|
||||
axes[0].set_xlim(xmin, xmax)
|
||||
axes[1].set_xlim(xmin, xmax)
|
||||
if args.yrange:
|
||||
ymin, ymax = list(map(float, args.yrange.split(":")))
|
||||
axes[0].set_ylim(ymin, ymax)
|
||||
axes[1].set_ylim(ymin, ymax)
|
||||
|
||||
# Setting axes 0
|
||||
axes[0].set_title(etitle, fontsize=fontsize)
|
||||
axes[0].set_xlabel(
|
||||
rlabel, fontname=font, fontsize=fontsize
|
||||
) # xlabel name, size 30pts
|
||||
axes[0].set_ylabel(
|
||||
elabel, fontname=font, fontsize=fontsize
|
||||
) # ylabel name, size 30pts
|
||||
axes[0].tick_params(
|
||||
axis="both", which="major", labelsize=fontsize
|
||||
) # Biggers ticks, bigger tick labels!
|
||||
|
||||
# Setting axes 1
|
||||
axes[1].set_title(ftitle, fontsize=fontsize)
|
||||
axes[1].legend(frameon=False, fontsize=fontsize) # Fat font, no frame
|
||||
axes[1].set_xlabel(
|
||||
rlabel, fontname=font, fontsize=fontsize
|
||||
)
|
||||
axes[1].set_ylabel(
|
||||
flabel, fontname=font, fontsize=fontsize
|
||||
)
|
||||
axes[1].tick_params(
|
||||
axis="both", which="major", labelsize=fontsize
|
||||
)
|
||||
|
||||
figManager = plt.get_current_fig_manager()
|
||||
figManager.window.showMaximized()
|
||||
plt.show()
|
||||
|
||||
if args.extract:
|
||||
for tp in toplot:
|
||||
outfile = "".join([tp[3], ".plot"])
|
||||
logger.info("Writing file {}".format(outfile))
|
||||
with open(outfile, "w") as f:
|
||||
f.write("# {} force extracted from {}\n".format(tp[3], infile))
|
||||
f.write("# {:^20} {:^20} {:^20}\n".format("r", "energy", "force"))
|
||||
for a, b, c in zip(tp[0], tp[1], tp[2]):
|
||||
f.write("{:>18.16e} {:>18.16e} {:>18.16e}\n".format(a, b, c))
|
||||
return
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
try:
|
||||
main()
|
||||
except KeyboardInterrupt:
|
||||
raise SystemExit("User interruption.")
|
||||
106
unittest/force-styles/tests/mol-pair-lj_pirani.yaml
Normal file
106
unittest/force-styles/tests/mol-pair-lj_pirani.yaml
Normal file
@ -0,0 +1,106 @@
|
||||
---
|
||||
lammps_version: 2 Apr 2025
|
||||
tags: generated
|
||||
date_generated: Fri Apr 18 00:29:06 2025
|
||||
epsilon: 5e-14
|
||||
skip_tests:
|
||||
prerequisites: ! |
|
||||
atom full
|
||||
pair lj/pirani
|
||||
pre_commands: ! ""
|
||||
post_commands: ! |
|
||||
pair_modify shift yes
|
||||
input_file: in.fourmol
|
||||
pair_style: lj/pirani 8.0
|
||||
pair_coeff: ! |
|
||||
1 1 1.0 12.0 6.0 2.5 0.02
|
||||
1 2 1.5 9.0 6.0 1.75 0.01
|
||||
1 3 2.0 10.0 4.0 2.85 0.2
|
||||
1 4 0.0 12.0 6.0 2.8 1.173205
|
||||
1 5 0.0 12.0 6.0 2.8 1.173205
|
||||
2 2 3.0 12.0 6.0 1.0 0.005
|
||||
2 3 3.0 12.0 6.0 2.1 0.01
|
||||
2 4 4.0 12.0 6.0 0.5 0.005
|
||||
2 5 4.0 12.0 6.0 0.5 0.005
|
||||
3 3 0.0 12.0 6.0 3.2 0.02
|
||||
3 4 2.0 10.0 8.0 3.15 0.0173205
|
||||
3 5 0.0 12.0 6.0 3.15 0.0173205
|
||||
4 4 3.0 12.0 7.0 3.1 0.015
|
||||
4 5 1.0 11.0 6.0 3.1 0.015
|
||||
5 5 1.0 10.0 6.0 3.1 0.015
|
||||
extract: |
|
||||
alpha 2
|
||||
beta 2
|
||||
gamma 2
|
||||
rm 2
|
||||
epsilon 2
|
||||
natoms: 29
|
||||
init_vdwl: 5331.956709045601
|
||||
init_coul: 0
|
||||
init_stress: ! |2-
|
||||
3.7750714745234441e+03 3.5537150597026862e+03 5.6902162916903799e+04 -2.8072005838112923e+03 -9.7491364848488483e+03 1.2729560029569215e+04
|
||||
init_forces: ! |2
|
||||
1 -2.3166507408966854e+00 1.1814050990372530e+02 1.2954641659287341e+02
|
||||
2 5.8363749910673519e+01 4.8024381552510725e+01 -6.8694028688889247e+01
|
||||
3 -5.9503921864571154e+01 -1.6071817439843306e+02 -5.6597405190854595e+01
|
||||
4 -7.1134649049699927e-01 1.8613975279282308e-01 -5.0783713414591092e-01
|
||||
5 -2.4747264716871775e-01 -4.3787394242558070e-01 1.0743542961236709e+00
|
||||
6 -1.3941574974998655e+03 5.9741882401547746e+03 2.8500324441306744e+04
|
||||
7 1.1203760354219273e+03 -5.7518534790958884e+03 -2.8703199690479934e+04
|
||||
8 3.7622461947067976e+01 -2.8555055127577102e+01 1.4150262351684700e+02
|
||||
9 2.9060018894398823e+01 3.1299476554394232e+01 1.2859478553122568e+02
|
||||
10 1.9944480887666327e+02 -2.2491652641790279e+02 -8.6942871327890643e+01
|
||||
11 -2.0289808517921609e-01 -5.1391164971132541e-01 -8.4255333034028079e-01
|
||||
12 8.0482288086259217e+00 4.0940567086154154e+00 -4.6411440719924455e+00
|
||||
13 7.3840236567206285e-01 -2.9549777720168802e-01 -1.1519406459700064e-02
|
||||
14 -3.1304672652449900e-01 6.0809222534746646e-02 -7.8343451899926286e-01
|
||||
15 -2.9160840534461819e-02 7.6806966144169408e-01 2.7736307268099775e-01
|
||||
16 7.6601931092015529e+03 -5.3441050851312375e+03 -2.0311855298358816e+04
|
||||
17 -7.6559350002813917e+03 5.3350894691675485e+03 2.0332161066375302e+04
|
||||
18 -2.2850238529967973e-01 -6.8692273157008765e-01 5.2668630845652842e-01
|
||||
19 1.3993973537234641e-04 -7.4609137078560217e-05 6.5478121693019603e-04
|
||||
20 -3.8798038866996498e-04 -3.9611965614972961e-04 1.1277297053050416e-04
|
||||
21 -2.2335101557381928e-01 2.5953861388473182e-01 2.6125582506722322e-01
|
||||
22 3.1945020233767546e-04 2.3371631910053088e-04 2.7451639031626763e-04
|
||||
23 -4.6519356154706869e-04 -1.2233408621505282e-04 -5.6364589599674433e-05
|
||||
24 -5.5133701213992095e-02 -9.4107862424568284e-02 -1.3530576594828972e-01
|
||||
25 -8.3346433834854592e-04 -5.7684616458937175e-04 -1.8678312599956254e-03
|
||||
26 -4.1125298930647845e-04 -2.9178102339592801e-04 -2.9005455019886876e-04
|
||||
27 7.8624077687895158e-02 6.7068875684824725e-02 -5.6605077000657762e-02
|
||||
28 4.0630624646001851e-04 1.1596433235006506e-04 1.5226442697779881e-04
|
||||
29 -2.2503045892777234e-04 -1.4024118468637903e-05 -2.7955865290279905e-04
|
||||
run_vdwl: 2025.7316256875408
|
||||
run_coul: 0
|
||||
run_stress: ! |2-
|
||||
1.8959076601557997e+03 1.7586201688489862e+03 2.0716275757392985e+04 -1.3557176846491059e+03 -3.8359015598221822e+03 4.6552290404369523e+03
|
||||
run_forces: ! |2
|
||||
1 -1.8167639162581390e+00 1.1771418779705203e+02 1.2815251304772025e+02
|
||||
2 5.7615154471191410e+01 4.7476701906758862e+01 -6.7602559000311004e+01
|
||||
3 -6.0020160870473553e+01 -1.5903864235776544e+02 -5.6495701885917882e+01
|
||||
4 -7.0791746457219973e-01 1.8499145789736801e-01 -5.0659181307529255e-01
|
||||
5 -2.4638395816098327e-01 -4.3656328614875212e-01 1.0715255453838386e+00
|
||||
6 -6.6236838764715912e+02 1.9782901051305839e+03 8.2691846158655808e+03
|
||||
7 3.3632789638183766e+02 -1.7087851653075204e+03 -8.4826596480934113e+03
|
||||
8 9.2107904962862747e+01 -7.7141915019739500e+01 1.5412942131381334e+02
|
||||
9 2.8696841862882369e+01 3.0788008201635915e+01 1.2685653408680490e+02
|
||||
10 1.9840914201813521e+02 -2.2423605001860076e+02 -8.4220175413918682e+01
|
||||
11 -2.0136130409518435e-01 -5.0871695527363625e-01 -8.3420174587161200e-01
|
||||
12 7.9486965233487901e+00 4.1183091926914734e+00 -4.6935508545099953e+00
|
||||
13 7.3529063529646599e-01 -2.9290091034111942e-01 -1.1219670866151606e-02
|
||||
14 -3.0970372301220633e-01 5.8993643040875808e-02 -7.7531072022570446e-01
|
||||
15 -3.0754487700089653e-02 7.6978878292200181e-01 2.7960992063325796e-01
|
||||
16 2.9425585704800874e+03 -2.0582404076291327e+03 -7.7977754078264770e+03
|
||||
17 -2.9382680038667982e+03 2.0497188081055310e+03 7.8153023929895626e+03
|
||||
18 -2.2072907650908857e-01 -6.7402978043350725e-01 5.2241217519808170e-01
|
||||
19 1.3982600561567973e-04 -6.9464132281322043e-05 6.4718566866913177e-04
|
||||
20 -3.8223457087076577e-04 -3.9078559145284481e-04 1.1283430745541179e-04
|
||||
21 -2.2903028812118761e-01 2.6391221253513358e-01 2.7018653722712460e-01
|
||||
22 3.1793078542262739e-04 2.3386893987825232e-04 2.7595044210987341e-04
|
||||
23 -4.6582345719296879e-04 -1.2113347815506628e-04 -5.4088842128117200e-05
|
||||
24 -5.6206303605705644e-02 -9.5139911882907352e-02 -1.3714387491875332e-01
|
||||
25 -8.3067285016685810e-04 -5.7449793978904009e-04 -1.8648669696730458e-03
|
||||
26 -4.1342324759342764e-04 -2.9391975760244239e-04 -2.9182720701375784e-04
|
||||
27 7.7361153033234667e-02 6.6839436744107361e-02 -5.6399328691578474e-02
|
||||
28 4.0714285537424614e-04 1.1626525958259480e-04 1.5291178721883419e-04
|
||||
29 -2.2832773147555020e-04 -1.5023852976718537e-05 -2.7935291755901582e-04
|
||||
...
|
||||
@ -48,7 +48,11 @@ TEST(JSON, serialize_deserialize)
|
||||
j1["nothing"] = nullptr;
|
||||
|
||||
std::string expected = "{\"happy\":true,\"name\":\"Niels\",\"nothing\":null,\"pi\":3.141}";
|
||||
std::string dumped = j1.dump();
|
||||
std::string dumped = j1.dump(-1);
|
||||
ASSERT_THAT(expected, Eq(dumped));
|
||||
|
||||
expected = "{\n \"happy\": true,\n \"name\": \"Niels\",\n \"nothing\": null,\n \"pi\": 3.141\n}";
|
||||
dumped = j1.dump(2, ' ');
|
||||
ASSERT_THAT(expected, Eq(dumped));
|
||||
|
||||
json j2 = json::parse(expected);
|
||||
|
||||
Reference in New Issue
Block a user