diff --git a/cmake/Modules/Packages/ML-PACE.cmake b/cmake/Modules/Packages/ML-PACE.cmake index ce8f02f5f4..248b8eea76 100644 --- a/cmake/Modules/Packages/ML-PACE.cmake +++ b/cmake/Modules/Packages/ML-PACE.cmake @@ -1,33 +1,40 @@ -set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.10.04.tar.gz" CACHE STRING "URL for PACE evaluator library sources") +set(PACELIB_URL "https://github.com/ICAMS/lammps-user-pace/archive/refs/tags/v.2023.11.25.fix.tar.gz" CACHE STRING "URL for PACE evaluator library sources") -set(PACELIB_MD5 "70ff79f4e59af175e55d24f3243ad1ff" CACHE STRING "MD5 checksum of PACE evaluator library tarball") +set(PACELIB_MD5 "b45de9a633f42ed65422567e3ce56f9f" CACHE STRING "MD5 checksum of PACE evaluator library tarball") mark_as_advanced(PACELIB_URL) mark_as_advanced(PACELIB_MD5) GetFallbackURL(PACELIB_URL PACELIB_FALLBACK) -# download library sources to build folder -if(EXISTS ${CMAKE_BINARY_DIR}/libpace.tar.gz) - file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) -endif() -if(NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}") - message(STATUS "Downloading ${PACELIB_URL}") - file(DOWNLOAD ${PACELIB_URL} ${CMAKE_BINARY_DIR}/libpace.tar.gz STATUS DL_STATUS SHOW_PROGRESS) - file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) - if((NOT DL_STATUS EQUAL 0) OR (NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}")) - message(WARNING "Download from primary URL ${PACELIB_URL} failed\nTrying fallback URL ${PACELIB_FALLBACK}") - file(DOWNLOAD ${PACELIB_FALLBACK} ${CMAKE_BINARY_DIR}/libpace.tar.gz EXPECTED_HASH MD5=${PACELIB_MD5} SHOW_PROGRESS) - endif() +# LOCAL_ML-PACE points to top-level dir with local lammps-user-pace repo, +# to make it easier to check local build without going through the public github releases +if(LOCAL_ML-PACE) + set(lib-pace "${LOCAL_ML-PACE}") else() - message(STATUS "Using already downloaded archive ${CMAKE_BINARY_DIR}/libpace.tar.gz") -endif() + # download library sources to build folder + if(EXISTS ${CMAKE_BINARY_DIR}/libpace.tar.gz) + file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) + endif() + if(NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}") + message(STATUS "Downloading ${PACELIB_URL}") + file(DOWNLOAD ${PACELIB_URL} ${CMAKE_BINARY_DIR}/libpace.tar.gz STATUS DL_STATUS SHOW_PROGRESS) + file(MD5 ${CMAKE_BINARY_DIR}/libpace.tar.gz DL_MD5) + if((NOT DL_STATUS EQUAL 0) OR (NOT "${DL_MD5}" STREQUAL "${PACELIB_MD5}")) + message(WARNING "Download from primary URL ${PACELIB_URL} failed\nTrying fallback URL ${PACELIB_FALLBACK}") + file(DOWNLOAD ${PACELIB_FALLBACK} ${CMAKE_BINARY_DIR}/libpace.tar.gz EXPECTED_HASH MD5=${PACELIB_MD5} SHOW_PROGRESS) + endif() + else() + message(STATUS "Using already downloaded archive ${CMAKE_BINARY_DIR}/libpace.tar.gz") + endif() -# uncompress downloaded sources -execute_process( - COMMAND ${CMAKE_COMMAND} -E remove_directory lammps-user-pace* - COMMAND ${CMAKE_COMMAND} -E tar xzf libpace.tar.gz - WORKING_DIRECTORY ${CMAKE_BINARY_DIR} -) -get_newest_file(${CMAKE_BINARY_DIR}/lammps-user-pace-* lib-pace) + + # uncompress downloaded sources + execute_process( + COMMAND ${CMAKE_COMMAND} -E remove_directory lammps-user-pace* + COMMAND ${CMAKE_COMMAND} -E tar xzf libpace.tar.gz + WORKING_DIRECTORY ${CMAKE_BINARY_DIR} + ) + get_newest_file(${CMAKE_BINARY_DIR}/lammps-user-pace-* lib-pace) +endif() add_subdirectory(${lib-pace} build-pace) set_target_properties(pace PROPERTIES CXX_EXTENSIONS ON OUTPUT_NAME lammps_pace${LAMMPS_MACHINE}) diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index dbd6b58ce7..07f5100ea4 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -115,6 +115,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`property/grid ` * :doc:`property/local ` * :doc:`ptm/atom ` + * :doc:`rattlers/atom ` * :doc:`rdf ` * :doc:`reduce ` * :doc:`reduce/chunk ` diff --git a/doc/src/Commands_fix.rst b/doc/src/Commands_fix.rst index 7301d1345e..3283561c41 100644 --- a/doc/src/Commands_fix.rst +++ b/doc/src/Commands_fix.rst @@ -122,6 +122,7 @@ OPT. * :doc:`mvv/tdpd ` * :doc:`neb ` * :doc:`neb/spin ` + * :doc:`nonaffine/displacement ` * :doc:`nph (ko) ` * :doc:`nph/asphere (o) ` * :doc:`nph/body ` diff --git a/doc/src/Howto_body.rst b/doc/src/Howto_body.rst index 19ae10b9a0..0588e079d0 100644 --- a/doc/src/Howto_body.rst +++ b/doc/src/Howto_body.rst @@ -355,7 +355,7 @@ faces are listed, so that M = 6 + 3\*N + 1. The integer line has three values: number of vertices (N), number of edges (E) and number of faces (F). The floating point line(s) list 6 moments of inertia followed by the coordinates of the N vertices (x1 -to zN) as 3N values, followed by 2N vertex indices corresponding to +to zN) as 3N values, followed by 2E vertex indices corresponding to the end points of the E edges, then 4\*F vertex indices defining F faces. The last value is the diameter value = the rounded diameter of the sphere that surrounds each vertex. The diameter value can be diff --git a/doc/src/compute.rst b/doc/src/compute.rst index 6737203618..a99d7d7c9c 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -279,6 +279,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`property/grid ` - convert per-grid attributes to per-grid vectors/arrays * :doc:`property/local ` - convert local attributes to local vectors/arrays * :doc:`ptm/atom ` - determines the local lattice structure based on the Polyhedral Template Matching method +* :doc:`rattlers/atom ` - identify under-coordinated rattler atoms * :doc:`rdf ` - radial distribution function :math:`g(r)` histogram of group of atoms * :doc:`reduce ` - combine per-atom quantities into a single global value * :doc:`reduce/chunk ` - reduce per-atom quantities within each chunk diff --git a/doc/src/compute_contact_atom.rst b/doc/src/compute_contact_atom.rst index 31aa24aa60..b7ed062ff6 100644 --- a/doc/src/compute_contact_atom.rst +++ b/doc/src/compute_contact_atom.rst @@ -36,6 +36,9 @@ sum of the radii of the two particles. The value of the contact number will be 0.0 for atoms not in the specified compute group. +The optional *group2-ID* argument allows to specify from which group atoms +contribute to the coordination number. Default setting is group 'all'. + Output info """"""""""" @@ -47,9 +50,6 @@ overview of LAMMPS output options. The per-atom vector values will be a number :math:`\ge 0.0`, as explained above. -The optional *group2-ID* argument allows to specify from which group atoms -contribute to the coordination number. Default setting is group 'all.' - Restrictions """""""""""" @@ -69,6 +69,3 @@ Default """"""" *group2-ID* = all - - -none diff --git a/doc/src/compute_rattlers_atom.rst b/doc/src/compute_rattlers_atom.rst new file mode 100644 index 0000000000..a69d091466 --- /dev/null +++ b/doc/src/compute_rattlers_atom.rst @@ -0,0 +1,92 @@ +.. index:: compute rattlers/atom + +compute rattlers/atom command +======================== + +Syntax +"""""" + +.. parsed-literal:: + + compute ID group-ID rattlers/atom cutoff zmin ntries + +* ID, group-ID are documented in :doc:`compute ` command +* rattlers/atom = style name of this compute command +* cutoff = *type* or *radius* + + .. parsed-literal:: + + *type* = cutoffs determined based on atom types + *radius* = cutoffs determined based on atom diameters (atom style sphere) + +* zmin = minimum coordination for a non-rattler atom +* ntries = maximum number of iterations to remove rattlers + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all rattlers/atom type 4 10 + +Description +""""""""""" + +.. versionadded:: TBD + +Define a compute that identifies rattlers in a system. Rattlers are often +identified in granular or glassy packings as undercoordinated atoms that +do not have the required number of contacts to constrain their translational +degrees of freedom. Such atoms are not considered rigid and can often freely +rattle around in the system. This compute identifies rattlers which can be +helpful for excluding them from analysis or providing extra damping forces +to accelerate relaxation processes. + +Rattlers are identified using an interactive approach. The coordination +number of all atoms is first calculated. The *type* and *radius* settings +are used to select whether interaction cutoffs are determined by atom +types or by the sum of atomic radii (atom style sphere), respectively. +Rattlers are then identified as atoms with a coordination number less +than *zmin* and are removed from consideration. Atomic coordination +numbers are then recalculated, excluding previously identified rattlers, +to identify a new set of rattlers. This process is iterated up to a maximum +of *ntries* or until no new rattlers are identified and the remaining +atoms form a stable network of contacts. + +In dense homogeneous systems where the average atom coordination number +is expected to be larger than *zmin*, this process usually only takes a few +iterations and a value of *ntries* around ten may be sufficient. In systems +with significant heterogeneity or average coordination numbers less than +*zmin*, an appropriate value of *ntries* depends heavily on the specific +system. For instance, a linear chain of N rattler atoms with a *zmin* of 2 +would take N/2 iterations to identify that all the atoms are rattlers. + +Output info +""""""""""" + +This compute calculates a per-atom vector and a global scalar. The vector +designates which atoms are rattlers, indicated by a value 1. Non-rattlers +have a value of 0. The global scalar returns the total number of rattlers +in the system. See the :doc:`Howto output ` page for an +overview of LAMMPS output options. + +Restrictions +"""""""""""" + +This compute is part of the EXTRA-COMPUTE package. It is only enabled if +LAMMPS was built with that package. See the +:doc:`Build package ` page for more info. + +The *radius* cutoff option requires that atoms store a radius as defined by the +:doc:`atom_style sphere ` or similar commands. + +Related commands +"""""""""""""""" + +:doc:`compute coord/atom ` +:doc:`compute contact/atom ` + +Default +""""""" + +none diff --git a/doc/src/fix.rst b/doc/src/fix.rst index 0889fe281f..69a7212487 100644 --- a/doc/src/fix.rst +++ b/doc/src/fix.rst @@ -287,6 +287,7 @@ accelerated styles exist. * :doc:`mvv/tdpd ` - constant temperature DPD using the modified velocity-Verlet algorithm * :doc:`neb ` - nudged elastic band (NEB) spring forces * :doc:`neb/spin ` - nudged elastic band (NEB) spring forces for spins +* :doc:`nonaffine/displacement ` - calculate nonaffine displacement of atoms * :doc:`nph ` - constant NPH time integration via Nose/Hoover * :doc:`nph/asphere ` - NPH for aspherical particles * :doc:`nph/body ` - NPH for body particles diff --git a/doc/src/fix_nonaffine_displacement.rst b/doc/src/fix_nonaffine_displacement.rst new file mode 100644 index 0000000000..363b0a747a --- /dev/null +++ b/doc/src/fix_nonaffine_displacement.rst @@ -0,0 +1,133 @@ +.. index:: fix nonaffine/displacement + +fix nonaffine/displacement command +================================== + +Syntax +"""""" + +.. parsed-literal:: + + fix ID group nonaffine/displacement style args reference/style nstep + +* ID, group are documented in :doc:`fix ` command +* nonaffine/displacement = style name of this fix command +* nevery = calculate nonaffine displacement every this many timesteps +* style = *d2min* or *integrated* + + .. parsed-literal:: + + *d2min* args = cutoff args + cutoff = *type* or *radius* or *custom* + *type* args = none, cutoffs determined by atom types + *radius* args = none, cutoffs determined based on atom diameters (atom style sphere) + *custom* args = *rmax*, cutoff set by a constant numeric value *rmax* (distance units) + *integrated* args = none + +* reference/style = *fixed* or *update* or *offset* + + .. parsed-literal:: + + *fixed* = use a fixed reference frame at *nstep* + *update* = update the reference frame every *nstep* timesteps + *offset* = update the reference frame *nstep* timesteps before calculating the nonaffine displacement + +Examples +"""""""" + +.. code-block:: LAMMPS + + fix 1 all nonaffine/displacement 100 integrated update 100 + fix 1 all nonaffine/displacement 1000 d2min type fixed 0 + fix 1 all nonaffine/displacement 1000 d2min custom 2.0 offset 100 + +Description +""""""""""" + +.. versionadded:: TBD + +This fix computes different metrics of the nonaffine displacement of +particles. The first metric, *d2min* calculates the :math:`D^2_\mathrm{min}` +nonaffine displacement by Falk and Langer in :ref:`(Falk) `. +For each atom, the fix computes the two tensors + +.. math:: + + X = \sum_{\mathrm{neighbors}} \vec{r} \left(\vec{r}_{0} \right)^T + +and + +.. math:: + + Y = \sum_{\mathrm{neighbors}} \vec{r}_0 \left(\vec{r}_{0} \right)^T + +where the neighbors include all other atoms within the distance criterion +set by the cutoff option, discussed below, :math:`\vec{r}` is the current +displacement between particles, and :math:`\vec{r}_0` is the reference +displacement. A deformation gradient tensor is then calculated as +:math:`F = X Y^{-1}` from which + +.. math:: + + D^2_\mathrm{min} = \sum_{\mathrm{neighbors}} \left| \vec{r} - F \vec{r}_0 \right|^2 + +and a strain tensor is calculated :math:`E = F F^{T} - I` where :math:`I` +is the identity tensor. This calculation is only performed on timesteps that +are a multiple of *nevery* (including timestep zero). Data accessed before +this occurs will simply be zeroed. + +The *integrated* style simply integrates the velocity of particles +every timestep to calculate a displacement. This style only works if +used in conjunction with another fix that deforms the box and displaces +atom positions such as :doc:`fix deform ` with remap x, +:doc:`fix press/berendsen `, or :doc:`fix nh `. + +Both of these methods require defining a reference state. With the *fixed* reference +style, the user picks a specific timestep *nstep* at which particle positions are saved. +If peratom data is accessed from this compute prior to this timestep, it will simply be +zeroed. The *update* reference style implies the reference state will be updated every +*nstep* timesteps. The *offset* reference only applies to the *d2min* metric and will +update the reference state *nstep* timesteps before a multiple of *nevery* timesteps. + + +---------- + +Restart, fix_modify, output, run start/stop, minimize info +""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" + +The reference state is saved to :doc:`binary restart files `. + +None of the :doc:`fix_modify ` options are relevant to this +fix. + +This fix computes a peratom array with 3 columns, which can be accessed +by indices 1-3 using any command that uses per-atom values from a fix +as input. + +For the *integrated* style, the three columns are the nonaffine +displacements in the x, y, and z directions. For the *d2min* style, +the three columns are the calculated :math:`\sqrt{D^2_\mathrm{min}}`, the +volumetric strain, and the deviatoric strain. + +Restrictions +"""""""""""" + +This compute is part of the EXTRA-FIX package. It is only enabled if +LAMMPS was built with that package. See the +:doc:`Build package ` page for more info. + +Related commands +"""""""""""""""" + +none + +Default +""""""" + +none + +---------- + +.. _d2min-Falk: + +**(Falk)** Falk and Langer PRE, 57, 7192 (1998). diff --git a/doc/src/fix_rigid.rst b/doc/src/fix_rigid.rst index 2192f9ed1f..0277793d0f 100644 --- a/doc/src/fix_rigid.rst +++ b/doc/src/fix_rigid.rst @@ -80,7 +80,7 @@ Syntax groupID1, groupID2, ... = list of N group IDs * zero or more keyword/value pairs may be appended -* keyword = *langevin* or *reinit* or *temp* or *iso* or *aniso* or *x* or *y* or *z* or *couple* or *tparam* or *pchain* or *dilate* or *force* or *torque* or *infile* or *gravity* +* keyword = *langevin* or *reinit* or *temp* or *mol* or *iso* or *aniso* or *x* or *y* or *z* or *couple* or *tparam* or *pchain* or *dilate* or *force* or *torque* or *infile* or *gravity* .. parsed-literal:: @@ -92,6 +92,8 @@ Syntax *temp* values = Tstart Tstop Tdamp Tstart,Tstop = desired temperature at start/stop of run (temperature units) Tdamp = temperature damping parameter (time units) + *mol* value = template-ID + template-ID = ID of molecule template specified in a separate :doc:`molecule ` command *iso* or *aniso* values = Pstart Pstop Pdamp Pstart,Pstop = scalar external pressure at start/end of run (pressure units) Pdamp = pressure damping parameter (time units) diff --git a/doc/src/molecule.rst b/doc/src/molecule.rst index 480e175e7b..1de905886c 100644 --- a/doc/src/molecule.rst +++ b/doc/src/molecule.rst @@ -154,21 +154,25 @@ These are the recognized header keywords. Header lines can come in any order. The numeric value(s) are read from the beginning of the line. The keyword should appear at the end of the line. All these settings have default values, as explained below. A line need only -appear if the value(s) are different than the default. +appear if the value(s) are different than the default, except when +defining a *body* particle, which requires setting the number of +*atoms* to 1, and setting the *inertia* in a specific section (see below). * N *atoms* = # of atoms N in molecule, default = 0 * Nb *bonds* = # of bonds Nb in molecule, default = 0 * Na *angles* = # of angles Na in molecule, default = 0 * Nd *dihedrals* = # of dihedrals Nd in molecule, default = 0 * Ni *impropers* = # of impropers Ni in molecule, default = 0 -* Nf *fragments* = # of fragments in molecule, default = 0 +* Nf *fragments* = # of fragments Nf in molecule, default = 0 +* Ninteger Ndouble *body* = # of integer and floating-point values in body +particle, default = 0 * Mtotal *mass* = total mass of molecule * Xc Yc Zc *com* = coordinates of center-of-mass of molecule * Ixx Iyy Izz Ixy Ixz Iyz *inertia* = 6 components of inertia tensor of molecule For *mass*, *com*, and *inertia*, the default is for LAMMPS to calculate this quantity itself if needed, assuming the molecules -consists of a set of point particles or finite-size particles (with a +consist of a set of point particles or finite-size particles (with a non-zero diameter) that do not overlap. If finite-size particles in the molecule do overlap, LAMMPS will not account for the overlap effects when calculating any of these 3 quantities, so you should @@ -188,6 +192,7 @@ These are the allowed section keywords for the body of the file. * *Bonds, Angles, Dihedrals, Impropers* = molecular topology sections * *Special Bond Counts, Special Bonds* = special neighbor info * *Shake Flags, Shake Atoms, Shake Bond Types* = SHAKE info +* *Body Integers, Body Doubles* = body-property sections For the Types, Bonds, Angles, Dihedrals, and Impropers sections, each atom/bond/angle/etc type can be specified either as a number (numeric @@ -515,6 +520,67 @@ of SHAKE clusters. ---------- +*Body Integers* section: + +* one line +* line syntax: N E F +* N = number of sub-particles or number or vertices +* E,F = number of edges and faces + +This section is only needed when the molecule is a body particle. the other +Body section must also appear in the file. + +The total number of values that must appear is determined by the body style, and +must be equal to the Ninteger value given in the *body* header. + +For *nparticle* and *rounded/polygon*, only the number of sub-particles or +vertices N is required, and Ninteger should have a value of 1. + +For *rounded/polyhedron*, the number of edges E and faces F is required, and +Ninteger should have a value of 3. + +See the :doc:`Howto body ` page for a further description of +the file format. + +---------- + +*Body Doubles* section: + +* first line +* line syntax: Ixx Iyy Izz Ixy Ixz Iyz +* Ixx Iyy Izz Ixy Ixz Iyz = 6 components of inertia tensor of body particle +* one line per sub-particle or vertex +* line syntax: x y z +* x, y, z = coordinates of sub-particle or vertex +* one line per edge +* line syntax: N1 N2 +* N1, N2 = vertex indices +* one line per face +* line syntax: N1 N2 N3 N4 +* N1, N2, N3, N4 = vertex indices +* last line +* line syntax: diam +* diam = rounded diameter that surrounds each vertex + +This section is only needed when the molecule is a body particle. the other +Body section must also appear in the file. + +The total number of values that must appear is determined by the body style, and +must be equal to the Ndouble value given in the *body* header. The 6 moments of +inertia and the 3N coordinates of the sub-particles or vertices are required +for all body styles. + +For *rounded/polygon*, the E = 6 + 3*N + 1 edges are automatically determined +from the vertices. + +For *rounded/polyhedron*, the 2E vertex indices for the end points of the edges +and 4F vertex indices defining the faces are required. + +See the :doc:`Howto body ` page for a further description of +the file format. + +---------- + Restrictions """""""""""" diff --git a/doc/src/pair_pace.rst b/doc/src/pair_pace.rst index d815f663fe..001214370c 100644 --- a/doc/src/pair_pace.rst +++ b/doc/src/pair_pace.rst @@ -40,6 +40,9 @@ Examples pair_style pace product chunksize 2048 pair_coeff * * Cu-PBE-core-rep.ace Cu + pair_style pace + pair_coeff * * Cu.yaml Cu + pair_style pace/extrapolation pair_coeff * * Cu.yaml Cu.asi Cu @@ -64,7 +67,7 @@ specifies an ACE coefficient file followed by N additional arguments specifying the mapping of ACE elements to LAMMPS atom types, where N is the number of LAMMPS atom types: -* ACE coefficient file +* ACE coefficient file (.yaml or .yace/.ace format) * N element names = mapping of ACE elements to atom types Only a single pair_coeff command is used with the *pace* style which @@ -136,6 +139,22 @@ product B-basis evaluator is always used and only *linear* ASI is supported. See the :doc:`pair_coeff ` page for alternate ways to specify the path for the ACE coefficient file. +Core repulsion +""""""""""""""""""" +The ACE potential can be configured to initiate core-repulsion from an inner cutoff, +seamlessly transitioning from ACE to ZBL. The core repulsion factor can be accessed +as a per-atom quantity, as demonstrated in the example below: + +.. code-block:: LAMMPS + + pair_style pace + pair_coeff * * CuNi.yaml Cu Ni + + fix pace_corerep all pair 1 pace corerep 1 + +In this case, per-atom `f_pace_corerep` quantities represent the fraction of ZBL +core-repulsion for each atom. + Mixing, shift, table, tail correction, restart, rRESPA info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 072282df4d..63564fe2e3 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1087,6 +1087,7 @@ facesets factorizable factorizations Fahrenberger +Falk Faken Farago Fasolino @@ -1835,6 +1836,7 @@ Lanczos Lande Landron Landsgesell +Langer langevin Langevin Langston @@ -2512,6 +2514,7 @@ noforce noguess Noid nolib +nonaffine nonequilibrium nongauss nonGaussian @@ -2584,6 +2587,7 @@ nthreads ntimestep Ntptask Ntriples +ntries ntris Ntype ntypes diff --git a/lib/pace/Install.py b/lib/pace/Install.py index 8d31852e44..fcd9497937 100644 --- a/lib/pace/Install.py +++ b/lib/pace/Install.py @@ -18,11 +18,11 @@ from install_helpers import fullpath, geturl, checkmd5sum, getfallback # settings thisdir = fullpath('.') -version ='v.2023.10.04' +version ='v.2023.11.25.fix' # known checksums for different PACE versions. used to validate the download. checksums = { \ - 'v.2023.10.04': '70ff79f4e59af175e55d24f3243ad1ff' + 'v.2023.11.25.fix': 'b45de9a633f42ed65422567e3ce56f9f' } parser = ArgumentParser(prog='Install.py', description="LAMMPS library build wrapper script") diff --git a/src/.gitignore b/src/.gitignore index 3ee771e139..e2c9fa6d5b 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -629,6 +629,8 @@ /compute_pressure_grem.h /compute_ptm_atom.cpp /compute_ptm_atom.h +/compute_rattlers_atom.cpp +/compute_rattlers_atom.h /compute_rigid_local.cpp /compute_rigid_local.h /compute_smd_triangle_vertices.cpp @@ -912,6 +914,8 @@ /fix_nvt_sllod_eff.h /fix_nve_tri.cpp /fix_nve_tri.h +/fix_nonaffine_displacement.cpp +/fix_nonaffine_displacement.h /fix_oneway.cpp /fix_oneway.h /fix_orient_bcc.cpp diff --git a/src/EXTRA-COMPUTE/compute_rattlers_atom.cpp b/src/EXTRA-COMPUTE/compute_rattlers_atom.cpp new file mode 100644 index 0000000000..602923b58a --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_rattlers_atom.cpp @@ -0,0 +1,312 @@ +// clang-format off +/* ---------------------------------------------------------------------- + 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 authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL) +------------------------------------------------------------------------- */ + +#include "compute_rattlers_atom.h" + +#include "atom.h" +#include "comm.h" +#include "error.h" +#include "force.h" +#include "memory.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "pair.h" +#include "update.h" + +#include +#include + +using namespace LAMMPS_NS; + +enum { TYPE, RADIUS }; + +/* ---------------------------------------------------------------------- */ + +ComputeRattlersAtom::ComputeRattlersAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), ncontacts(nullptr), rattler(nullptr) +{ + if (narg != 6) error->all(FLERR, "Illegal compute rattlers/atom command"); + + if (strcmp(arg[3], "type") == 0) + cutstyle = TYPE; + else if (strcmp(arg[3], "radius") == 0) + cutstyle = RADIUS; + else + error->all(FLERR, "Illegal compute rattlers/atom command"); + + if (cutstyle == RADIUS && !atom->radius_flag) + error->all(FLERR, "Compute rattlers/atom radius style requires atom attribute radius"); + + ncontacts_rattler = utils::inumeric(FLERR, arg[4], false, lmp); + max_tries = utils::inumeric(FLERR, arg[5], false, lmp); + + nmax = 0; + invoked_peratom = -1; + + scalar_flag = 1; + extscalar = 1; + peratom_flag = 1; + size_peratom_cols = 0; + comm_forward = 1; + comm_reverse = 1; +} + +/* ---------------------------------------------------------------------- */ + +ComputeRattlersAtom::~ComputeRattlersAtom() +{ + memory->destroy(ncontacts); + memory->destroy(rattler); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlersAtom::init() +{ + if (force->pair == nullptr) error->all(FLERR, "No pair style is defined for compute rattlers"); + + // Cannot calculate distance from radii for JKR/DMT + if (force->pair->beyond_contact) + error->all(FLERR, "Compute rattlers does not currently support pair styles that extend beyond contact"); + + // need an occasional half neighbor list + // set size to same value as request made by force->pair + // this should enable it to always be a copy list (e.g. for granular pstyle) + + auto pairrequest = neighbor->find_request(force->pair); + if (pairrequest && pairrequest->get_size()) + neighbor->add_request(this, NeighConst::REQ_SIZE | NeighConst::REQ_OCCASIONAL); + else + neighbor->add_request(this, NeighConst::REQ_OCCASIONAL); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlersAtom::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlersAtom::compute_peratom() +{ + if (invoked_peratom == update->ntimestep) return; + invoked_peratom = update->ntimestep; + + int i, j, ii, jj, inum, jnum, itype, jtype, tmp_flag; + tagint itag, jtag; + double xtmp, ytmp, ztmp, delx, dely, delz; + double rsq, radsum; + + if (nmax < atom->nmax) { + nmax = atom->nmax; + memory->destroy(ncontacts); + memory->destroy(rattler); + memory->create(ncontacts, nmax, "rattlers:ncontacts"); + memory->create(rattler, nmax, "rattlers:rattler"); + vector_atom = rattler; + } + + for (i = 0; i < nmax; i++) rattler[i] = 0; + + int *ilist, *jlist, *numneigh, **firstneigh; + + double **x = atom->x; + double *radius = atom->radius; + tagint *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + // invoke half neighbor list (will copy or build if necessary) + neighbor->build_one(list); + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + Pair *pair = force->pair; + double **cutsq = force->pair->cutsq; + + int change_flag = 1; + int ntry = 0; + while (ntry < max_tries) { + change_flag = 0; + + for (i = 0; i < nmax; i++) ncontacts[i] = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + if (rattler[i] == 1) continue; + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itag = tag[i]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + if (rattler[j] == 1) continue; + + // itag = jtag is possible for long cutoffs that include images of self + + if (newton_pair == 0 && j >= nlocal) { + jtag = tag[j]; + if (itag > jtag) { + if ((itag + jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag + jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + } + + jtype = type[j]; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + + if (cutstyle == TYPE) { + if (rsq >= cutsq[itype][jtype]) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq >= radsum * radsum) continue; + } + ncontacts[i] += 1; + if (newton_pair || j < nlocal) + ncontacts[j] += 1; + } + } + + // add contributions from ghosts + if (force->newton_pair) comm->reverse_comm(this); + + // Set flags for rattlers + for (i = 0; i < atom->nlocal; i++) { + if (ncontacts[i] < ncontacts_rattler && rattler[i] == 0) { + rattler[i] = 1; + change_flag = 1; + } + } + + comm->forward_comm(this); + + MPI_Allreduce(&change_flag, &tmp_flag, 1, MPI_INT, MPI_MAX, world); + change_flag = tmp_flag; + if (change_flag == 0) break; + + ntry += 1; + } + + if (change_flag == 1) + error->warning(FLERR, "Rattler calculation failed to converge within max tries"); +} + +/* ---------------------------------------------------------------------- */ + +double ComputeRattlersAtom::compute_scalar() +{ + if (invoked_peratom != update->ntimestep) + compute_peratom(); + + invoked_scalar = update->ntimestep; + + double total_rattlers = 0; + for (int i = 0; i < atom->nlocal; i++) { + if (rattler[i] == 1) { + total_rattlers += 1; + } + } + + //Total across processors + MPI_Allreduce(&total_rattlers, &scalar, 1, MPI_DOUBLE, MPI_SUM, world); + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRattlersAtom::pack_reverse_comm(int n, int first, double *buf) +{ + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = ubuf(ncontacts[i]).d; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlersAtom::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i, j, m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + ncontacts[j] += (int) ubuf(buf[m++]).i; + } +} + +/* ---------------------------------------------------------------------- */ + +int ComputeRattlersAtom::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i, j, m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = rattler[j]; + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeRattlersAtom::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m, last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + rattler[i] = buf[m++]; + } +} diff --git a/src/EXTRA-COMPUTE/compute_rattlers_atom.h b/src/EXTRA-COMPUTE/compute_rattlers_atom.h new file mode 100644 index 0000000000..257bae8374 --- /dev/null +++ b/src/EXTRA-COMPUTE/compute_rattlers_atom.h @@ -0,0 +1,52 @@ +/* -*- 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 COMPUTE_CLASS +// clang-format off +ComputeStyle(rattlers/atom,ComputeRattlersAtom); +// clang-format on +#else + +#ifndef LMP_COMPUTE_RATTLERS_ATOM_H +#define LMP_COMPUTE_RATTLERS_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeRattlersAtom : public Compute { + public: + ComputeRattlersAtom(class LAMMPS *, int, char **); + ~ComputeRattlersAtom() override; + void init() override; + void init_list(int, class NeighList *) override; + void compute_peratom() override; + double compute_scalar() override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; + + private: + int pstyle, cutstyle; + int ncontacts_rattler, max_tries, nmax, invoked_peratom; + int *ncontacts; + double *rattler; + class NeighList *list; + +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.cpp b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp new file mode 100644 index 0000000000..ef5481601f --- /dev/null +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.cpp @@ -0,0 +1,736 @@ +// clang-format off +/* ---------------------------------------------------------------------- + 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 authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL) +------------------------------------------------------------------------- */ + +#include "fix_nonaffine_displacement.h" + +#include "atom.h" +#include "citeme.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "fix_store_atom.h" +#include "force.h" +#include "group.h" +#include "math_extra.h" +#include "memory.h" +#include "modify.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "neighbor.h" +#include "pair.h" +#include "update.h" + +#include + +using namespace LAMMPS_NS; +using namespace FixConst; +using namespace MathExtra; + +enum { TYPE, RADIUS, CUSTOM }; +enum { INTEGRATED, D2MIN }; +enum { FIXED, OFFSET, UPDATE }; + +static const char cite_nonaffine_d2min[] = + "@article{PhysRevE.57.7192,\n" + " title = {Dynamics of viscoplastic deformation in amorphous solids},\n" + " author = {Falk, M. L. and Langer, J. S.},\n" + " journal = {Phys. Rev. E},\n" + " volume = {57},\n" + " issue = {6},\n" + " pages = {7192--7205},\n" + " numpages = {0},\n" + " year = {1998},\n" + " month = {Jun},\n" + " publisher = {American Physical Society},\n" + " doi = {10.1103/PhysRevE.57.7192},\n" + "url = {https://link.aps.org/doi/10.1103/PhysRevE.57.7192}\n" + "}\n\n"; + +/* ---------------------------------------------------------------------- */ + +FixNonaffineDisplacement::FixNonaffineDisplacement(LAMMPS *lmp, int narg, char **arg) : + Fix(lmp, narg, arg), id_fix(nullptr), X(nullptr), Y(nullptr), F(nullptr), norm(nullptr) +{ + if (narg < 4) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + + nevery = utils::inumeric(FLERR, arg[3], false, lmp); + if (nevery <= 0) error->all(FLERR,"Illegal nevery value {} in fix nonaffine/displacement", nevery); + + int iarg = 4; + if (strcmp(arg[iarg], "integrated") == 0) { + nad_style = INTEGRATED; + nevery = 1; + iarg += 1; + } else if (strcmp(arg[iarg], "d2min") == 0) { + if (iarg + 1 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + nad_style = D2MIN; + if (strcmp(arg[iarg + 1], "type") == 0) { + cut_style = TYPE; + } else if (strcmp(arg[iarg + 1], "radius") == 0) { + cut_style = RADIUS; + } else if (strcmp(arg[iarg + 1], "custom") == 0) { + if (iarg + 2 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + cut_style = CUSTOM; + cutoff_custom = utils::numeric(FLERR, arg[iarg + 2], false, lmp); + cutsq_custom = cutoff_custom * cutoff_custom; + if (cutoff_custom <= 0) + error->all(FLERR, "Illegal custom cutoff length {}", arg[iarg + 2]); + iarg += 1; + } else error->all(FLERR,"Illegal cutoff style {} in fix nonaffine/displacement", arg[iarg + 1]); + iarg += 2; + } else error->all(FLERR,"Illegal nonaffine displacement style {} in fix nonaffine/displacement", arg[iarg]); + + if (iarg + 2 > narg) error->all(FLERR,"Illegal fix nonaffine/displacement command"); + if (strcmp(arg[iarg], "fixed") == 0) { + reference_style = FIXED; + reference_timestep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if (update_timestep < 0) + error->all(FLERR, "Illegal reference timestep {} in fix nonaffine/displacement", arg[iarg + 1]); + } else if (strcmp(arg[iarg], "update") == 0) { + reference_style = UPDATE; + update_timestep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if (update_timestep < 0) + error->all(FLERR, "Illegal update timestep {} in fix nonaffine/displacement", arg[iarg + 1]); + } else if (strcmp(arg[iarg], "offset") == 0) { + reference_style = OFFSET; + offset_timestep = utils::inumeric(FLERR, arg[iarg + 1], false, lmp); + if ((offset_timestep <= 0) || (offset_timestep > nevery)) + error->all(FLERR, "Illegal offset timestep {} in fix nonaffine/displacement", arg[iarg + 1]); + } else error->all(FLERR,"Illegal reference style {} in fix nonaffine/displacement", arg[iarg]); + + if (nad_style == D2MIN) + if (cut_style == RADIUS && (!atom->radius_flag)) + error->all(FLERR, "Fix nonaffine/displacement radius style requires atom attribute radius"); + + if (nad_style == INTEGRATED && reference_style == OFFSET) + error->all(FLERR, "Fix nonaffine/displacement cannot use the integrated style with an offset reference state"); + + peratom_flag = 1; + peratom_freq = nevery; + nmax = -1; + reference_saved = 0; + restart_global = 1; + + size_peratom_cols = 3; + comm_reverse = 0; + comm_forward = 0; + if (nad_style == D2MIN) { + comm_reverse = 18; + comm_forward = 9; + } + + if (nad_style == D2MIN && lmp->citeme) lmp->citeme->add(cite_nonaffine_d2min); +} + +/* ---------------------------------------------------------------------- */ + +FixNonaffineDisplacement::~FixNonaffineDisplacement() +{ + if (id_fix && modify->nfix) modify->delete_fix(id_fix); + delete[] id_fix; + + if (nad_style == D2MIN) { + memory->destroy(X); + memory->destroy(Y); + memory->destroy(F); + memory->destroy(norm); + memory->destroy(array_atom); + } +} + +/* ---------------------------------------------------------------------- */ + +int FixNonaffineDisplacement::setmask() +{ + int mask = 0; + mask |= POST_FORCE; + return mask; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::post_constructor() +{ + // Create persistent peratom storage for either an integrated velocity or reference position + // Ghost atoms need reference coordinates for D2min + std::string ghost_status = "0"; + if (nad_style == D2MIN) ghost_status = "1"; + + id_fix = utils::strdup(id + std::string("_FIX_PA")); + fix = dynamic_cast(modify->add_fix(fmt::format("{} {} STORE/ATOM 3 0 {} 1", id_fix, group->names[igroup], ghost_status))); + + if (nad_style == INTEGRATED) + array_atom = fix->astore; + + if (nad_style == D2MIN) + grow_arrays(atom->nmax); + + for (int i = 0; i < atom->nlocal; i++) + for (int j = 0; j < 3; j++) array_atom[i][j] = 0.0; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::init() +{ + dtv = update->dt; + + if ((!reference_saved) && (reference_style == FIXED) && (update->ntimestep > reference_timestep)) + error->all(FLERR, "Initial timestep exceeds that of the reference state in fix nonaffine/displacement"); + + if (nad_style == D2MIN) { + if ((!force->pair) && (cut_style == TYPE)) + error->all(FLERR,"Fix nonaffine/displacement D2Min option requires a pair style be defined " + "or cutoff specified"); + + // need an occasional half neighbor list + + if (cut_style == RADIUS) { + auto req = neighbor->add_request(this, NeighConst::REQ_SIZE | NeighConst::REQ_OCCASIONAL); + } else { + auto req = neighbor->add_request(this, NeighConst::REQ_OCCASIONAL); + if (cut_style == CUSTOM) { + double skin = neighbor->skin; + mycutneigh = cutoff_custom + skin; + + double cutghost; // as computed by Neighbor and Comm + if (force->pair) + cutghost = MAX(force->pair->cutforce + skin, comm->cutghostuser); + else + cutghost = comm->cutghostuser; + + if (mycutneigh > cutghost) + error->all(FLERR,"Fix nonaffine/displacement D2Min option cutoff exceeds ghost atom range - use comm_modify cutoff command"); + + req->set_cutoff(mycutneigh); + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::init_list(int /*id*/, NeighList *ptr) +{ + list = ptr; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::setup(int vflag) +{ + post_force(0); // Save state if needed before starting the 1st timestep +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::post_force(int /*vflag*/) +{ + if (reference_saved && (!update->setupflag)) { + if (nad_style == INTEGRATED) { + integrate_velocity(); + } else { + if ((update->ntimestep % nevery) == 0) calculate_D2Min(); + } + } + + if (reference_style == FIXED) + if (update->ntimestep == reference_timestep) + save_reference_state(); + + if (reference_style == UPDATE) + if ((update->ntimestep % update_timestep) == 0) + save_reference_state(); + + if (reference_style == OFFSET) + if (((update->ntimestep + offset_timestep) % nevery) == 0) + save_reference_state(); +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::write_restart(FILE *fp) +{ + if (comm->me == 0) { + int size = sizeof(int); + fwrite(&size, sizeof(int), 1, fp); + fwrite(&reference_saved, sizeof(int), 1, fp); + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::restart(char *buf) +{ + reference_saved = (int) ubuf(buf[0]).i; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::integrate_velocity() +{ + int i,n; + dtv = update->dt; + + double **v = atom->v; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int m = 0; m < 3; m++) { + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) { + array_atom[i][m] += dtv * v[i][m]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::save_reference_state() +{ + int i, n; + double **x = atom->x; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + int nall = nlocal + atom->nghost; + + if (nad_style == D2MIN) { + for (int m = 0; m < 3; m++) { + for (int i = 0; i < nall; i++) { + if (mask[i] & groupbit) array_atom[i][m] = x[i][m]; + } + } + } else { + for (int m = 0; m < 3; m++) { + for (int i = 0; i < nall; i++) { + if (mask[i] & groupbit) array_atom[i][m] = 0.0; + } + } + } + + if (nad_style == D2MIN) { + xprd0 = domain->xprd; + yprd0 = domain->yprd; + zprd0 = domain->zprd; + xprd0_half = domain->xprd_half; + yprd0_half = domain->yprd_half; + zprd0_half = domain->zprd_half; + xy0 = domain->xy; + xz0 = domain->xz; + yz0 = domain->yz; + } + + reference_saved = 1; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::calculate_D2Min() +{ + // invoke half neighbor list (will copy or build if necessary) + neighbor->build_one(list); + + if (atom->nmax > nmax) + grow_arrays(atom->nmax); + + int i, j, k, l, ii, jj, inum, jnum, itype, jtype; + double evol, j2, edev; + double r[3], r0[3], rsq, rsq0, radsum, temp[3]; + double X_tmp[3][3], Y_tmp[3][3], F_tmp[3][3], E[3][3]; + double Y_inv[3][3] = {0.0}; // Zero for 2d since not all entries used + int *ilist, *jlist, *numneigh, **firstneigh; + + double **x = atom->x; + double **x0 = array_atom; + double *radius = atom->radius; + tagint *tag = atom->tag; + int *type = atom->type; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + int dim = domain->dimension; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + Pair *pair = force->pair; + double **cutsq; + if (pair) cutsq = force->pair->cutsq; + + for (i = 0; i < nmax; i++) { + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l++) { + X[i][k][l] = 0.0; + Y[i][k][l] = 0.0; + } + } + norm[i] = 0; + array_atom[i][0] = 0; + } + + // First loop through neighbors + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + jtype = type[j]; + r[0] = x[i][0] - x[j][0]; + r[1] = x[i][1] - x[j][1]; + r[2] = x[i][2] - x[j][2]; + rsq = lensq3(r); + + // Only include contributions from atoms that are CURRENTLY neighbors + if (cut_style == TYPE) { + if (rsq > cutsq[itype][jtype]) continue; + } else if (cut_style == CUSTOM) { + if (rsq > cutsq_custom) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq > (radsum * radsum)) continue; + } + + r0[0] = x0[i][0] - x0[j][0]; + r0[1] = x0[i][1] - x0[j][1]; + r0[2] = x0[i][2] - x0[j][2]; + minimum_image0(r0); + + // Using notation from Falk & Langer 1998 + outer3(r, r0, X_tmp); + outer3(r0, r0, Y_tmp); + + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l++) { + X[i][k][l] += X_tmp[k][l]; + Y[i][k][l] += Y_tmp[k][l]; + } + } + + if (newton_pair || j < nlocal) { + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l++) { + X[j][k][l] += X_tmp[k][l]; + Y[j][k][l] += Y_tmp[k][l]; + } + } + } + } + } + + comm_flag = 0; + if (newton_pair) comm->reverse_comm(this, 18); + + // Calculate contributions to strain tensor + double denom; + for (i = 0; i < nlocal; i++) { + if (!(mask[i] & groupbit)) continue; + for (j = 0; j < 3; j++) { + for (k = 0; k < 3; k++) { + Y_tmp[j][k] = Y[i][j][k]; + X_tmp[j][k] = X[i][j][k]; + } + } + + if (dim == 3) { + invert3(Y_tmp, Y_inv); + } else { + denom = Y_tmp[0][0] * Y_tmp[1][1] - Y_tmp[0][1] * Y_tmp[1][0]; + if (denom != 0.0) denom = 1.0 / denom; + Y_inv[0][0] = Y_tmp[1][1] * denom; + Y_inv[0][1] = -Y_tmp[0][1] * denom; + Y_inv[1][0] = -Y_tmp[1][0] * denom; + Y_inv[1][1] = Y_tmp[0][0] * denom; + } + + times3(X_tmp, Y_inv, F_tmp); + + for (j = 0; j < 3; j++) { + for (k = 0; k < 3; k++) { + F[i][j][k] = F_tmp[j][k]; + } + } + } + + comm->forward_comm(this); + + // Second loop through neighbors + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (!(mask[i] & groupbit)) continue; + + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + if (!(mask[j] & groupbit)) continue; + + jtype = type[j]; + r[0] = x[i][0] - x[j][0]; + r[1] = x[i][1] - x[j][1]; + r[2] = x[i][2] - x[j][2]; + rsq = lensq3(r); + + // Only include contributions from atoms that are CURRENTLY neighbors + if (cut_style == TYPE) { + if (rsq >= cutsq[itype][jtype]) continue; + } else if (cut_style == CUSTOM) { + if (rsq >= cutsq_custom) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq >= radsum * radsum) continue; + } + + r0[0] = x0[i][0] - x0[j][0]; + r0[1] = x0[i][1] - x0[j][1]; + r0[2] = x0[i][2] - x0[j][2]; + minimum_image0(r0); + + // E * r0 + for (k = 0; k < 3; k++) { + temp[k] = 0.0; + for (l = 0; l < 3; l++) + temp[k] += F[i][k][l] * r0[l]; + } + + sub3(r, temp, temp); + array_atom[i][0] += lensq3(temp); + norm[i] += 1; + + if (newton_pair || j < nlocal) { + for (k = 0; k < 3; k++) { + temp[k] = 0.0; + for (l = 0; l < 3; l++) + temp[k] += F[j][k][l] * r0[l]; + } + + sub3(r, temp, temp); + array_atom[j][0] += lensq3(temp); + norm[j] += 1; + } + } + } + + comm_flag = 1; + if (newton_pair) comm->reverse_comm(this, 2); + + for (i = 0; i < nlocal; i++) { + if (!(mask[i] & groupbit)) continue; + + if (norm[i] != 0) + array_atom[i][0] /= norm[i]; + else + array_atom[i][0] = 0.0; + array_atom[i][0] = sqrt(array_atom[i][0]); + + for (j = 0; j < 3; j++) + for (k = 0; k < 3; k++) + F_tmp[j][k] = F[i][j][k]; + + transpose_times3(F_tmp, F_tmp, E); + for (j = 0; j < dim; j++) E[j][j] -= 1.0; + + evol = (E[0][0] + E[1][1] + E[2][2]) / dim; + + // Calculate deviatoric strain + for (j = 0; j < dim; j++) E[j][j] -= evol; + j2 = 0.0; + for (j = 0; j < 3; j++) + for (k = 0; k < 3; k++) + j2 += E[j][k] * E[j][k]; + + edev = sqrt(0.5 * j2); + + array_atom[i][1] = evol; + array_atom[i][2] = edev; + } +} + +/* ---------------------------------------------------------------------- */ + +int FixNonaffineDisplacement::pack_reverse_comm(int n, int first, double *buf) +{ + int i, m, last, k, l; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + if (comm_flag == 0) { + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l++) { + buf[m++] = X[i][k][l]; + buf[m++] = Y[i][k][l]; + } + } + } else { + buf[m++] = array_atom[i][0]; + buf[m++] = ubuf(norm[i]).d; + } + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i, j, m, k, l; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + if (comm_flag == 0) { + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l++) { + X[j][k][l] += buf[m++]; + Y[j][k][l] += buf[m++]; + } + } + } else { + array_atom[j][0] += buf[m++]; + norm[j] += (int) ubuf(buf[m++]).i; + } + } +} + +/* ---------------------------------------------------------------------- */ + +int FixNonaffineDisplacement::pack_forward_comm(int n, int *list, double *buf, + int /*pbc_flag*/, int * /*pbc*/) +{ + int i, j, m, k, l; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l ++) { + buf[m++] = F[j][k][l]; + } + } + } + + return m; +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::unpack_forward_comm(int n, int first, double *buf) +{ + int i, m, last, k, l; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + for (k = 0; k < 3; k++) { + for (l = 0; l < 3; l ++) { + F[i][k][l] = buf[m++]; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::minimum_image0(double *delta) +{ + if (domain->triclinic == 0) { + if (domain->xperiodic) { + while (fabs(delta[0]) > xprd0_half) { + if (delta[0] < 0.0) delta[0] += xprd0; + else delta[0] -= xprd0; + } + } + if (domain->yperiodic) { + while (fabs(delta[1]) > yprd0_half) { + if (delta[1] < 0.0) delta[1] += yprd0; + else delta[1] -= yprd0; + } + } + if (domain->zperiodic) { + while (fabs(delta[2]) > zprd0_half) { + if (delta[2] < 0.0) delta[2] += zprd0; + else delta[2] -= zprd0; + } + } + + } else { + if (domain->zperiodic) { + while (fabs(delta[2]) > zprd0_half) { + if (delta[2] < 0.0) { + delta[2] += zprd0; + delta[1] += yz0; + delta[0] += xz0; + } else { + delta[2] -= zprd0; + delta[1] -= yz0; + delta[0] -= xz0; + } + } + } + if (domain->yperiodic) { + while (fabs(delta[1]) > yprd0_half) { + if (delta[1] < 0.0) { + delta[1] += yprd0; + delta[0] += xy0; + } else { + delta[1] -= yprd0; + delta[0] -= xy0; + } + } + } + if (domain->xperiodic) { + while (fabs(delta[0]) > xprd0_half) { + if (delta[0] < 0.0) delta[0] += xprd0; + else delta[0] -= xprd0; + } + } + } +} + +/* ---------------------------------------------------------------------- */ + +void FixNonaffineDisplacement::grow_arrays(int nmax_new) +{ + nmax = nmax_new; + memory->destroy(X); + memory->destroy(Y); + memory->destroy(F); + memory->destroy(norm); + memory->create(X, nmax, 3, 3, "fix_nonaffine_displacement:X"); + memory->create(Y, nmax, 3, 3, "fix_nonaffine_displacement:Y"); + memory->create(F, nmax, 3, 3, "fix_nonaffine_displacement:F"); + memory->create(norm, nmax, "fix_nonaffine_displacement:norm"); +} diff --git a/src/EXTRA-FIX/fix_nonaffine_displacement.h b/src/EXTRA-FIX/fix_nonaffine_displacement.h new file mode 100644 index 0000000000..0a195dc08e --- /dev/null +++ b/src/EXTRA-FIX/fix_nonaffine_displacement.h @@ -0,0 +1,71 @@ +/* -*- 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 FIX_CLASS +// clang-format off +FixStyle(nonaffine/displacement,FixNonaffineDisplacement) +// clang-format on +#else + +#ifndef LMP_FIX_NONAFFINE_DISPLACEMENT_H +#define LMP_FIX_NONAFFINE_DISPLACEMENT_H + +#include "fix.h" + +namespace LAMMPS_NS { + +class FixNonaffineDisplacement : public Fix { + public: + FixNonaffineDisplacement(class LAMMPS *, int, char **); + ~FixNonaffineDisplacement() override; + int setmask() override; + void post_constructor() override; + void init() override; + void init_list(int, class NeighList *) override; + void setup(int); + void post_force(int) override; + void write_restart(FILE *fp) override; + void restart(char *buf) override; + int pack_forward_comm(int, int *, double *, int, int *) override; + void unpack_forward_comm(int, int, double *) override; + int pack_reverse_comm(int, int, double *) override; + void unpack_reverse_comm(int, int *, double *) override; + + private: + double dtv; + char *id_fix; + class FixStoreAtom *fix; + int nmax, comm_flag; + int nad_style, cut_style; + int reference_style, offset_timestep, reference_timestep, update_timestep; + int reference_saved; + double cutoff_custom, cutsq_custom, mycutneigh; + double xprd0, yprd0, zprd0, xprd0_half, yprd0_half, zprd0_half, xy0, xz0, yz0; + + double ***X, ***Y, ***F; + int *norm; + + class NeighList *list; // half neighbor list + + + void integrate_velocity(); + void calculate_D2Min(); + void save_reference_state(); + void minimum_image0(double *); + void grow_arrays(int); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/GRANULAR/compute_fabric.cpp b/src/GRANULAR/compute_fabric.cpp index fb95a8b446..adaf242c92 100644 --- a/src/GRANULAR/compute_fabric.cpp +++ b/src/GRANULAR/compute_fabric.cpp @@ -11,6 +11,10 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + Contributing authors: Joel Clemmer (SNL), Ishan Srivastava (LBNL) +------------------------------------------------------------------------- */ + #include "compute_fabric.h" #include "atom.h" diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp index 61722bf62d..0980ad776d 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.cpp @@ -123,6 +123,10 @@ void PairPACEExtrapolationKokkos::grow(int natom, int maxneigh) // hard-core repulsion MemKK::realloc_kokkos(rho_core, "pace:rho_core", natom); MemKK::realloc_kokkos(dF_drho_core, "pace:dF_drho_core", natom); + MemKK::realloc_kokkos(dF_dfcut, "pace:dF_dfcut", natom); + MemKK::realloc_kokkos(d_d_min, "pace:r_min_pair", natom); + MemKK::realloc_kokkos(d_jj_min, "pace:j_min_pair", natom); + MemKK::realloc_kokkos(d_corerep, "pace:corerep", natom); // per-atom corerep MemKK::realloc_kokkos(dB_flatten, "pace:dB_flatten", natom, idx_ms_combs_max, basis_set->rankmax); @@ -219,6 +223,24 @@ void PairPACEExtrapolationKokkos::copy_pertype() Kokkos::deep_copy(d_wpre, h_wpre); Kokkos::deep_copy(d_mexp, h_mexp); + + + // ZBL core-rep + MemKK::realloc_kokkos(d_cut_in, "pace:d_cut_in", nelements, nelements); + MemKK::realloc_kokkos(d_dcut_in, "pace:d_dcut_in", nelements, nelements); + auto h_cut_in = Kokkos::create_mirror_view(d_cut_in); + auto h_dcut_in = Kokkos::create_mirror_view(d_dcut_in); + + for (int mu_i = 0; mu_i < nelements; ++mu_i) { + for (int mu_j = 0; mu_j < nelements; ++mu_j) { + h_cut_in(mu_i,mu_j) = basis_set->map_bond_specifications.at({mu_i,mu_j}).rcut_in; + h_dcut_in(mu_i,mu_j) = basis_set->map_bond_specifications.at({mu_i,mu_j}).dcut_in; + } + } + Kokkos::deep_copy(d_cut_in, h_cut_in); + Kokkos::deep_copy(d_dcut_in, h_dcut_in); + + is_zbl = basis_set->radial_functions->inner_cutoff_type == "zbl"; } /* ---------------------------------------------------------------------- */ @@ -573,13 +595,20 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in d_vatom = k_vatom.view(); } - if (gamma_flag && atom->nlocal > nmax) { + if (flag_compute_extrapolation_grade && atom->nlocal > nmax) { memory->destroy(extrapolation_grade_gamma); nmax = atom->nlocal; memory->create(extrapolation_grade_gamma, nmax, "pace/atom:gamma"); //zeroify array memset(extrapolation_grade_gamma, 0, nmax * sizeof(*extrapolation_grade_gamma)); } + if (flag_corerep_factor && atom->nlocal > nmax_corerep) { + memory->destroy(corerep_factor); + nmax_corerep = atom->nlocal; + memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep"); + //zeroify array + memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor)); + } copymode = 1; if (!force->newton_pair) @@ -631,8 +660,13 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in Kokkos::deep_copy(A_rank1, 0.0); Kokkos::deep_copy(rhos, 0.0); + Kokkos::deep_copy(rho_core, 0.0); + Kokkos::deep_copy(d_d_min, PairPACEExtrapolation::aceimpl->basis_set->cutoffmax); + Kokkos::deep_copy(d_jj_min, -1); + Kokkos::deep_copy(projections, 0.0); Kokkos::deep_copy(d_gamma, 0.0); + Kokkos::deep_copy(d_corerep, 0.0); EV_FLOAT ev_tmp; @@ -696,7 +730,7 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in } //ComputeGamma - if (gamma_flag) { + if (flag_compute_extrapolation_grade) { typename Kokkos::RangePolicy policy_gamma(0,chunk_size); Kokkos::parallel_for("ComputeGamma",policy_gamma,*this); } @@ -738,12 +772,17 @@ void PairPACEExtrapolationKokkos::compute(int eflag_in, int vflag_in } ev += ev_tmp; - //if gamma_flag - copy current d_gamma to extrapolation_grade_gamma - if (gamma_flag){ + //if flag_compute_extrapolation_grade - copy current d_gamma to extrapolation_grade_gamma + if (flag_compute_extrapolation_grade){ h_gamma = Kokkos::create_mirror_view(d_gamma); Kokkos::deep_copy(h_gamma, d_gamma); memcpy(extrapolation_grade_gamma+chunk_offset, (void *) h_gamma.data(), sizeof(double)*chunk_size); } + if (flag_corerep_factor) { + h_corerep = Kokkos::create_mirror_view(d_corerep); + Kokkos::deep_copy(h_corerep,d_corerep); + memcpy(corerep_factor+chunk_offset, (void *) h_corerep.data(), sizeof(double)*chunk_size); + } chunk_offset += chunk_size; } // end while @@ -799,6 +838,7 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeNeig const X_FLOAT ytmp = x(i,1); const X_FLOAT ztmp = x(i,2); const int jnum = d_numneigh[i]; + const int mu_i = d_map(type(i)); // get a pointer to scratch memory // This is used to cache whether or not an atom is within the cutoff @@ -858,6 +898,36 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeNeig } offset++; }); + + if (is_zbl) { + //adapted from https://www.osti.gov/servlets/purl/1429450 + if (ncount > 0) { + using minloc_value_type=Kokkos::MinLoc::value_type; + minloc_value_type djjmin; + djjmin.val=1e20; + djjmin.loc=-1; + Kokkos::MinLoc reducer_scalar(djjmin); + // loop over ncount (actual neighbours withing cutoff) rather than jnum (total number of neigh in cutoff+skin) + Kokkos::parallel_reduce(Kokkos::TeamThreadRange(team, ncount), + [&](const int offset, minloc_value_type &min_d_dist) { + int j = d_nearest(ii,offset); + j &= NEIGHMASK; + const int jtype = type(j); + auto r = d_rnorms(ii,offset); + const int mu_j = d_map(type(j)); + const F_FLOAT d = r - (d_cut_in(mu_i, mu_j) - d_dcut_in(mu_i, mu_j)); + if (d < min_d_dist.val) { + min_d_dist.val = d; + min_d_dist.loc = offset; + } + }, reducer_scalar); + d_d_min(ii) = djjmin.val; + d_jj_min(ii) = djjmin.loc;// d_jj_min should be NOT in 0..jnum range, but in 0..d_ncount(<=jnum) + } else { + d_d_min(ii) = 1e20; + d_jj_min(ii) = -1; + } + } } /* ---------------------------------------------------------------------- */ @@ -998,7 +1068,7 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeRho, //gamma_i - if (gamma_flag) + if (flag_compute_extrapolation_grade) Kokkos::atomic_add(&projections(ii, func_ind), d_gen_cgs(mu_i, idx_ms_comb) * A_cur); } else { // rank > 1 @@ -1037,7 +1107,7 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeRho, Kokkos::atomic_add(&rhos(ii, p), B.real_part_product(d_coeffs(mu_i, func_ind, p) * d_gen_cgs(mu_i, idx_ms_comb))); } //gamma_i - if (gamma_flag) + if (flag_compute_extrapolation_grade) Kokkos::atomic_add(&projections(ii, func_ind), B.real_part_product(d_gen_cgs(mu_i, idx_ms_comb))); } } @@ -1056,23 +1126,45 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeFS, const int ndensity = d_ndensity(mu_i); double evdwl, fcut, dfcut; + double evdwl_cut; evdwl = fcut = dfcut = 0.0; inner_cutoff(rho_core(ii), rho_cut, drho_cut, fcut, dfcut); FS_values_and_derivatives(ii, evdwl, mu_i); - dF_drho_core(ii) = evdwl * dfcut + 1; + if (is_zbl) { + if (d_jj_min(ii) != -1) { + const int mu_jmin = d_mu(ii,d_jj_min(ii)); + F_FLOAT dcutin = d_dcut_in(mu_i, mu_jmin); + F_FLOAT transition_coordinate = dcutin - d_d_min(ii); // == cutin - r_min + cutoff_func_poly(transition_coordinate, dcutin, dcutin, fcut, dfcut); + dfcut = -dfcut; // invert, because rho_core = cutin - r_min + } else { + // no neighbours + fcut = 1; + dfcut = 0; + } + evdwl_cut = evdwl * fcut + rho_core(ii) * (1 - fcut); // evdwl * fcut + rho_core_uncut - rho_core_uncut* fcut + dF_drho_core(ii) = 1 - fcut; + dF_dfcut(ii) = evdwl * dfcut - rho_core(ii) * dfcut; + } else { + inner_cutoff(rho_core(ii), rho_cut, drho_cut, fcut, dfcut); + dF_drho_core(ii) = evdwl * dfcut + 1; + evdwl_cut = evdwl * fcut + rho_core(ii); + } for (int p = 0; p < ndensity; ++p) dF_drho(ii, p) *= fcut; // tally energy contribution if (eflag) { - double evdwl_cut = evdwl * fcut + rho_core(ii); // E0 shift evdwl_cut += d_E0vals(mu_i); e_atom(ii) = evdwl_cut; } + + if (flag_corerep_factor) + d_corerep(ii) = 1-fcut; } /* ---------------------------------------------------------------------- */ @@ -1240,6 +1332,15 @@ void PairPACEExtrapolationKokkos::operator() (TagPairPACEComputeDeri f_ij(ii, jj, 0) = scale * f_ji[0] + fpair * r_hat[0]; f_ij(ii, jj, 1) = scale * f_ji[1] + fpair * r_hat[1]; f_ij(ii, jj, 2) = scale * f_ji[2] + fpair * r_hat[2]; + + if (is_zbl) { + if (jj==d_jj_min(ii)) { + // DCRU = 1.0 + f_ij(ii, jj, 0) += dF_dfcut(ii) * r_hat[0]; + f_ij(ii, jj, 1) += dF_dfcut(ii) * r_hat[1]; + f_ij(ii, jj, 2) += dF_dfcut(ii) * r_hat[2]; + } + } } /* ---------------------------------------------------------------------- */ @@ -1777,6 +1878,8 @@ double PairPACEExtrapolationKokkos::memory_usage() bytes += MemKK::memory_usage(weights_rank1); bytes += MemKK::memory_usage(rho_core); bytes += MemKK::memory_usage(dF_drho_core); + bytes += MemKK::memory_usage(dF_dfcut); + bytes += MemKK::memory_usage(d_corerep); bytes += MemKK::memory_usage(dB_flatten); bytes += MemKK::memory_usage(fr); bytes += MemKK::memory_usage(dfr); @@ -1794,6 +1897,8 @@ double PairPACEExtrapolationKokkos::memory_usage() bytes += MemKK::memory_usage(d_mu); bytes += MemKK::memory_usage(d_rhats); bytes += MemKK::memory_usage(d_rnorms); + bytes += MemKK::memory_usage(d_d_min); + bytes += MemKK::memory_usage(d_jj_min); bytes += MemKK::memory_usage(d_nearest); bytes += MemKK::memory_usage(f_ij); bytes += MemKK::memory_usage(d_rho_core_cutoff); @@ -1842,9 +1947,10 @@ double PairPACEExtrapolationKokkos::memory_usage() template void *PairPACEExtrapolationKokkos::extract(const char *str, int &dim) { - //check if str=="gamma_flag" then compute extrapolation grades on this iteration dim = 0; - if (strcmp(str, "gamma_flag") == 0) return (void *) &gamma_flag; + //check if str=="flag_compute_extrapolation_grade" then compute extrapolation grades on this iteration + if (strcmp(str, "gamma_flag") == 0) return (void *) &flag_compute_extrapolation_grade; + if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor; dim = 2; if (strcmp(str, "scale") == 0) return (void *) scale; @@ -1867,6 +1973,10 @@ void *PairPACEExtrapolationKokkos::extract_peratom(const char *str, ncol = 0; return (void *) extrapolation_grade_gamma; } + if (strcmp(str, "corerep") == 0) { + ncol = 0; + return (void *) corerep_factor; + } return nullptr; } diff --git a/src/KOKKOS/pair_pace_extrapolation_kokkos.h b/src/KOKKOS/pair_pace_extrapolation_kokkos.h index 55bcf4fead..aa6c49c36d 100644 --- a/src/KOKKOS/pair_pace_extrapolation_kokkos.h +++ b/src/KOKKOS/pair_pace_extrapolation_kokkos.h @@ -106,7 +106,6 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { protected: int inum, maxneigh, chunk_size, chunk_offset, idx_ms_combs_max, total_num_functions_max; int host_flag; - int gamma_flag; int eflag, vflag; @@ -130,6 +129,7 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { tdual_fparams k_cutsq, k_scale; typedef Kokkos::View t_fparams; t_fparams d_cutsq, d_scale; + t_fparams d_cut_in, d_dcut_in; // inner cutoff typename AT::t_int_1d d_map; @@ -234,12 +234,16 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { t_ace_2d rhos; t_ace_2d dF_drho; + t_ace_3c dB_flatten; + // hard-core repulsion t_ace_1d rho_core; - t_ace_3c dB_flatten; t_ace_2d cr; t_ace_2d dcr; t_ace_1d dF_drho_core; + t_ace_1d dF_dfcut; + t_ace_1d d_corerep; + th_ace_1d h_corerep; // radial functions t_ace_4d fr; @@ -282,6 +286,11 @@ class PairPACEExtrapolationKokkos : public PairPACEExtrapolation { t_ace_3d3 d_rhats; t_ace_2i d_nearest; + // for ZBL core-rep implementation + t_ace_1d d_d_min; // [i] -> min-d for atom ii, d=d = r - (cut_in(mu_i, mu_j) - dcut_in(mu_i, mu_j)) + t_ace_1i d_jj_min; // [i] -> jj-index of nearest neigh (by r-(cut_in-dcut_in) criterion) + bool is_zbl; + // per-type t_ace_1i d_ndensity; t_ace_1i d_npoti; diff --git a/src/KOKKOS/pair_pace_kokkos.cpp b/src/KOKKOS/pair_pace_kokkos.cpp index 153a6d0333..805d7f68bb 100644 --- a/src/KOKKOS/pair_pace_kokkos.cpp +++ b/src/KOKKOS/pair_pace_kokkos.cpp @@ -121,6 +121,11 @@ void PairPACEKokkos::grow(int natom, int maxneigh) // hard-core repulsion MemKK::realloc_kokkos(rho_core, "pace:rho_core", natom); MemKK::realloc_kokkos(dF_drho_core, "pace:dF_drho_core", natom); + MemKK::realloc_kokkos(dF_dfcut, "pace:dF_dfcut", natom); + MemKK::realloc_kokkos(d_d_min, "pace:r_min_pair", natom); + MemKK::realloc_kokkos(d_jj_min, "pace:j_min_pair", natom); + MemKK::realloc_kokkos(d_corerep, "pace:corerep", natom); // per-atom corerep + MemKK::realloc_kokkos(dB_flatten, "pace:dB_flatten", natom, idx_rho_max, basis_set->rankmax); } @@ -212,6 +217,23 @@ void PairPACEKokkos::copy_pertype() Kokkos::deep_copy(d_wpre, h_wpre); Kokkos::deep_copy(d_mexp, h_mexp); + + // ZBL core-rep + MemKK::realloc_kokkos(d_cut_in, "pace:d_cut_in", nelements, nelements); + MemKK::realloc_kokkos(d_dcut_in, "pace:d_dcut_in", nelements, nelements); + auto h_cut_in = Kokkos::create_mirror_view(d_cut_in); + auto h_dcut_in = Kokkos::create_mirror_view(d_dcut_in); + + for (int mu_i = 0; mu_i < nelements; ++mu_i) { + for (int mu_j = 0; mu_j < nelements; ++mu_j) { + h_cut_in(mu_i,mu_j) = basis_set->map_bond_specifications.at({mu_i,mu_j}).rcut_in; + h_dcut_in(mu_i,mu_j) = basis_set->map_bond_specifications.at({mu_i,mu_j}).dcut_in; + } + } + Kokkos::deep_copy(d_cut_in, h_cut_in); + Kokkos::deep_copy(d_dcut_in, h_dcut_in); + + is_zbl = basis_set->radial_functions->inner_cutoff_type == "zbl"; } /* ---------------------------------------------------------------------- */ @@ -535,6 +557,13 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) memoryKK->create_kokkos(k_vatom,vatom,maxvatom,"pair:vatom"); d_vatom = k_vatom.view(); } + if (flag_corerep_factor && atom->nlocal > nmax_corerep) { + memory->destroy(corerep_factor); + nmax_corerep = atom->nlocal; + memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep"); + //zeroify array + memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor)); + } copymode = 1; if (!force->newton_pair) @@ -588,6 +617,9 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) Kokkos::deep_copy(A_rank1, 0.0); Kokkos::deep_copy(rhos, 0.0); Kokkos::deep_copy(rho_core, 0.0); + Kokkos::deep_copy(d_d_min, PairPACE::aceimpl->basis_set->cutoffmax); + Kokkos::deep_copy(d_jj_min, -1); + Kokkos::deep_copy(d_corerep, 0.0); EV_FLOAT ev_tmp; @@ -686,6 +718,13 @@ void PairPACEKokkos::compute(int eflag_in, int vflag_in) } } ev += ev_tmp; + + if (flag_corerep_factor) { + h_corerep = Kokkos::create_mirror_view(d_corerep); + Kokkos::deep_copy(h_corerep,d_corerep); + memcpy(corerep_factor+chunk_offset, (void *) h_corerep.data(), sizeof(double)*chunk_size); + } + chunk_offset += chunk_size; } // end while @@ -741,6 +780,7 @@ void PairPACEKokkos::operator() (TagPairPACEComputeNeigh,const typen const X_FLOAT ytmp = x(i,1); const X_FLOAT ztmp = x(i,2); const int jnum = d_numneigh[i]; + const int mu_i = d_map(type(i)); // get a pointer to scratch memory // This is used to cache whether or not an atom is within the cutoff @@ -800,6 +840,36 @@ void PairPACEKokkos::operator() (TagPairPACEComputeNeigh,const typen } offset++; }); + + if (is_zbl) { + //adapted from https://www.osti.gov/servlets/purl/1429450 + if(ncount>0) { + using minloc_value_type=Kokkos::MinLoc::value_type; + minloc_value_type djjmin; + djjmin.val=1e20; + djjmin.loc=-1; + Kokkos::MinLoc reducer_scalar(djjmin); + // loop over ncount (actual neighbours withing cutoff) rather than jnum (total number of neigh in cutoff+skin) + Kokkos::parallel_reduce(Kokkos::TeamThreadRange(team, ncount), + [&](const int offset, minloc_value_type &min_d_dist) { + int j = d_nearest(ii,offset); + j &= NEIGHMASK; + const int jtype = type(j); + auto r = d_rnorms(ii,offset); + const int mu_j = d_map(type(j)); + const F_FLOAT d = r - (d_cut_in(mu_i, mu_j) - d_dcut_in(mu_i, mu_j)); + if (d < min_d_dist.val) { + min_d_dist.val = d; + min_d_dist.loc = offset; + } + }, reducer_scalar); + d_d_min(ii) = djjmin.val; + d_jj_min(ii) = djjmin.loc;// d_jj_min should be NOT in 0..jnum range, but in 0..d_ncount(<=jnum) + } else { + d_d_min(ii) = 1e20; + d_jj_min(ii) = -1; + } + } } /* ---------------------------------------------------------------------- */ @@ -990,23 +1060,42 @@ void PairPACEKokkos::operator() (TagPairPACEComputeFS, const int& ii const int ndensity = d_ndensity(mu_i); double evdwl, fcut, dfcut; + double evdwl_cut; evdwl = fcut = dfcut = 0.0; - inner_cutoff(rho_core(ii), rho_cut, drho_cut, fcut, dfcut); FS_values_and_derivatives(ii, evdwl, mu_i); - - dF_drho_core(ii) = evdwl * dfcut + 1; + if (is_zbl) { + if (d_jj_min(ii) != -1) { + const int mu_jmin = d_mu(ii,d_jj_min(ii)); + F_FLOAT dcutin = d_dcut_in(mu_i, mu_jmin); + F_FLOAT transition_coordinate = dcutin - d_d_min(ii); // == cutin - r_min + cutoff_func_poly(transition_coordinate, dcutin, dcutin, fcut, dfcut); + dfcut = -dfcut; // invert, because rho_core = cutin - r_min + } else { + // no neighbours + fcut = 1; + dfcut = 0; + } + evdwl_cut = evdwl * fcut + rho_core(ii) * (1 - fcut); // evdwl * fcut + rho_core_uncut - rho_core_uncut* fcut + dF_drho_core(ii) = 1 - fcut; + dF_dfcut(ii) = evdwl * dfcut - rho_core(ii) * dfcut; + } else { + inner_cutoff(rho_core(ii), rho_cut, drho_cut, fcut, dfcut); + dF_drho_core(ii) = evdwl * dfcut + 1; + evdwl_cut = evdwl * fcut + rho_core(ii); + } for (int p = 0; p < ndensity; ++p) - dF_drho(ii, p) *= fcut; - + dF_drho(ii, p) *= fcut; // tally energy contribution if (eflag) { - double evdwl_cut = evdwl * fcut + rho_core(ii); - // E0 shift - evdwl_cut += d_E0vals(mu_i); - e_atom(ii) = evdwl_cut; + // E0 shift + evdwl_cut += d_E0vals(mu_i); + e_atom(ii) = evdwl_cut; } + + if (flag_corerep_factor) + d_corerep(ii) = 1-fcut; } /* ---------------------------------------------------------------------- */ @@ -1146,6 +1235,15 @@ void PairPACEKokkos::operator() (TagPairPACEComputeDerivative, const f_ij(ii, jj, 0) = scale * f_ji[0] + fpair * r_hat[0]; f_ij(ii, jj, 1) = scale * f_ji[1] + fpair * r_hat[1]; f_ij(ii, jj, 2) = scale * f_ji[2] + fpair * r_hat[2]; + + if (is_zbl) { + if (jj==d_jj_min(ii)) { + // DCRU = 1.0 + f_ij(ii, jj, 0) += dF_dfcut(ii) * r_hat[0]; + f_ij(ii, jj, 1) += dF_dfcut(ii) * r_hat[1]; + f_ij(ii, jj, 2) += dF_dfcut(ii) * r_hat[2]; + } + } } /* ---------------------------------------------------------------------- */ @@ -1683,6 +1781,8 @@ double PairPACEKokkos::memory_usage() bytes += MemKK::memory_usage(weights_rank1); bytes += MemKK::memory_usage(rho_core); bytes += MemKK::memory_usage(dF_drho_core); + bytes += MemKK::memory_usage(dF_dfcut); + bytes += MemKK::memory_usage(d_corerep); bytes += MemKK::memory_usage(dB_flatten); bytes += MemKK::memory_usage(fr); bytes += MemKK::memory_usage(dfr); @@ -1700,6 +1800,8 @@ double PairPACEKokkos::memory_usage() bytes += MemKK::memory_usage(d_mu); bytes += MemKK::memory_usage(d_rhats); bytes += MemKK::memory_usage(d_rnorms); + bytes += MemKK::memory_usage(d_d_min); + bytes += MemKK::memory_usage(d_jj_min); bytes += MemKK::memory_usage(d_nearest); bytes += MemKK::memory_usage(f_ij); bytes += MemKK::memory_usage(d_rho_core_cutoff); diff --git a/src/KOKKOS/pair_pace_kokkos.h b/src/KOKKOS/pair_pace_kokkos.h index 39cfd100f8..36486f8628 100644 --- a/src/KOKKOS/pair_pace_kokkos.h +++ b/src/KOKKOS/pair_pace_kokkos.h @@ -121,6 +121,7 @@ class PairPACEKokkos : public PairPACE { tdual_fparams k_cutsq, k_scale; typedef Kokkos::View t_fparams; t_fparams d_cutsq, d_scale; + t_fparams d_cut_in, d_dcut_in; // inner cutoff typename AT::t_int_1d d_map; @@ -209,6 +210,8 @@ class PairPACEKokkos : public PairPACE { typedef Kokkos::View t_ace_4c; typedef Kokkos::View t_ace_4c3; + typedef typename Kokkos::View::HostMirror th_ace_1d; + t_ace_3d A_rank1; t_ace_4c A; @@ -222,12 +225,16 @@ class PairPACEKokkos : public PairPACE { t_ace_2d rhos; t_ace_2d dF_drho; + t_ace_3c dB_flatten; + // hard-core repulsion t_ace_1d rho_core; - t_ace_3c dB_flatten; t_ace_2d cr; t_ace_2d dcr; t_ace_1d dF_drho_core; + t_ace_1d dF_dfcut; + t_ace_1d d_corerep; + th_ace_1d h_corerep; // radial functions t_ace_4d fr; @@ -265,6 +272,11 @@ class PairPACEKokkos : public PairPACE { t_ace_3d3 d_rhats; t_ace_2i d_nearest; + // for ZBL core-rep implementation + t_ace_1d d_d_min; // [i] -> min-d for atom ii, d=d = r - (cut_in(mu_i, mu_j) - dcut_in(mu_i, mu_j)) + t_ace_1i d_jj_min; // [i] -> jj-index of nearest neigh (by r-(cut_in-dcut_in) criterion) + bool is_zbl; + // per-type t_ace_1i d_ndensity; t_ace_1i d_npoti; diff --git a/src/ML-PACE/pair_pace.cpp b/src/ML-PACE/pair_pace.cpp index 57f12597d1..e9bd25f9d7 100644 --- a/src/ML-PACE/pair_pace.cpp +++ b/src/ML-PACE/pair_pace.cpp @@ -45,6 +45,7 @@ Copyright 2021 Yury Lysogorskiy^1, Cas van der Oord^2, Anton Bochkarev^1, #include "ace-evaluator/ace_evaluator.h" #include "ace-evaluator/ace_recursive.h" #include "ace-evaluator/ace_version.h" +#include "ace/ace_b_basis.h" namespace LAMMPS_NS { struct ACEImpl { @@ -87,6 +88,10 @@ PairPACE::PairPACE(LAMMPS *lmp) : Pair(lmp) one_coeff = 1; manybody_flag = 1; + nmax_corerep = 0; + flag_corerep_factor = 0; + corerep_factor = nullptr; + aceimpl = new ACEImpl; recursive = false; @@ -109,6 +114,7 @@ PairPACE::~PairPACE() memory->destroy(setflag); memory->destroy(cutsq); memory->destroy(scale); + memory->destroy(corerep_factor); } } @@ -143,10 +149,18 @@ void PairPACE::compute(int eflag, int vflag) // the pointer to the list of neighbors of "i" firstneigh = list->firstneigh; + if (flag_corerep_factor && atom->nlocal > nmax_corerep) { + memory->destroy(corerep_factor); + nmax_corerep = atom->nlocal; + memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep_factor"); + //zeroify array + memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor)); + } + //determine the maximum number of neighbours int max_jnum = 0; int nei = 0; - for (ii = 0; ii < list->inum; ii++) { + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jnum = numneigh[i]; nei = nei + jnum; @@ -156,7 +170,7 @@ void PairPACE::compute(int eflag, int vflag) aceimpl->ace->resize_neighbours_cache(max_jnum); //loop over atoms - for (ii = 0; ii < list->inum; ii++) { + for (ii = 0; ii < inum; ii++) { i = list->ilist[ii]; const int itype = type[i]; @@ -181,6 +195,9 @@ void PairPACE::compute(int eflag, int vflag) error->one(FLERR, e.what()); } + if (flag_corerep_factor) + corerep_factor[i] = 1 - aceimpl->ace->ace_fcut; + // 'compute_atom' will update the `aceimpl->ace->e_atom` and `aceimpl->ace->neighbours_forces(jj, alpha)` arrays for (jj = 0; jj < jnum; jj++) { @@ -287,7 +304,14 @@ void PairPACE::coeff(int narg, char **arg) //load potential file delete aceimpl->basis_set; if (comm->me == 0) utils::logmesg(lmp, "Loading {}\n", potential_file_name); - aceimpl->basis_set = new ACECTildeBasisSet(potential_file_name); + // if potential is in ACEBBasisSet (YAML) format, then convert to ACECTildeBasisSet automatically + if (utils::strmatch(potential_file_name,".*\\.yaml$")) { + ACEBBasisSet bBasisSet = ACEBBasisSet(potential_file_name); + ACECTildeBasisSet cTildeBasisSet = bBasisSet.to_ACECTildeBasisSet(); + aceimpl->basis_set = new ACECTildeBasisSet(cTildeBasisSet); + } else { + aceimpl->basis_set = new ACECTildeBasisSet(potential_file_name); + } if (comm->me == 0) { utils::logmesg(lmp, "Total number of basis functions\n"); @@ -374,7 +398,29 @@ double PairPACE::init_one(int i, int j) ---------------------------------------------------------------------- */ void *PairPACE::extract(const char *str, int &dim) { + dim = 0; + //check if str=="corerep_flag" then compute extrapolation grades on this iteration + if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor; + dim = 2; if (strcmp(str, "scale") == 0) return (void *) scale; return nullptr; } + +/* ---------------------------------------------------------------------- + peratom requests from FixPair + return ptr to requested data + also return ncol = # of quantites per atom + 0 = per-atom vector + 1 or more = # of columns in per-atom array + return NULL if str is not recognized +---------------------------------------------------------------------- */ +void *PairPACE::extract_peratom(const char *str, int &ncol) +{ + if (strcmp(str, "corerep") == 0) { + ncol = 0; + return (void *) corerep_factor; + } + + return nullptr; +} diff --git a/src/ML-PACE/pair_pace.h b/src/ML-PACE/pair_pace.h index 94649ecaab..9b5d2c5480 100644 --- a/src/ML-PACE/pair_pace.h +++ b/src/ML-PACE/pair_pace.h @@ -48,11 +48,15 @@ class PairPACE : public Pair { double init_one(int, int) override; void *extract(const char *, int &) override; + void *extract_peratom(const char *, int &) override; protected: struct ACEImpl *aceimpl; + int nmax_corerep = 0; virtual void allocate(); + double *corerep_factor; //per-atom core-rep factor (= 1 - fcut) + int flag_corerep_factor; double **scale; bool recursive; // "recursive" option for ACERecursiveEvaluator diff --git a/src/ML-PACE/pair_pace_extrapolation.cpp b/src/ML-PACE/pair_pace_extrapolation.cpp index dc0fb1848b..d9b8d3588a 100644 --- a/src/ML-PACE/pair_pace_extrapolation.cpp +++ b/src/ML-PACE/pair_pace_extrapolation.cpp @@ -93,11 +93,14 @@ PairPACEExtrapolation::PairPACEExtrapolation(LAMMPS *lmp) : Pair(lmp) manybody_flag = 1; nmax = 0; + nmax_corerep = 0; aceimpl = new ACEALImpl; scale = nullptr; flag_compute_extrapolation_grade = 0; extrapolation_grade_gamma = nullptr; + flag_corerep_factor = 0; + corerep_factor = nullptr; chunksize = 4096; } @@ -118,6 +121,7 @@ PairPACEExtrapolation::~PairPACEExtrapolation() memory->destroy(scale); memory->destroy(map); memory->destroy(extrapolation_grade_gamma); + memory->destroy(corerep_factor); } } @@ -166,11 +170,18 @@ void PairPACEExtrapolation::compute(int eflag, int vflag) //zeroify array memset(extrapolation_grade_gamma, 0, nmax * sizeof(*extrapolation_grade_gamma)); } + if (flag_corerep_factor && atom->nlocal > nmax_corerep) { + memory->destroy(corerep_factor); + nmax_corerep = atom->nlocal; + memory->create(corerep_factor, nmax_corerep, "pace/atom:corerep_factor"); + //zeroify array + memset(corerep_factor, 0, nmax_corerep * sizeof(*corerep_factor)); + } //determine the maximum number of neighbours int max_jnum = 0; int nei = 0; - for (ii = 0; ii < list->inum; ii++) { + for (ii = 0; ii < inum; ii++) { i = ilist[ii]; jnum = numneigh[i]; nei = nei + jnum; @@ -183,7 +194,7 @@ void PairPACEExtrapolation::compute(int eflag, int vflag) aceimpl->rec_ace->resize_neighbours_cache(max_jnum); //loop over atoms - for (ii = 0; ii < list->inum; ii++) { + for (ii = 0; ii < inum; ii++) { i = list->ilist[ii]; const int itype = type[i]; @@ -216,6 +227,11 @@ void PairPACEExtrapolation::compute(int eflag, int vflag) if (flag_compute_extrapolation_grade) extrapolation_grade_gamma[i] = aceimpl->ace->max_gamma_grade; + if (flag_corerep_factor) { + corerep_factor[i] = 1 - (flag_compute_extrapolation_grade ? aceimpl->ace->ace_fcut + : aceimpl->rec_ace->ace_fcut); + } + Array2D &neighbours_forces = (flag_compute_extrapolation_grade ? aceimpl->ace->neighbours_forces : aceimpl->rec_ace->neighbours_forces); @@ -437,9 +453,11 @@ double PairPACEExtrapolation::init_one(int i, int j) ---------------------------------------------------------------------- */ void *PairPACEExtrapolation::extract(const char *str, int &dim) { - //check if str=="gamma_flag" then compute extrapolation grades on this iteration dim = 0; + //check if str=="gamma_flag" then compute extrapolation grades on this iteration if (strcmp(str, "gamma_flag") == 0) return (void *) &flag_compute_extrapolation_grade; + //check if str=="corerep_flag" then compute extrapolation grades on this iteration + if (strcmp(str, "corerep_flag") == 0) return (void *) &flag_corerep_factor; dim = 2; if (strcmp(str, "scale") == 0) return (void *) scale; @@ -461,5 +479,10 @@ void *PairPACEExtrapolation::extract_peratom(const char *str, int &ncol) return (void *) extrapolation_grade_gamma; } + if (strcmp(str, "corerep") == 0) { + ncol = 0; + return (void *) corerep_factor; + } + return nullptr; } diff --git a/src/ML-PACE/pair_pace_extrapolation.h b/src/ML-PACE/pair_pace_extrapolation.h index 6f7eeb279e..2dcec04d4b 100644 --- a/src/ML-PACE/pair_pace_extrapolation.h +++ b/src/ML-PACE/pair_pace_extrapolation.h @@ -47,13 +47,15 @@ class PairPACEExtrapolation : public Pair { protected: struct ACEALImpl *aceimpl; - int nmax; + int nmax = 0, nmax_corerep = 0; virtual void allocate(); std::vector element_names; // list of elements (used by dump pace/extrapolation) - double *extrapolation_grade_gamma; //per-atom gamma value + double *extrapolation_grade_gamma = nullptr; //per-atom gamma value + double *corerep_factor = nullptr; //per-atom core-rep factor (= 1 - fcut) - int flag_compute_extrapolation_grade; + int flag_compute_extrapolation_grade = 0; + int flag_corerep_factor = 0; double **scale;