Merge branch 'lammps:develop' into fortran-expansion

This commit is contained in:
hammondkd
2022-10-27 14:15:47 -05:00
committed by GitHub
26 changed files with 7346 additions and 40 deletions

View File

@ -7,6 +7,10 @@ cmake_minimum_required(VERSION 3.10)
if(POLICY CMP0074)
cmake_policy(SET CMP0074 NEW)
endif()
# set policy to silence warnings about ignoring ${CMAKE_REQUIRED_LIBRARIES} but use it
if(POLICY CMP0075)
cmake_policy(SET CMP0075 NEW)
endif()
# set policy to silence warnings about missing executable permissions in
# pythonx.y-config when cross-compiling. review occasionally if it may be set to NEW
if(POLICY CMP0109)
@ -385,9 +389,9 @@ pkg_depends(EXTRA-MOLECULE MOLECULE)
# detect if we may enable OpenMP support by default
set(BUILD_OMP_DEFAULT OFF)
find_package(OpenMP QUIET)
if(OpenMP_FOUND)
check_include_file_cxx(omp.h HAVE_OMP_H_INCLUDE)
find_package(OpenMP COMPONENTS CXX QUIET)
if(OpenMP_CXX_FOUND)
check_omp_h_include()
if(HAVE_OMP_H_INCLUDE)
set(BUILD_OMP_DEFAULT ON)
endif()
@ -396,8 +400,8 @@ endif()
option(BUILD_OMP "Build with OpenMP support" ${BUILD_OMP_DEFAULT})
if(BUILD_OMP)
find_package(OpenMP REQUIRED)
check_include_file_cxx(omp.h HAVE_OMP_H_INCLUDE)
find_package(OpenMP COMPONENTS CXX REQUIRED)
check_omp_h_include()
if(NOT HAVE_OMP_H_INCLUDE)
message(FATAL_ERROR "Cannot find the 'omp.h' header file required for full OpenMP support")
endif()

View File

@ -24,6 +24,21 @@ function(validate_option name values)
endif()
endfunction(validate_option)
# helper function to check for usable omp.h header
function(check_omp_h_include)
find_package(OpenMP COMPONENTS CXX QUIET)
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})
set(CMAKE_REQUIRED_LIBRARIES ${OpenMP_CXX_LIBRARIES})
check_include_file_cxx(omp.h _have_omp_h)
else()
set(_have_omp_h FALSE)
endif()
set(HAVE_OMP_H_INCLUDE ${_have_omp_h} PARENT_SCOPE)
endfunction()
# helper function for getting the most recently modified file or folder from a glob pattern
function(get_newest_file path variable)
file(GLOB _dirs ${path})

View File

@ -15,8 +15,9 @@ if(Kokkos_ENABLE_OPENMP)
if(NOT BUILD_OMP)
message(FATAL_ERROR "Must enable BUILD_OMP with Kokkos_ENABLE_OPENMP")
else()
if(LAMMPS_OMP_COMPAT_LEVEL LESS 4)
message(FATAL_ERROR "Compiler must support OpenMP 4.0 or later with Kokkos_ENABLE_OPENMP")
# NVHPC does not seem to provide a detectable OpenMP version, but is far beyond version 3.1
if((OpenMP_CXX_VERSION VERSION_LESS 3.1) AND NOT (CMAKE_CXX_COMPILER_ID STREQUAL "NVHPC"))
message(FATAL_ERROR "Compiler must support OpenMP 3.1 or later with Kokkos_ENABLE_OPENMP")
endif()
endif()
endif()

View File

@ -28,10 +28,3 @@ set(MPI_CXX "clang++" CACHE STRING "" FORCE)
set(MPI_CXX_COMPILER "mpicxx" CACHE STRING "" FORCE)
unset(HAVE_OMP_H_INCLUDE CACHE)
set(OpenMP_C "clang" CACHE STRING "" FORCE)
set(OpenMP_C_FLAGS "-fopenmp" CACHE STRING "" FORCE)
set(OpenMP_C_LIB_NAMES "omp" CACHE STRING "" FORCE)
set(OpenMP_CXX "clang++" CACHE STRING "" FORCE)
set(OpenMP_CXX_FLAGS "-fopenmp" CACHE STRING "" FORCE)
set(OpenMP_CXX_LIB_NAMES "omp" CACHE STRING "" FORCE)
set(OpenMP_omp_LIBRARY "libomp.so" CACHE PATH "" FORCE)

View File

@ -19,11 +19,3 @@ set(CMAKE_Fortran_FLAGS_DEBUG "-Wall -Og -g -std=f2003" CACHE STRING "" FORCE)
set(CMAKE_Fortran_FLAGS_RELWITHDEBINFO "-g -O2 -DNDEBUG -std=f2003" CACHE STRING "" FORCE)
set(CMAKE_Fortran_FLAGS_RELEASE "-O3 -DNDEBUG -std=f2003" CACHE STRING "" FORCE)
unset(HAVE_OMP_H_INCLUDE CACHE)
set(OpenMP_C "gcc" CACHE STRING "" FORCE)
set(OpenMP_C_FLAGS "-fopenmp" CACHE STRING "" FORCE)
set(OpenMP_C_LIB_NAMES "gomp" CACHE STRING "" FORCE)
set(OpenMP_CXX "g++" CACHE STRING "" FORCE)
set(OpenMP_CXX_FLAGS "-fopenmp" CACHE STRING "" FORCE)
set(OpenMP_CXX_LIB_NAMES "gomp" CACHE STRING "" FORCE)
set(OpenMP_omp_LIBRARY "libgomp.so" CACHE PATH "" FORCE)

View File

@ -0,0 +1,9 @@
# preset that will enable Nvidia HPC SDK compilers with support for MPI and OpenMP (on Linux boxes)
set(CMAKE_CXX_COMPILER "nvc++" CACHE STRING "" FORCE)
set(CMAKE_C_COMPILER "nvc" CACHE STRING "" FORCE)
set(CMAKE_Fortran_COMPILER "nvfortran" CACHE STRING "" FORCE)
set(MPI_CXX "nvc++" CACHE STRING "" FORCE)
set(MPI_CXX_COMPILER "mpicxx" CACHE STRING "" FORCE)
unset(HAVE_OMP_H_INCLUDE CACHE)

View File

@ -1,4 +1,4 @@
# preset that will restore gcc/g++ with support for MPI and OpenMP (on Linux boxes)
# preset that will set gcc/g++ with extra warnings enabled and support for MPI and OpenMP (on Linux boxes)
set(CMAKE_CXX_COMPILER "g++" CACHE STRING "" FORCE)
set(CMAKE_C_COMPILER "gcc" CACHE STRING "" FORCE)
@ -17,10 +17,3 @@ set(MPI_Fortran "gfortran" CACHE STRING "" FORCE)
set(MPI_Fortran_COMPILER "mpifort" CACHE STRING "" FORCE)
unset(HAVE_OMP_H_INCLUDE CACHE)
set(OpenMP_C "gcc" CACHE STRING "" FORCE)
set(OpenMP_C_FLAGS "-fopenmp" CACHE STRING "" FORCE)
set(OpenMP_C_LIB_NAMES "gomp" CACHE STRING "" FORCE)
set(OpenMP_CXX "g++" CACHE STRING "" FORCE)
set(OpenMP_CXX_FLAGS "-fopenmp" CACHE STRING "" FORCE)
set(OpenMP_CXX_LIB_NAMES "gomp" CACHE STRING "" FORCE)
set(OpenMP_omp_LIBRARY "libgomp.so" CACHE PATH "" FORCE)

View File

@ -7,10 +7,3 @@ set(MPI_CXX "pgc++" CACHE STRING "" FORCE)
set(MPI_CXX_COMPILER "mpicxx" CACHE STRING "" FORCE)
unset(HAVE_OMP_H_INCLUDE CACHE)
set(OpenMP_C "pgcc" CACHE STRING "" FORCE)
set(OpenMP_C_FLAGS "-mp" CACHE STRING "" FORCE)
set(OpenMP_C_LIB_NAMES "omp" CACHE STRING "" FORCE)
set(OpenMP_CXX "pgc++" CACHE STRING "" FORCE)
set(OpenMP_CXX_FLAGS "-mp" CACHE STRING "" FORCE)
set(OpenMP_CXX_LIB_NAMES "omp" CACHE STRING "" FORCE)
set(OpenMP_omp_LIBRARY "libomp.so" CACHE PATH "" FORCE)

View File

@ -295,6 +295,7 @@ OPT.
* :doc:`vashishta (gko) <pair_vashishta>`
* :doc:`vashishta/table (o) <pair_vashishta>`
* :doc:`wf/cut <pair_wf_cut>`
* :doc:`ylz <pair_ylz>`
* :doc:`yukawa (gko) <pair_yukawa>`
* :doc:`yukawa/colloid (go) <pair_yukawa_colloid>`
* :doc:`zbl (gko) <pair_zbl>`

View File

@ -200,6 +200,7 @@ particle models including ellipsoids, 2d lines, and 3d triangles.
* :doc:`Howto spherical <Howto_spherical>`
* :doc:`pair_style gayberne <pair_gayberne>`
* :doc:`pair_style resquared <pair_resquared>`
* :doc:`pair_style ylz <pair_ylz>`
* `doc/PDF/pair_gayberne_extra.pdf <PDF/pair_gayberne_extra.pdf>`_
* `doc/PDF/pair_resquared_extra.pdf <PDF/pair_resquared_extra.pdf>`_
* examples/ASPHERE

View File

@ -374,6 +374,7 @@ accelerated styles exist.
* :doc:`vashishta <pair_vashishta>` - Vashishta 2-body and 3-body potential
* :doc:`vashishta/table <pair_vashishta>` -
* :doc:`wf/cut <pair_wf_cut>` - Wang-Frenkel Potential for short-ranged interactions
* :doc:`ylz <pair_ylz>` - Yuan-Li-Zhang Potential for anisotropic interactions
* :doc:`yukawa <pair_yukawa>` - Yukawa potential
* :doc:`yukawa/colloid <pair_yukawa_colloid>` - screened Yukawa potential for finite-size particles
* :doc:`zbl <pair_zbl>` - Ziegler-Biersack-Littmark potential

177
doc/src/pair_ylz.rst Normal file
View File

@ -0,0 +1,177 @@
.. index:: pair_style ylz
pair_style ylz command
===========================
Syntax
""""""
.. code-block:: LAMMPS
pair_style ylz cutoff
* cutoff = global cutoff for interactions (distance units)
Examples
""""""""
.. code-block:: LAMMPS
pair_style ylz 2.6
pair_coeff * * 1.0 1.0 4 3 0.0 2.6
Description
"""""""""""
.. versionadded:: TBD
The *ylz* (Yuan-Li-Zhang) style computes an anisotropic interaction
between pairs of coarse-grained particles considering the relative
particle orientations. This potential was originally developed as a
particle-based solvent-free model for biological membranes
:ref:`(Yuan2010a) <Yuan>`. Unlike :doc:`pair_style gayberne
<pair_gayberne>`, whose orientation dependence is strictly derived from
the closest distance between two ellipsoidal rigid bodies, the
orientation-dependence of this pair style is mathematically defined such
that the particles can self-assemble into one-particle-thick fluid
membranes. The potential of this pair style is described by:
.. math::
U ( \mathbf{r}_{ij}, \mathbf{n}_i, \mathbf{n}_j ) =\left\{\begin{matrix} {u}_R(r)+\left [ 1-\phi (\mathbf{\hat{r}}_{ij}, \mathbf{n}_i, \mathbf{n}_j ) \right ]\epsilon, ~~ r<{r}_{min} \\ {u}_A(r)\phi (\mathbf{\hat{r}}_{ij}, \mathbf{n}_i, \mathbf{n}_j ),~~ {r}_{min}<r<{r}_{c} \\ \end{matrix}\right.\\\\ \phi (\mathbf{\hat{r}}_{ij}, \mathbf{n}_i, \mathbf{n}_j )=1+\left [ \mu (a(\mathbf{\hat{r}}_{ij}, \mathbf{n}_i, \mathbf{n}_j )-1) \right ] \\\\a(\mathbf{\hat{r}}_{ij}, \mathbf{n}_i, \mathbf{n}_j )=(\mathbf{n}_i\times\mathbf{\hat{r}}_{ij} )\cdot (\mathbf{n}_j\times\mathbf{\hat{r}}_{ij} )+{\beta}(\mathbf{n}_i-\mathbf{n}_j)\cdot \mathbf{\hat{r}}_{ij}-\beta^{2}\\\\ {u}_R(r)=\epsilon \left [ \left ( \frac{{r}_{min}}{r} \right )^{4}-2\left ( \frac{{r}_{min}}{r}\right )^{2} \right ] \\\\ {u}_A(r)=-\epsilon\;cos^{2\zeta }\left [ \frac{\pi}{2}\frac{\left ( {r}-{r}_{min} \right )}{\left ( {r}_{c}-{r}_{min} \right )} \right ]\\
where :math:`\mathbf{r}_{i}` and :math:`\mathbf{r}_{j}` are the center
position vectors of particles i and j, respectively,
:math:`\mathbf{r}_{ij}=\mathbf{r}_{i}-\mathbf{r}_{j}` is the
inter-particle distance vector, :math:`r=\left|\mathbf{r}_{ij} \right|`
and :math:`{\hat{\mathbf{r}}}_{ij}=\mathbf{r}_{ij}/r`. The unit vectors
:math:`\mathbf{n}_{i}` and :math:`\mathbf{n}_{j}` represent the axes of
symmetry of particles i and j, respectively, :math:`u_R` and :math:`u_A`
are the repulsive and attractive potentials, :math:`\phi` is an angular
function which depends on the relative orientation between pair
particles, :math:`\mu` is the parameter related to the bending rigidity
of the membrane, :math:`\beta` is the parameter related to the
spontaneous curvature, and :math:`\epsilon` is the energy unit,
respectively. The :math:`\zeta` controls the slope of the attractive
branch and hence the diffusivity of the particles in the in-plane
direction of the membrane. :math:`{r}_{c}` is the cutoff radius,
:math:`r_{min}` is the distance which minimizes the potential energy
:math:`u_{A}(r)` and :math:`r_{min}=2^{1/6}\sigma`, where :math:`\sigma`
is the length unit.
This pair style is suited for solvent-free coarse-grained simulations of
biological systems involving lipid bilayer membranes, such as vesicle
shape transformations :ref:`(Yuan2010b) <Yuan>`, nanoparticle
endocytosis :ref:`(Huang) <Huang>`, modeling of red blood cell membranes
:ref:`(Fu) <Fu>`, :ref:`(Appshaw) <Appshaw>`, and modeling of cell
elasticity :ref:`(Becton) <Becton>`.
Use of this pair style requires the NVE, NVT, or NPT fixes with the
*asphere* extension (e.g. :doc:`fix nve/asphere <fix_nve_asphere>`) in
order to integrate particle rotation. Additionally, :doc:`atom_style
ellipsoid <atom_style>` should be used since it defines the rotational
state of each particle.
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, or by
mixing as described below:
* :math:`\epsilon` = well depth (energy units)
* :math:`\sigma` = minimum effective particle radii (distance units)
* :math:`\zeta` = tune parameter for the slope of the attractive branch
* :math:`\mu` = parameter related to bending rigidity
* :math:`\beta` = parameter related to the spontaneous curvature
* cutoff (distance units)
The last coefficient is optional. If not specified, the global
cutoff specified in the pair_style command is used.
----------
Mixing, shift, table, tail correction, restart, rRESPA info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
For atom type pairs I,J and I != J, the epsilon and sigma coefficients
and cutoff distance for this pair style can be mixed. The default mix
value is *geometric*\ . See the "pair_modify" command for details.
The :doc:`pair_modify <pair_modify>` table option is 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 can only be used via the *pair* keyword of the
:doc:`run_style respa <run_style>` command. It does not support the
*inner*, *middle*, *outer* keywords.
----------
Restrictions
""""""""""""
The *ylz* style is part of the ASPHERE package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` page for more info.
This pair style requires that atoms store torque and a quaternion to
represent their orientation, as defined by the :doc:`atom_style
<atom_style>`. It also requires they store a per-atom :doc:`shape
<set>`. The particles cannot store a per-particle diameter. To avoid
being mistakenly considered as point particles, the shape parameters ought
to be non-spherical, like [1 0.99 0.99]. Unlike the :doc:`resquared
<pair_resquared>` pair style for which the shape directly determines the
mathematical expressions of the potential, the shape parameters for this
pair style is only involved in the computation of the moment of inertia
and thus only influences the rotational dynamics of individual
particles.
This pair style requires that **all** atoms are ellipsoids as defined by
the :doc:`atom_style ellipsoid <atom_style>` command.
Related commands
""""""""""""""""
:doc:`pair_coeff <pair_coeff>`, :doc:`fix nve/asphere
:doc:<fix_nve_asphere>`, `compute temp/asphere <compute_temp_asphere>`,
:doc::doc:`pair_style resquared <pair_resquared>`, :doc:`pair_style
:doc:gayberne <pair_gayberne>`
Default
"""""""
none
----------
.. _Yuan:
**(Yuan2010a)** Yuan, Huang, Li, Lykotrafitis, Zhang, Phys. Rev. E, 82, 011905(2010).
**(Yuan2010b)** Yuan, Huang, Zhang, Soft. Matter, 6, 4571(2010).
.. _Huang:
**(Huang)** Huang, Zhang, Yuan, Gao, Zhang, Nano Lett. 13, 4546(2013).
.. _Fu:
**(Fu)** Fu, Peng, Yuan, Kfoury, Young, Comput. Phys. Commun, 210, 193-203(2017).
.. _Appshaw:
**(Appshaw)** Appshaw, Seddon, Hanna, Soft. Matter,18, 1747(2022).
.. _Becton:
**(Becton)** Becton, Averett, Wang, Biomech. Model. Mechanobiology, 18, 425-433(2019).

View File

@ -123,6 +123,7 @@ Antonelli
api
Apoorva
Appl
Appshaw
apptainer
Apu
arallel
@ -181,6 +182,7 @@ Avalos
avalue
aveforce
Avendano
Averett
avi
AVX
awpmd
@ -237,6 +239,7 @@ bcolor
bdiam
bdw
Beckman
Becton
Belak
Bellott
bem
@ -270,6 +273,7 @@ bilayers
binsize
binstyle
binutils
Biomech
biomolecular
biomolecule
Biomolecules
@ -924,6 +928,7 @@ Emmrich
emol
eN
endian
endocytosis
energetics
energyCorr
eng
@ -1155,6 +1160,7 @@ ftm
ftol
fuer
fugacity
Fu
Fumi
func
funcs
@ -1306,6 +1312,7 @@ Hamaker
Hamel
Hammerschmidt
Hanley
Hanna
haptic
Haque
Hara
@ -1657,6 +1664,7 @@ kepler
keV
Keyes
keyfile
Kfoury
Khersonskii
Khrapak
Khvostov
@ -1912,6 +1920,7 @@ lwsock
lx
ly
Lybrand
Lykotrafitis
lyon
Lysogorskiy
Lyulin
@ -2022,6 +2031,7 @@ MEAM
meamf
meanDist
mech
Mechanobiology
Mecke
mediumaquamarine
mediumblue
@ -3143,6 +3153,7 @@ seagreen
Secor
sectoring
sed
Seddon
segmental
Seifert
Seleson
@ -3899,6 +3910,7 @@ yhi
yi
ylat
ylo
ylz
ymax
ymin
yml

View File

@ -29,8 +29,10 @@ stochastic rotation dynamics solvent, using the fix srd command.
box = 2d rigid box particles in SRDs, self-diffusion and viscosity
dimer = 2d rigid dimers in SRDs, self-diffusion and viscosity
ellipsoid = 2d ellipsoids in SRDs, self-diffusion and viscosity
flat_membrane = 2d flat membrane in NVE, self-diffusion and viscosity
line = 2d line-faceted rigid particles, NEMD shearing for viscosity,
implicit and in SRDs
poly = 2d poly-disperse spheroids in SRDs, self-diffusion and viscosity
star = 2d rigid star particles in SRDs, self-diffusion and viscosity
tri = 3d triangle-faceted rigid particles in SRDs, NEMD shearing for viscosity
vesicle = 3d vesicle in NVE, self-diffusionand viscosity

View File

@ -0,0 +1,55 @@
# flat membrane demo
variable r0 equal 0.97
variable d1 equal ${r0}
variable d2 equal sqrt(3.0)*${r0}
variable d3 equal 3.0*${r0}
variable ro equal 2./${d1}/${d2}/${d3}
variable T equal 0.23
variable LD equal 1.0
units lj
atom_style ellipsoid
boundary p p p
lattice custom ${ro} a1 ${d1} 0.0 0.0 a2 0.0 ${d2} 0.0 &
a3 0.0 0.0 ${d3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
region box block 0 40 0 24 -20 20
create_box 1 box
region membrane block 0 40 0 24 -0.5 0.5
create_atoms 1 region membrane
group membrane region membrane
set type 1 mass 1.0
set type 1 shape 1 0.99 0.99
set group all quat 0 -1 0 90
#compute memb all temp/com
#compute rot all temp/asphere bias memb
velocity all create ${T} 87287 loop geom
pair_style ylz 2.6
pair_coeff * * 1.0 1.0 4 3 0.0 2.6
neighbor 1.0 bin
thermo_style custom step temp press pxx pyy
thermo 200
timestep 0.01
#dump 1 all atom 10 dump_onlymembrane.lammpstrj
fix 1 all langevin ${T} ${T} ${LD} 48279
fix 2 all nve/asphere
run 3000

View File

@ -0,0 +1,159 @@
LAMMPS (15 Sep 2022)
using 1 OpenMP thread(s) per MPI task
# flat membrane demo
variable r0 equal 0.97
variable d1 equal ${r0}
variable d1 equal 0.97
variable d2 equal sqrt(3.0)*${r0}
variable d2 equal sqrt(3.0)*0.97
variable d3 equal 3.0*${r0}
variable d3 equal 3.0*0.97
variable ro equal 2./${d1}/${d2}/${d3}
variable ro equal 2./0.97/${d2}/${d3}
variable ro equal 2./0.97/1.68008928334181/${d3}
variable ro equal 2./0.97/1.68008928334181/2.91
variable T equal 0.23
variable LD equal 1.0
units lj
atom_style ellipsoid
boundary p p p
lattice custom ${ro} a1 ${d1} 0.0 0.0 a2 0.0 ${d2} 0.0 a3 0.0 0.0 ${d3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
lattice custom 0.421728460751825 a1 ${d1} 0.0 0.0 a2 0.0 ${d2} 0.0 a3 0.0 0.0 ${d3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
lattice custom 0.421728460751825 a1 0.97 0.0 0.0 a2 0.0 ${d2} 0.0 a3 0.0 0.0 ${d3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
lattice custom 0.421728460751825 a1 0.97 0.0 0.0 a2 0.0 1.68008928334181 0.0 a3 0.0 0.0 ${d3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
lattice custom 0.421728460751825 a1 0.97 0.0 0.0 a2 0.0 1.68008928334181 0.0 a3 0.0 0.0 2.91 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
Lattice spacing in x,y,z = 0.97 1.6800893 2.91
region box block 0 40 0 24 -20 20
create_box 1 box
Created orthogonal box = (0 0 -58.2) to (38.8 40.322143 58.2)
1 by 1 by 1 MPI processor grid
region membrane block 0 40 0 24 -0.5 0.5
create_atoms 1 region membrane
Created 1920 atoms
using lattice units in orthogonal box = (0 0 -58.2) to (38.8 40.322143 58.2)
create_atoms CPU = 0.001 seconds
group membrane region membrane
1920 atoms in group membrane
set type 1 mass 1.0
Setting atom values ...
1920 settings made for mass
set type 1 shape 1 1 1
Setting atom values ...
1920 settings made for shape
set group all quat 0 -1 0 90
Setting atom values ...
1920 settings made for quat
#compute memb all temp/com
#compute rot all temp/asphere bias memb
velocity all create ${T} 87287 loop geom
velocity all create 0.23 87287 loop geom
pair_style ylz 2.6
pair_coeff * * 1.0 1.0 4 3 0.0 2.6
neighbor 1.0 bin
thermo_style custom step temp press pxx pyy
thermo 200
timestep 0.01
#dump 1 all atom 10 dump_onlymembrane.lammpstrj
fix 1 all langevin ${T} ${T} ${LD} 48279
fix 1 all langevin 0.23 ${T} ${LD} 48279
fix 1 all langevin 0.23 0.23 ${LD} 48279
fix 1 all langevin 0.23 0.23 1 48279
fix 2 all nve/asphere
run 3000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- pair ylz command:
@Article{Yuan10,
author = {H. Yuan, C. Huang, J. Li, G. Lykotrafitis, and S. Zhang},
title = {One-particle-thick, solvent-free, coarse-grained model for biological and biomimetic fluid membranes},
journal = {Phys. Rev. E},
year = 2010,
volume = 82,
pages = {011905}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 10 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.6
ghost atom cutoff = 3.6
binsize = 1.8, bins = 22 23 65
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair ylz, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.024 | 5.024 | 5.024 Mbytes
Step Temp Press Pxx Pyy
0 0.23 -0.0073508785 -0.012283389 -0.012234574
200 0.20903886 -0.0010605951 -0.0011885957 -0.00198842
400 0.21898026 -0.00069250685 -0.0013217981 -0.00073225707
600 0.22689361 -0.00057919328 -0.00076880503 -0.0010242283
800 0.22983221 -0.00032145682 -0.00051928834 -0.00059337525
1000 0.23819392 -0.00027969126 -0.00088082301 -5.2666567e-05
1200 0.22053795 -0.00029571329 -0.0004446455 -0.00035529929
1400 0.22285021 -0.0002690371 -0.00068896571 -3.6258442e-05
1600 0.22687044 2.8599875e-05 -0.00032651798 0.0004056081
1800 0.23356905 -2.28742e-05 -0.00027073251 0.00025081131
2000 0.22499821 8.8230586e-06 -7.5750159e-05 0.0001988705
2200 0.23162995 -9.026855e-05 -0.00025832535 5.4904927e-05
2400 0.22920223 0.00016700455 3.5283125e-05 0.00034955857
2600 0.2260299 5.3095557e-05 0.00025691786 0.00013353467
2800 0.2296401 0.00043234854 0.00058344966 0.00063645193
3000 0.22564577 2.6423111e-05 8.9918406e-05 0.00022146229
Loop time of 6.76659 on 1 procs for 3000 steps with 1920 atoms
Performance: 383058.431 tau/day, 443.355 timesteps/s
99.8% CPU use with 1 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 5.7968 | 5.7968 | 5.7968 | 0.0 | 85.67
Neigh | 0.086077 | 0.086077 | 0.086077 | 0.0 | 1.27
Comm | 0.034761 | 0.034761 | 0.034761 | 0.0 | 0.51
Output | 0.00038014 | 0.00038014 | 0.00038014 | 0.0 | 0.01
Modify | 0.8268 | 0.8268 | 0.8268 | 0.0 | 12.22
Other | | 0.02181 | | | 0.32
Nlocal: 1920 ave 1920 max 1920 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Nghost: 772 ave 772 max 772 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Neighs: 46804 ave 46804 max 46804 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 46804
Ave neighs/atom = 24.377083
Neighbor list builds = 99
Dangerous builds = 0
Total wall time: 0:00:06

View File

@ -0,0 +1,159 @@
LAMMPS (15 Sep 2022)
using 1 OpenMP thread(s) per MPI task
# flat membrane demo
variable r0 equal 0.97
variable d1 equal ${r0}
variable d1 equal 0.97
variable d2 equal sqrt(3.0)*${r0}
variable d2 equal sqrt(3.0)*0.97
variable d3 equal 3.0*${r0}
variable d3 equal 3.0*0.97
variable ro equal 2./${d1}/${d2}/${d3}
variable ro equal 2./0.97/${d2}/${d3}
variable ro equal 2./0.97/1.68008928334181/${d3}
variable ro equal 2./0.97/1.68008928334181/2.91
variable T equal 0.23
variable LD equal 1.0
units lj
atom_style ellipsoid
boundary p p p
lattice custom ${ro} a1 ${d1} 0.0 0.0 a2 0.0 ${d2} 0.0 a3 0.0 0.0 ${d3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
lattice custom 0.421728460751825 a1 ${d1} 0.0 0.0 a2 0.0 ${d2} 0.0 a3 0.0 0.0 ${d3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
lattice custom 0.421728460751825 a1 0.97 0.0 0.0 a2 0.0 ${d2} 0.0 a3 0.0 0.0 ${d3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
lattice custom 0.421728460751825 a1 0.97 0.0 0.0 a2 0.0 1.68008928334181 0.0 a3 0.0 0.0 ${d3} basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
lattice custom 0.421728460751825 a1 0.97 0.0 0.0 a2 0.0 1.68008928334181 0.0 a3 0.0 0.0 2.91 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0
Lattice spacing in x,y,z = 0.97 1.6800893 2.91
region box block 0 40 0 24 -20 20
create_box 1 box
Created orthogonal box = (0 0 -58.2) to (38.8 40.322143 58.2)
1 by 1 by 4 MPI processor grid
region membrane block 0 40 0 24 -0.5 0.5
create_atoms 1 region membrane
Created 1920 atoms
using lattice units in orthogonal box = (0 0 -58.2) to (38.8 40.322143 58.2)
create_atoms CPU = 0.001 seconds
group membrane region membrane
1920 atoms in group membrane
set type 1 mass 1.0
Setting atom values ...
1920 settings made for mass
set type 1 shape 1 1 1
Setting atom values ...
1920 settings made for shape
set group all quat 0 -1 0 90
Setting atom values ...
1920 settings made for quat
#compute memb all temp/com
#compute rot all temp/asphere bias memb
velocity all create ${T} 87287 loop geom
velocity all create 0.23 87287 loop geom
pair_style ylz 2.6
pair_coeff * * 1.0 1.0 4 3 0.0 2.6
neighbor 1.0 bin
thermo_style custom step temp press pxx pyy
thermo 200
timestep 0.01
#dump 1 all atom 10 dump_onlymembrane.lammpstrj
fix 1 all langevin ${T} ${T} ${LD} 48279
fix 1 all langevin 0.23 ${T} ${LD} 48279
fix 1 all langevin 0.23 0.23 ${LD} 48279
fix 1 all langevin 0.23 0.23 1 48279
fix 2 all nve/asphere
run 3000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- pair ylz command:
@Article{Yuan10,
author = {H. Yuan, C. Huang, J. Li, G. Lykotrafitis, and S. Zhang},
title = {One-particle-thick, solvent-free, coarse-grained model for biological and biomimetic fluid membranes},
journal = {Phys. Rev. E},
year = 2010,
volume = 82,
pages = {011905}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 10 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.6
ghost atom cutoff = 3.6
binsize = 1.8, bins = 22 23 65
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair ylz, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.182 | 4.794 | 5.472 Mbytes
Step Temp Press Pxx Pyy
0 0.23 -0.0073508785 -0.012283389 -0.012234574
200 0.20647718 -0.0012368074 -0.0021167303 -0.0015343502
400 0.21648371 -0.00085695085 -0.0015627331 -0.0011177093
600 0.22929515 -0.00050218657 -0.0008332 -0.00062622609
800 0.22062664 -0.00049172378 -0.000611884 -0.00075089294
1000 0.22422425 -0.00039405068 -0.00037600355 -0.00070786572
1200 0.2298767 -0.00025939082 -0.00021616578 -0.00053125505
1400 0.2335927 5.8028332e-05 0.00017530192 -3.1675138e-05
1600 0.22884878 -0.0001733902 -0.0008056431 0.00014276754
1800 0.22813498 0.00019873459 0.00051040124 5.8860949e-05
2000 0.2273166 -3.3595127e-05 0.0001705632 -0.00026498213
2200 0.2251643 -2.4517311e-05 -4.0618888e-05 1.066658e-05
2400 0.22460629 -4.5661259e-05 -0.00019144039 -1.6649099e-05
2600 0.23085675 0.00014029405 0.00017983536 0.00017895001
2800 0.22364591 4.2999164e-05 -0.00011000466 0.00024363243
3000 0.23421357 0.00023505702 0.00020752013 0.00053567111
Loop time of 4.68577 on 4 procs for 3000 steps with 1920 atoms
Performance: 553164.568 tau/day, 640.237 timesteps/s
95.6% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.00015072 | 1.6029 | 3.8573 | 131.7 | 34.21
Neigh | 0.00055747 | 0.025423 | 0.065858 | 17.0 | 0.54
Comm | 0.0052259 | 0.48173 | 1.624 | 96.5 | 10.28
Output | 0.0003894 | 0.023428 | 0.047223 | 15.0 | 0.50
Modify | 0.00037337 | 0.2141 | 0.44595 | 46.3 | 4.57
Other | | 2.338 | | | 49.90
Nlocal: 480 ave 1011 max 0 min
Histogram: 2 0 0 0 0 0 0 0 1 1
Nghost: 860 ave 1771 max 0 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs: 11697.8 ave 30095 max 0 min
Histogram: 2 0 0 0 0 1 0 0 0 1
Total # of neighbors = 46791
Ave neighs/atom = 24.370313
Neighbor list builds = 99
Dangerous builds = 0
Total wall time: 0:00:04

View File

@ -0,0 +1,33 @@
# Vesicle demo
variable T equal 0.2
variable LD equal 1.0
units lj
atom_style ellipsoid
boundary p p p
read_data read_data.vesicle1026
compute ali all temp/com
compute rott all temp/asphere bias ali
velocity all create ${T} 87287 loop geom
pair_style ylz 2.6
pair_coeff * * 1.0 1.0 4 3 0.0 2.6
neighbor 1.0 bin
thermo_style custom step temp press pxx pyy
thermo 200
timestep 0.001
#dump 1 all atom 10 onlymembrane2.lammpstrj
fix 1 all langevin ${T} ${T} ${LD} 48279
fix 2 all nve/asphere
run 3000

View File

@ -0,0 +1,122 @@
LAMMPS (15 Sep 2022)
using 1 OpenMP thread(s) per MPI task
# Vesicle demo
variable T equal 0.2
variable LD equal 1.0
units lj
atom_style ellipsoid
boundary p p p
read_data read_data.vesicle1026
Reading data file ...
orthogonal box = (-35 -35 -35) to (35 35 35)
1 by 1 by 1 MPI processor grid
reading atoms ...
2938 atoms
2938 ellipsoids
read_data CPU = 0.010 seconds
compute ali all temp/com
compute rott all temp/asphere bias ali
velocity all create ${T} 87287 loop geom
velocity all create 0.2 87287 loop geom
pair_style ylz 2.6
pair_coeff * * 1.0 1.0 4 3 0.0 2.6
neighbor 1.0 bin
thermo_style custom step temp press pxx pyy
thermo 200
timestep 0.001
#dump 1 all atom 10 onlymembrane2.lammpstrj
fix 1 all langevin ${T} ${T} ${LD} 48279
fix 1 all langevin 0.2 ${T} ${LD} 48279
fix 1 all langevin 0.2 0.2 ${LD} 48279
fix 1 all langevin 0.2 0.2 1 48279
fix 2 all nve/asphere
run 3000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- pair ylz command:
@Article{Yuan10,
author = {H. Yuan, C. Huang, J. Li, G. Lykotrafitis, and S. Zhang},
title = {One-particle-thick, solvent-free, coarse-grained model for biological and biomimetic fluid membranes},
journal = {Phys. Rev. E},
year = 2010,
volume = 82,
pages = {011905}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 10 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.6
ghost atom cutoff = 3.6
binsize = 1.8, bins = 39 39 39
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair ylz, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 5.091 | 5.091 | 5.091 Mbytes
Step Temp Press Pxx Pyy
0 0.2 -0.0054891414 -0.0052713616 -0.0053641136
200 0.12816247 -0.0051288861 -0.0048542514 -0.0049847561
400 0.1377632 -0.0048071582 -0.0045651263 -0.0048444087
600 0.14983781 -0.0045305725 -0.0043305994 -0.0046127777
800 0.16205271 -0.0041176346 -0.0040534483 -0.0041351779
1000 0.17462122 -0.0037000069 -0.0034938539 -0.0037915494
1200 0.18335448 -0.0032674318 -0.0032790248 -0.0031967931
1400 0.19195613 -0.0029332101 -0.0030823703 -0.0028287799
1600 0.19261762 -0.0025263447 -0.0025152249 -0.0026205398
1800 0.19758674 -0.0021087725 -0.001981333 -0.002309048
2000 0.19748896 -0.0017662369 -0.0019316344 -0.0016696035
2200 0.20196986 -0.0013363214 -0.0015581191 -0.0013384961
2400 0.20109248 -0.0009190831 -0.0010331417 -0.0010664316
2600 0.20228664 -0.00053590675 -0.00071808747 -0.00050218533
2800 0.20512955 -0.00030845899 -0.00016244901 -0.00047877516
3000 0.19855941 -7.9520073e-05 -0.00014969215 -5.4724226e-06
Loop time of 9.6866 on 1 procs for 3000 steps with 2938 atoms
Performance: 26758.610 tau/day, 309.706 timesteps/s
99.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 | 8.3097 | 8.3097 | 8.3097 | 0.0 | 85.79
Neigh | 0.039037 | 0.039037 | 0.039037 | 0.0 | 0.40
Comm | 0.0021766 | 0.0021766 | 0.0021766 | 0.0 | 0.02
Output | 0.00048628 | 0.00048628 | 0.00048628 | 0.0 | 0.01
Modify | 1.3043 | 1.3043 | 1.3043 | 0.0 | 13.47
Other | | 0.0309 | | | 0.32
Nlocal: 2938 ave 2938 max 2938 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: 73242 ave 73242 max 73242 min
Histogram: 1 0 0 0 0 0 0 0 0 0
Total # of neighbors = 73242
Ave neighs/atom = 24.929204
Neighbor list builds = 26
Dangerous builds = 0
Total wall time: 0:00:09

View File

@ -0,0 +1,122 @@
LAMMPS (15 Sep 2022)
using 1 OpenMP thread(s) per MPI task
# Vesicle demo
variable T equal 0.2
variable LD equal 1.0
units lj
atom_style ellipsoid
boundary p p p
read_data read_data.vesicle1026
Reading data file ...
orthogonal box = (-35 -35 -35) to (35 35 35)
1 by 2 by 2 MPI processor grid
reading atoms ...
2938 atoms
2938 ellipsoids
read_data CPU = 0.137 seconds
compute ali all temp/com
compute rott all temp/asphere bias ali
velocity all create ${T} 87287 loop geom
velocity all create 0.2 87287 loop geom
pair_style ylz 2.6
pair_coeff * * 1.0 1.0 4 3 0.0 2.6
neighbor 1.0 bin
thermo_style custom step temp press pxx pyy
thermo 200
timestep 0.001
#dump 1 all atom 10 onlymembrane2.lammpstrj
fix 1 all langevin ${T} ${T} ${LD} 48279
fix 1 all langevin 0.2 ${T} ${LD} 48279
fix 1 all langevin 0.2 0.2 ${LD} 48279
fix 1 all langevin 0.2 0.2 1 48279
fix 2 all nve/asphere
run 3000
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- pair ylz command:
@Article{Yuan10,
author = {H. Yuan, C. Huang, J. Li, G. Lykotrafitis, and S. Zhang},
title = {One-particle-thick, solvent-free, coarse-grained model for biological and biomimetic fluid membranes},
journal = {Phys. Rev. E},
year = 2010,
volume = 82,
pages = {011905}
}
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Generated 0 of 1 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update: every = 1 steps, delay = 10 steps, check = yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 3.6
ghost atom cutoff = 3.6
binsize = 1.8, bins = 39 39 39
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair ylz, perpetual
attributes: half, newton on
pair build: half/bin/atomonly/newton
stencil: half/bin/3d
bin: standard
Per MPI rank memory allocation (min/avg/max) = 4.899 | 4.9 | 4.902 Mbytes
Step Temp Press Pxx Pyy
0 0.2 -0.0054891414 -0.0052713616 -0.0053641136
200 0.12893798 -0.0051492794 -0.0048734875 -0.0049624005
400 0.13798694 -0.004875313 -0.0047071897 -0.0049305051
600 0.14725193 -0.0046349542 -0.004719983 -0.0045791451
800 0.16146954 -0.0042232199 -0.0040789193 -0.0043672895
1000 0.17268468 -0.0037146703 -0.0036270393 -0.0039169034
1200 0.18266242 -0.0032749755 -0.0032971704 -0.0033323855
1400 0.18500165 -0.0028179031 -0.0030659821 -0.0027519633
1600 0.19513132 -0.0023407512 -0.0025109801 -0.0023416835
1800 0.19645259 -0.0019995412 -0.0019064341 -0.0021757747
2000 0.19658104 -0.0015897919 -0.0015850523 -0.0016828478
2200 0.1989936 -0.0011794062 -0.0011779716 -0.0012070706
2400 0.20011525 -0.0009147432 -0.00094040885 -0.001073922
2600 0.2013975 -0.00059253676 -0.00051920304 -0.00075138934
2800 0.19715513 -0.00020995605 -0.00024386426 -0.0005475745
3000 0.1976782 -0.0001308553 5.693004e-05 -0.00034478469
Loop time of 8.46954 on 4 procs for 3000 steps with 2938 atoms
Performance: 30603.800 tau/day, 354.211 timesteps/s
77.9% CPU use with 4 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 2.1751 | 2.4162 | 2.6707 | 13.2 | 28.53
Neigh | 0.0090503 | 0.0098389 | 0.010503 | 0.6 | 0.12
Comm | 3.5807 | 4.1526 | 4.9283 | 24.3 | 49.03
Output | 0.00032165 | 0.0029648 | 0.010842 | 8.4 | 0.04
Modify | 0.34554 | 0.37442 | 0.39198 | 2.8 | 4.42
Other | | 1.514 | | | 17.87
Nlocal: 734.5 ave 739 max 730 min
Histogram: 1 0 1 0 0 0 0 1 0 1
Nghost: 420.25 ave 424 max 415 min
Histogram: 1 0 0 0 0 1 0 1 0 1
Neighs: 18304.5 ave 19839 max 16636 min
Histogram: 1 0 0 0 1 0 1 0 0 1
Total # of neighbors = 73218
Ave neighs/atom = 24.921035
Neighbor list builds = 25
Dangerous builds = 0
Total wall time: 0:00:08

File diff suppressed because it is too large Load Diff

2
src/.gitignore vendored
View File

@ -1336,6 +1336,8 @@
/pair_yukawa_colloid.h
/pair_wf_cut.cpp
/pair_wf_cut.h
/pair_ylz.cpp
/pair_ylz.h
/pair_momb.cpp
/pair_momb.h
/pppm.cpp

509
src/ASPHERE/pair_ylz.cpp Normal file
View File

@ -0,0 +1,509 @@
/* ----------------------------------------------------------------------
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: Mehdi Baghaee (SUSTech)
------------------------------------------------------------------------- */
#include "pair_ylz.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
using namespace LAMMPS_NS;
using MathConst::MY_4PI;
using MathConst::MY_PI2;
using MathConst::MY_TWOBYSIXTH;
static const char cite_pair_ylz[] =
"pair ylz command:\n\n"
"@Article{Yuan10,\n"
" author = {H. Yuan, C. Huang, J. Li, G. Lykotrafitis, and S. Zhang},\n"
" title = {One-particle-thick, solvent-free, coarse-grained model for biological and "
"biomimetic fluid membranes},\n"
" journal = {Phys. Rev. E},\n"
" year = 2010,\n"
" volume = 82,\n"
" pages = {011905}\n"
"}\n\n";
/* ---------------------------------------------------------------------- */
PairYLZ::PairYLZ(LAMMPS *lmp) :
Pair(lmp), epsilon(nullptr), sigma(nullptr), cut(nullptr), zeta(nullptr), mu(nullptr),
beta(nullptr)
{
if (lmp->citeme) lmp->citeme->add(cite_pair_ylz);
single_enable = 0;
writedata = 1;
}
/* ----------------------------------------------------------------------
free all arrays
------------------------------------------------------------------------- */
PairYLZ::~PairYLZ()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(epsilon);
memory->destroy(sigma);
memory->destroy(cut);
memory->destroy(zeta);
memory->destroy(mu);
memory->destroy(beta);
}
}
/* ---------------------------------------------------------------------- */
void PairYLZ::compute(int eflag, int vflag)
{
int i, j, ii, jj, inum, jnum, itype, jtype;
double evdwl, one_eng, rsq, factor_lj;
double fforce[3], ttor[3], rtor[3], r12[3];
double a1[3][3], a2[3][3];
int *ilist, *jlist, *numneigh, **firstneigh;
double *iquat, *jquat;
evdwl = 0.0;
ev_init(eflag, vflag);
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
int *ellipsoid = atom->ellipsoid;
double **x = atom->x;
double **f = atom->f;
double **tor = atom->torque;
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];
itype = type[i];
iquat = bonus[ellipsoid[i]].quat;
MathExtra::quat_to_mat_trans(iquat, a1);
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
j &= NEIGHMASK;
// r12 = center to center vector
r12[0] = x[j][0] - x[i][0];
r12[1] = x[j][1] - x[i][1];
r12[2] = x[j][2] - x[i][2];
rsq = MathExtra::dot3(r12, r12);
jtype = type[j];
// compute if less than cutoff
if (rsq < cutsq[itype][jtype]) {
jquat = bonus[ellipsoid[j]].quat;
MathExtra::quat_to_mat_trans(jquat, a2);
one_eng = ylz_analytic(i, j, a1, a2, r12, rsq, fforce, ttor, rtor);
fforce[0] *= factor_lj;
fforce[1] *= factor_lj;
fforce[2] *= factor_lj;
ttor[0] *= factor_lj;
ttor[1] *= factor_lj;
ttor[2] *= factor_lj;
f[i][0] += fforce[0];
f[i][1] += fforce[1];
f[i][2] += fforce[2];
tor[i][0] += ttor[0];
tor[i][1] += ttor[1];
tor[i][2] += ttor[2];
if (newton_pair || j < nlocal) {
rtor[0] *= factor_lj;
rtor[1] *= factor_lj;
rtor[2] *= factor_lj;
f[j][0] -= fforce[0];
f[j][1] -= fforce[1];
f[j][2] -= fforce[2];
tor[j][0] += rtor[0];
tor[j][1] += rtor[1];
tor[j][2] += rtor[2];
}
if (eflag) evdwl = factor_lj * one_eng;
if (evflag)
ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fforce[0], fforce[1], fforce[2],
-r12[0], -r12[1], -r12[2]);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairYLZ::allocate()
{
allocated = 1;
int np1 = atom->ntypes + 1;
memory->create(setflag, np1, np1, "pair:setflag");
for (int i = 1; i < np1; i++)
for (int j = i; j < np1; j++) setflag[i][j] = 0;
memory->create(cutsq, np1, np1, "pair:cutsq");
memory->create(epsilon, np1, np1, "pair:epsilon");
memory->create(sigma, np1, np1, "pair:sigma");
memory->create(cut, np1, np1, "pair:cut");
memory->create(zeta, np1, np1, "pair:zeta");
memory->create(mu, np1, np1, "pair:mu");
memory->create(beta, np1, np1, "pair:beta");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairYLZ::settings(int narg, char **arg)
{
if (narg != 1) error->all(FLERR, "Illegal pair_style command");
cut_global = utils::numeric(FLERR, arg[0], false, lmp);
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairYLZ::coeff(int narg, char **arg)
{
if (narg < 8 || 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 epsilon_one = utils::numeric(FLERR, arg[2], false, lmp);
double sigma_one = utils::numeric(FLERR, arg[3], false, lmp);
double zeta_one = utils::numeric(FLERR, arg[4], false, lmp);
double mu_one = utils::numeric(FLERR, arg[5], false, lmp);
double beta_one = utils::numeric(FLERR, arg[6], false, lmp);
double cut_one = utils::numeric(FLERR, arg[7], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one;
cut[i][j] = cut_one;
zeta[i][j] = zeta_one;
mu[i][j] = mu_one;
beta[i][j] = beta_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairYLZ::init_style()
{
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
if (!avec) error->all(FLERR, "Pair style ylz requires atom style ellipsoid");
neighbor->request(this, instance_me);
// check that all atoms are ellipsoids
int flag_all, flag = 0;
const int *ellipsoid = atom->ellipsoid;
const int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; ++i)
if (ellipsoid[i] < 0) flag = 1;
MPI_Allreduce(&flag, &flag_all, 1, MPI_INT, MPI_MAX, world);
if (flag_all) error->all(FLERR, "All atoms must be ellipsoids for pair style ylz");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairYLZ::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i], epsilon[j][j], sigma[i][i], sigma[j][j]);
sigma[i][j] = mix_distance(sigma[i][i], sigma[j][j]);
cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
}
epsilon[j][i] = epsilon[i][j];
sigma[j][i] = sigma[i][j];
zeta[j][i] = zeta[i][j];
mu[j][i] = mu[i][j];
beta[j][i] = beta[i][j];
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairYLZ::write_restart(FILE *fp)
{
write_restart_settings(fp);
for (int i = 1; i <= atom->ntypes; i++) {
for (int j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) {
fwrite(&epsilon[i][j], sizeof(double), 1, fp);
fwrite(&sigma[i][j], sizeof(double), 1, fp);
fwrite(&cut[i][j], sizeof(double), 1, fp);
fwrite(&zeta[i][j], sizeof(double), 1, fp);
fwrite(&mu[i][j], sizeof(double), 1, fp);
fwrite(&beta[i][j], sizeof(double), 1, fp);
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairYLZ::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
for (int i = 1; i <= atom->ntypes; i++) {
for (int j = i; j <= atom->ntypes; j++) {
if (comm->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 (comm->me == 0) {
utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &zeta[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &mu[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &beta[i][j], sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&zeta[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&mu[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&beta[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairYLZ::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 PairYLZ::read_restart_settings(FILE *fp)
{
if (comm->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 PairYLZ::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp, "%d %g %g %g %g %g %g\n", i, epsilon[i][i], sigma[i][i], cut[i][i], zeta[i][i],
mu[i][i], beta[i][i]);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairYLZ::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, epsilon[i][i], sigma[i][i], cut[i][j],
zeta[i][j], mu[i][j], beta[i][j]);
}
/* ----------------------------------------------------------------------
compute analytic energy, force (fforce), and torque (ttor & rtor)
based on rotation matrices a
if newton is off, rtor is not calculated for ghost atoms
------------------------------------------------------------------------- */
double PairYLZ::ylz_analytic(const int i, const int j, double a1[3][3], double a2[3][3],
double *r12, const double rsq, double *fforce, double *ttor,
double *rtor)
{
int *type = atom->type;
int newton_pair = force->newton_pair;
int nlocal = atom->nlocal;
double r12hat[3];
MathExtra::normalize3(r12, r12hat);
double r = sqrt(rsq);
double ni1[3], nj1[3], dphi_drhat[3], dUdrhat[3], dUdni1[3], dUdnj1[3];
double dphi_dni1[3], dphi_dnj1[3];
double t, t1, t2, t4, cos_t, U, uR, uA, dUdr, dUdphi;
const double energy_well = epsilon[type[i]][type[j]];
const double rmin = MY_TWOBYSIXTH * sigma[type[i]][type[j]];
const double rcut = cut[type[i]][type[j]];
const double zt = zeta[type[i]][type[j]];
const double muu = mu[type[i]][type[j]];
const double sint = beta[type[i]][type[j]];
ni1[0] = a1[0][0];
ni1[1] = a1[0][1];
ni1[2] = a1[0][2];
nj1[0] = a2[0][0];
nj1[1] = a2[0][1];
nj1[2] = a2[0][2];
const double ninj = MathExtra::dot3(ni1, nj1);
const double ni1rhat = MathExtra::dot3(ni1, r12hat);
const double nj1rhat = MathExtra::dot3(nj1, r12hat);
const double a = ninj + (sint - ni1rhat) * (sint + nj1rhat) - 2.0 * sint * sint;
const double phi = 1.0 + (a - 1.0) * muu;
dphi_drhat[0] = muu * ((sint - ni1rhat) * nj1[0] - ni1[0] * (sint + nj1rhat));
dphi_drhat[1] = muu * ((sint - ni1rhat) * nj1[1] - ni1[1] * (sint + nj1rhat));
dphi_drhat[2] = muu * ((sint - ni1rhat) * nj1[2] - ni1[2] * (sint + nj1rhat));
dphi_dni1[0] = muu * (nj1[0] - r12hat[0] * (sint + nj1rhat));
dphi_dni1[1] = muu * (nj1[1] - r12hat[1] * (sint + nj1rhat));
dphi_dni1[2] = muu * (nj1[2] - r12hat[2] * (sint + nj1rhat));
dphi_dnj1[0] = muu * (ni1[0] + r12hat[0] * (sint - ni1rhat));
dphi_dnj1[1] = muu * (ni1[1] + r12hat[1] * (sint - ni1rhat));
dphi_dnj1[2] = muu * (ni1[2] + r12hat[2] * (sint - ni1rhat));
if (r < rmin) {
t = rmin / r;
t2 = t * t;
t4 = t2 * t2;
uR = (t4 - 2.0 * t2) * energy_well;
U = uR + (1.0 - phi) * energy_well;
dUdr = 4.0 * (t2 - t4) / r * energy_well;
dUdphi = -energy_well;
} else {
t = MY_PI2 * (r - rmin) / (rcut - rmin);
cos_t = cos(t);
t1 = cos_t;
for (int k = 1; k <= 2 * zt - 2; k++) t1 *= cos_t; // get cos()^(2zt-1)
uA = -energy_well * t1 * cos_t;
U = uA * phi;
dUdr = MY_4PI / (rcut - rmin) * (t1) *sin(t) * phi * energy_well;
dUdphi = uA;
}
dUdrhat[0] = dUdphi * dphi_drhat[0];
dUdrhat[1] = dUdphi * dphi_drhat[1];
dUdrhat[2] = dUdphi * dphi_drhat[2];
double dUdrhatrhat = MathExtra::dot3(dUdrhat, r12hat);
fforce[0] = dUdr * r12hat[0] + (dUdrhat[0] - dUdrhatrhat * r12hat[0]) / r;
fforce[1] = dUdr * r12hat[1] + (dUdrhat[1] - dUdrhatrhat * r12hat[1]) / r;
fforce[2] = dUdr * r12hat[2] + (dUdrhat[2] - dUdrhatrhat * r12hat[2]) / r;
// torque i
dUdni1[0] = dUdphi * dphi_dni1[0];
dUdni1[1] = dUdphi * dphi_dni1[1];
dUdni1[2] = dUdphi * dphi_dni1[2];
MathExtra::cross3(dUdni1, ni1, ttor); //minus sign is replace by swapping ni1 and dUdni1
if (newton_pair || j < nlocal) {
dUdnj1[0] = dUdphi * dphi_dnj1[0];
dUdnj1[1] = dUdphi * dphi_dnj1[1];
dUdnj1[2] = dUdphi * dphi_dnj1[2];
MathExtra::cross3(dUdnj1, nj1, rtor); //minus sign is replace by swapping ni1 and dUdni1
}
return U;
}

56
src/ASPHERE/pair_ylz.h Normal file
View File

@ -0,0 +1,56 @@
/* -*- 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(ylz,PairYLZ);
// clang-format on
#else
#ifndef LMP_PAIR_YLZ_H
#define LMP_PAIR_YLZ_H
#include "pair.h"
namespace LAMMPS_NS {
class PairYLZ : public Pair {
public:
PairYLZ(LAMMPS *lmp);
virtual ~PairYLZ();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
virtual void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
protected:
double cut_global;
double **epsilon, **sigma, **cut, **zeta, **mu, **beta;
class AtomVecEllipsoid *avec;
void allocate();
double ylz_analytic(const int i, const int j, double a1[3][3], double a2[3][3], double *r12,
const double rsq, double *fforce, double *ttor, double *rtor);
};
} // namespace LAMMPS_NS
#endif
#endif

View File

@ -30,7 +30,8 @@ namespace MathConst {
static constexpr double MY_ISPI4 = 1.12837916709551257389; // 1/sqrt(pi/4)
static constexpr double MY_SQRT2 = 1.41421356237309504880; // sqrt(2)
static constexpr double MY_ISQRT2 = 0.707106781186547524401; // 1/sqrt(2)
static constexpr double MY_CUBEROOT2 = 1.25992104989487316476; // 2*(1/3)
static constexpr double MY_CUBEROOT2 = 1.25992104989487316476; // 2^(1/3)
static constexpr double MY_TWOBYSIXTH = 1.12246204830937298142; // 2^(1/6) //
static constexpr double DEG2RAD = MY_PI / 180.0; // degree to radians
static constexpr double RAD2DEG = 180.0 / MY_PI; // radians to degree
} // namespace MathConst

View File

@ -49,7 +49,7 @@ else()
target_compile_definitions(style_tests PRIVATE TEST_INPUT_FOLDER=${TEST_INPUT_FOLDER} YAML_DECLARE_STATIC)
endif()
target_include_directories(style_tests PRIVATE ${LAMMPS_SOURCE_DIR})
target_link_libraries(style_tests PUBLIC gmock Yaml::Yaml lammps)
target_link_libraries(style_tests PUBLIC GTest::GMock Yaml::Yaml lammps)
# propagate sanitizer options to test tools
if(ENABLE_SANITIZER AND (NOT (ENABLE_SANITIZER STREQUAL "none")))
@ -62,7 +62,7 @@ endif()
# unit test for error stats class
add_executable(test_error_stats test_error_stats.cpp)
target_include_directories(test_error_stats PRIVATE ${CMAKE_CURRENT_SOURCE_DIR} ${LAMMPS_SOURCE_DIR})
target_link_libraries(test_error_stats PRIVATE gtest_main)
target_link_libraries(test_error_stats PRIVATE GTest::GMockMain)
add_test(NAME ErrorStats COMMAND test_error_stats)
# pair style tester