lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'master' into prepare-clang-format 

# Conflicts:
#	src/MOLECULE/bond_fene.h
#	src/MOLECULE/bond_fene_expand.h
This commit is contained in:
Axel Kohlmeyer 2021-05-11 21:49:48 -04:00
parent d1941392d3 ee546e2fa3
commit d8291eea7b
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
205 changed files with 16072 additions and 8024 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitattributes vendored
View File

@ -1,3 +1,4 @@
.gitattributes export-ignore
.gitignore export-ignore
.github export-ignore
.lgtm.yml export-ignore

4
.github/codeql/cpp.yml vendored Normal file
View File

@ -0,0 +1,4 @@
paths:
- src
- lib
- tools

5
.github/codeql/python.yml vendored Normal file
View File

@ -0,0 +1,5 @@
paths:
- python/lammps
queries:
- uses: security-and-quality

8
.github/workflows/codeql-analysis.yml vendored
View File

@ -31,16 +31,18 @@ jobs:
uses: github/codeql-action/init@v1
with:
languages: ${{ matrix.language }}
config-file: ./.github/codeql/${{ matrix.language }}.yml
- name: Create Build Environment
run: cmake -E make_directory ${{github.workspace}}/build
if: ${{ matrix.language == 'cpp' }}
run: mkdir build
- name: Building LAMMPS via CMake
if: ${{ matrix.language == 'cpp' }}
shell: bash
working-directory: ${{github.workspace}}/build
working-directory: build
run: |
cmake -C $GITHUB_WORKSPACE/cmake/presets/most.cmake $GITHUB_WORKSPACE/cmake
cmake -C ../cmake/presets/most.cmake ../cmake
cmake --build . --parallel 2
- name: Perform CodeQL Analysis

32
.github/workflows/unittest-macos.yml vendored
View File

@ -10,6 +10,8 @@ jobs:
name: MacOS Unit Test
if: ${{ github.repository == 'lammps/lammps' }}
runs-on: macos-latest
env:
CCACHE_DIR: ${{ github.workspace }}/.ccache
steps:
- name: Checkout repository
@ -17,20 +19,36 @@ jobs:
with:
fetch-depth: 2
- name: Install ccache
run: brew install ccache
- name: Create Build Environment
run: cmake -E make_directory ${{github.workspace}}/build
run: mkdir build
- name: Set up ccache
uses: actions/cache@v2
with:
path: ${{ env.CCACHE_DIR }}
key: macos-ccache-${{ github.sha }}
restore-keys: macos-ccache-
- name: Building LAMMPS via CMake
shell: bash
working-directory: ${{github.workspace}}/build
working-directory: build
run: |
cmake -C $GITHUB_WORKSPACE/cmake/presets/clang.cmake \
-C $GITHUB_WORKSPACE/cmake/presets/most.cmake \
$GITHUB_WORKSPACE/cmake \
-DENABLE_TESTING=ON -DBUILD_SHARED_LIBS=ON -DLAMMPS_EXCEPTIONS=ON
ccache -z
cmake -C ../cmake/presets/clang.cmake \
-C ../cmake/presets/most.cmake \
-D CMAKE_CXX_COMPILER_LAUNCHER=ccache \
-D CMAKE_C_COMPILER_LAUNCHER=ccache \
-D ENABLE_TESTING=on \
-D BUILD_SHARED_LIBS=on \
-D LAMMPS_EXCEPTIONS=on \
../cmake
cmake --build . --parallel 2
ccache -s
- name: Run Tests
working-directory: ${{github.workspace}}/build
working-directory: build
shell: bash
run: ctest -V

14
.lgtm.yml Normal file
View File

@ -0,0 +1,14 @@
extraction:
cpp:
configure:
command:
- "mkdir build"
- "cd build"
- "cmake -G Ninja -C ../cmake/presets/minimal.cmake ../cmake"
index:
build_command:
- "cd build"
- "ninja"
python:
python_setup:
version: 3

3
cmake/Modules/GenerateOpenCLHeader.cmake Normal file
View File

@ -0,0 +1,3 @@
# utility script to call WriteOpenCLHeader function
include(${SOURCE_DIR}/Modules/OpenCLUtils.cmake)
WriteOpenCLHeader(${VARNAME} ${HEADER_FILE} ${SOURCE_FILES})

4
cmake/Modules/OpenCLLoader.cmake
View File

@ -1,6 +1,6 @@
message(STATUS "Downloading and building OpenCL loader library")
set(OPENCL_LOADER_URL "${LAMMPS_THIRDPARTY_URL}/opencl-loader-2020.12.18.tar.gz" CACHE STRING "URL for OpenCL loader tarball")
set(OPENCL_LOADER_MD5 "011cdcbd41030be94f3fced6d763a52a" CACHE STRING "MD5 checksum of OpenCL loader tarball")
set(OPENCL_LOADER_URL "${LAMMPS_THIRDPARTY_URL}/opencl-loader-2021.05.02.tar.gz" CACHE STRING "URL for OpenCL loader tarball")
set(OPENCL_LOADER_MD5 "29180b05056578afda92f0d394c3a373" CACHE STRING "MD5 checksum of OpenCL loader tarball")
mark_as_advanced(OPENCL_LOADER_URL)
mark_as_advanced(OPENCL_LOADER_MD5)

20
cmake/Modules/OpenCLUtils.cmake
View File

@ -1,10 +1,8 @@
function(GenerateOpenCLHeader varname outfile files)
message("Creating ${outfile}...")
function(WriteOpenCLHeader varname outfile files)
file(WRITE ${outfile} "const char * ${varname} = \n")
math(EXPR ARG_END "${ARGC}-1")
separate_arguments(files)
foreach(IDX RANGE 2 ${ARG_END})
list(GET ARGV ${IDX} filename)
foreach(filename ${files})
file(READ ${filename} content)
string(REGEX REPLACE "\\s*//[^\n]*\n" "\n" content "${content}")
string(REGEX REPLACE "\\\\" "\\\\\\\\" content "${content}")
@ -15,4 +13,16 @@ function(GenerateOpenCLHeader varname outfile files)
endforeach()
file(APPEND ${outfile} ";\n")
endfunction(WriteOpenCLHeader)
function(GenerateOpenCLHeader varname outfile files)
list(REMOVE_AT ARGV 0 1)
add_custom_command(OUTPUT ${outfile}
COMMAND ${CMAKE_COMMAND} -D SOURCE_DIR=${CMAKE_CURRENT_SOURCE_DIR}
-D VARNAME=${varname}
-D HEADER_FILE=${outfile}
-D SOURCE_FILES="${ARGV}"
-P ${CMAKE_CURRENT_SOURCE_DIR}/Modules/GenerateOpenCLHeader.cmake
DEPENDS ${ARGV}
COMMENT "Generating ${outfile}...")
endfunction(GenerateOpenCLHeader)

9
cmake/Modules/Packages/GPU.cmake
View File

@ -145,7 +145,14 @@ if(GPU_API STREQUAL "CUDA")
target_include_directories(nvc_get_devices PRIVATE ${CUDA_INCLUDE_DIRS})
elseif(GPU_API STREQUAL "OPENCL")
option(USE_STATIC_OPENCL_LOADER "Download and include a static OpenCL ICD loader" ON)
# the static OpenCL loader doesn't seem to work on macOS. use the system provided
# version by default instead (for as long as it will be available)
if("${CMAKE_SYSTEM_NAME}" STREQUAL "Darwin")
set(_opencl_static_default OFF)
else()
set(_opencl_static_default ON)
endif()
option(USE_STATIC_OPENCL_LOADER "Download and include a static OpenCL ICD loader" ${_opencl_static_default})
mark_as_advanced(USE_STATIC_OPENCL_LOADER)
if (USE_STATIC_OPENCL_LOADER)
include(OpenCLLoader)

12
doc/src/Library_config.rst
View File

@ -15,6 +15,8 @@ This section documents the following functions:
- :cpp:func:`lammps_config_package_count`
- :cpp:func:`lammps_config_package_name`
- :cpp:func:`lammps_config_accelerator`
- :cpp:func:`lammps_has_gpu_device`
- :cpp:func:`lammps_gpu_device_info`
- :cpp:func:`lammps_has_style`
- :cpp:func:`lammps_style_count`
- :cpp:func:`lammps_style_name`
@ -132,6 +134,16 @@ approach.
-----------------------
.. doxygenfunction:: lammps_has_gpu_device
:project: progguide
-----------------------
.. doxygenfunction:: lammps_get_gpu_device_info
:project: progguide
-----------------------
.. doxygenfunction:: lammps_has_style
:project: progguide

98
doc/src/Tools.rst
View File

@ -94,6 +94,7 @@ Miscellaneous tools
* :ref:`kate <kate>`
* :ref:`LAMMPS shell <lammps_shell>`
* :ref:`LAMMPS magic patterns for file(1) <magic>`
* :ref:`Offline build tool <offline>`
* :ref:`singularity <singularity_tool>`
* :ref:`SWIG interface <swig>`
* :ref:`vim <vim>`
@ -756,6 +757,103 @@ See the README file in the tools/msi2lmp folder for more information.
----------
.. _offline:
Scripts for building LAMMPS when offline
----------------------------------------
In some situations it might be necessary to build LAMMPS on a system
without direct internet access. The scripts in ``tools/offline`` folder
allow you to pre-load external dependencies for both the documentation
build and for building LAMMPS with CMake.
It does so by
#. downloading necessary ``pip`` packages,
#. cloning ``git`` repositories
#. downloading tarballs
to a designated cache folder.
As of April 2021, all of these downloads make up around 600MB. By
default, the offline scripts will download everything into the
``$HOME/.cache/lammps`` folder, but this can be changed by setting the
``LAMMPS_CACHING_DIR`` environment variable.
Once the caches have been initialized, they can be used for building the
LAMMPS documentation or compiling LAMMPS using CMake on an offline
system.
The ``use_caches.sh`` script must be sourced into the current shell
to initialize the offline build environment. Note that it must use
the same ``LAMMPS_CACHING_DIR``. This script does the following:
#. Set up environment variables that modify the behavior of both,
``pip`` and ``git``
#. Start a simple local HTTP server using Python to host files for CMake
Afterwards, it will print out instruction on how to modify the CMake
command line to make sure it uses the local HTTP server.
To undo the environment changes and shutdown the local HTTP server,
run the ``deactivate_caches`` command.
Examples
^^^^^^^^
For all of the examples below, you first need to create the cache, which
requires an internet connection.
.. code-block:: bash
./tools/offline/init_caches.sh
Afterwards, you can disconnect or copy the contents of the
``LAMMPS_CACHING_DIR`` folder to an offline system.
Documentation Build
^^^^^^^^^^^^^^^^^^^
The documentation build will create a new virtual environment that
typically first installs dependencies from ``pip``. With the offline
environment loaded, these installations will instead grab the necessary
packages from your local cache.
.. code-block:: bash
# if LAMMPS_CACHING_DIR is different from default, make sure to set it first
# export LAMMPS_CACHING_DIR=path/to/folder
source tools/offline/use_caches.sh
cd doc/
make html
deactivate_caches
CMake Build
^^^^^^^^^^^
When compiling certain packages with external dependencies, the CMake
build system will download necessary files or sources from the web. For
more flexibility the CMake configuration allows users to specify the URL
of each of these dependencies. What the ``init_caches.sh`` script does
is create a CMake "preset" file, which sets the URLs for all of the known
dependencies and redirects the download to the local cache.
.. code-block:: bash
# if LAMMPS_CACHING_DIR is different from default, make sure to set it first
# export LAMMPS_CACHING_DIR=path/to/folder
source tools/offline/use_caches.sh
mkdir build
cd build
cmake -D LAMMPS_DOWNLOADS_URL=${HTTP_CACHE_URL} -C "${LAMMPS_HTTP_CACHE_CONFIG}" -C ../cmake/presets/most.cmake ../cmake
make -j 8
deactivate_caches
----------
.. _phonon:
phonon tool

12
doc/src/pair_hybrid.rst
View File

@ -257,7 +257,7 @@ then overrides them with 0.0 only for CHARMM:
.. code-block:: LAMMPS
special_bonds amber
pair_hybrid lj/charmm/coul/long 8.0 10.0 lj/cut/coul/long 10.0
pair_style hybrid lj/charmm/coul/long 8.0 10.0 lj/cut/coul/long 10.0
pair_modify pair lj/charmm/coul/long special lj/coul 0.0 0.0 0.0
The this input achieves the same effect:
@ -265,7 +265,7 @@ The this input achieves the same effect:
.. code-block:: LAMMPS
special_bonds 0.0 0.0 0.1
pair_hybrid lj/charmm/coul/long 8.0 10.0 lj/cut/coul/long 10.0
pair_style hybrid lj/charmm/coul/long 8.0 10.0 lj/cut/coul/long 10.0
pair_modify pair lj/cut/coul/long special lj 0.0 0.0 0.5
pair_modify pair lj/cut/coul/long special coul 0.0 0.0 0.83333333
pair_modify pair lj/charmm/coul/long special lj/coul 0.0 0.0 0.0
@ -279,7 +279,7 @@ effectively *lj/coul 0.0 0.0 0.5* as required for OPLS/AA:
.. code-block:: LAMMPS
special_bonds lj/coul 1e-20 1e-20 0.5
pair_hybrid tersoff lj/cut/coul/long 12.0
pair_style hybrid tersoff lj/cut/coul/long 12.0
pair_modify pair tersoff special lj/coul 1.0 1.0 1.0
For use with the various :doc:`compute \*/tally <compute_tally>`
@ -417,7 +417,11 @@ assigned automatically to the sub-style defined for both I,I and J,J and
its coefficients generated by the mixing rule used by that sub-style.
For the *hybrid/overlay* and *hybrid/scaled* style, there is an
additional requirement that both the I,I and J,J pairs are assigned to a
single sub-style. See the :doc:`pair_modify <pair_modify>` command for
single sub-style. If this requirement is not met, no I,J coeffs will be
generated, even if the sub-styles support mixing, and I,J pair
coefficients must be explicitly defined.
See the :doc:`pair_modify <pair_modify>` command for
details of mixing rules. See the See the doc page for the sub-style to
see if allows for mixing.

46
doc/src/thermo_modify.rst
View File

@ -11,12 +11,13 @@ Syntax
thermo_modify keyword value ...
* one or more keyword/value pairs may be listed
* keyword = *lost* or *lost/bond* or *norm* or *flush* or *line* or *format* or *temp* or *press*
* keyword = *lost* or *lost/bond* or *warn* or *norm* or *flush* or *line* or *format* or *temp* or *press*
.. parsed-literal::
*lost* value = *error* or *warn* or *ignore*
*lost/bond* value = *error* or *warn* or *ignore*
*warn* value = *ignore* or *reset* or *default* or a number
*norm* value = *yes* or *no*
*flush* value = *yes* or *no*
*line* value = *one* or *multi*
@ -75,6 +76,43 @@ are drifting out of the box through a fixed boundary condition (see
the :doc:`boundary <boundary>` command). In this case one atom may be
deleted before the rest of the molecule is, on a later timestep.
The *warn* keyword allows you to control whether LAMMPS will print
warning messages and how many of them. Most warning messages are only
printed by MPI rank 0. They are usually pointing out important issues
that should be investigated, but LAMMPS cannot determine for
certain whether they are an indication of an error.
Some warning messages are printed during a run (or immediately before)
each time a specific MPI rank encounters the issue, e.g. bonds that are
stretched too far or dihedrals in extreme configurations. These number
of these can quickly blow up the size of the log file and screen output.
Thus a limit of 100 warning messages is applied by default. The warning
count is applied to the entire input unless reset with a ``thermo_modify
warn reset`` command. If there are more warnings than the limit, LAMMPS
will print one final warning that it will not print any additional
warning messages.
.. note::
The warning limit is enforced on either the per-processor count or
the total count across all processors. For efficiency reasons,
however, the total count is only updated at steps with thermodynamic
output. Thus when running on a large number of processors in
parallel, the total number of warnings printed can be significantly
larger than the given limit.
Any number after the keyword *warn* will change the warning limit
accordingly. With the value *ignore* all warnings will be suppressed,
with the value *always* no limit will be applied and warnings will
always be printed, with the value *reset* the internal warning counter
will be reset to zero, and with the value *default*, the counter is
reset and the limit set to 100. An example usage of either *reset* or
*default* would be to re-enable warnings that were disabled or have
reached the limit during equilibration, where the warnings would be
acceptable while the system is still adjusting, but then change
to all warnings for the production run, where they would indicate
problems that would require a closer look at what is causing them.
The *norm* keyword determines whether various thermodynamic output
values are normalized by the number of atoms or not, depending on
whether it is set to *yes* or *no*\ . Different unit styles have
@ -183,9 +221,9 @@ Related commands
Default
"""""""
The option defaults are lost = error, norm = yes for unit style of
*lj*\ , norm = no for unit style of *real* and *metal*\ , flush = no,
and temp/press = compute IDs defined by thermo_style.
The option defaults are lost = error, warn = 100, norm = yes for unit
style of *lj*\ , norm = no for unit style of *real* and *metal*\ ,
flush = no, and temp/press = compute IDs defined by thermo_style.
The defaults for the line and format options depend on the thermo
style. For styles "one" and "custom", the line and format defaults

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

@ -1024,6 +1024,7 @@ fmag
fmass
fmm
fmt
fmtlib
fmx
fmy
fmz

2
examples/python/in.fix_python_invoke_neighlist
View File

@ -33,7 +33,7 @@ def post_force_callback(lmp, v):
print(pid_prefix, "### POST_FORCE ###", t)
#mylist = L.numpy.get_neighlist(0)
idx = L.find_pair_neighlist("lj/cut", request=0)
idx = L.find_pair_neighlist("lj/cut", reqid=0)
mylist = L.numpy.get_neighlist(idx)
print(pid_prefix, mylist)
nlocal = L.extract_global("nlocal")

16
examples/python/py_nve.py
View File

@ -1,7 +1,5 @@
from __future__ import print_function
from lammps import lammps, LAMMPS_INT, LAMMPS_DOUBLE
import ctypes
import traceback
from lammps import lammps
import numpy as np
class LAMMPSFix(object):
@ -80,36 +78,30 @@ class NVE_Opt(LAMMPSFixMove):
self.ntypes = self.lmp.extract_global("ntypes")
self.dtv = dt
self.dtf = 0.5 * dt * ftm2v
self.mass = self.lmp.numpy.extract_atom("mass")
def initial_integrate(self, vflag):
nlocal = self.lmp.extract_global("nlocal")
mass = self.lmp.numpy.extract_atom("mass")
atype = self.lmp.numpy.extract_atom("type")
x = self.lmp.numpy.extract_atom("x")
v = self.lmp.numpy.extract_atom("v")
f = self.lmp.numpy.extract_atom("f")
dtf = self.dtf
dtv = self.dtv
mass = self.mass
dtfm = dtf / np.take(mass, atype)
dtfm.reshape((nlocal, 1))
for d in range(x.shape[1]):
v[:,d] += dtfm[:,0] * f[:,d]
v[:,d] += dtfm * f[:,d]
x[:,d] += dtv * v[:,d]
def final_integrate(self):
nlocal = self.lmp.extract_global("nlocal")
mass = self.lmp.numpy.extract_atom("mass")
atype = self.lmp.numpy.extract_atom("type")
v = self.lmp.numpy.extract_atom("v")
f = self.lmp.numpy.extract_atom("f")
dtf = self.dtf
mass = self.mass
dtfm = dtf / np.take(mass, atype)
dtfm.reshape((nlocal, 1))
for d in range(v.shape[1]):
v[:,d] += dtfm[:,0] * f[:,d]
v[:,d] += dtfm * f[:,d]

4
lib/gpu/geryon/nvd_device.h
View File

@ -24,10 +24,6 @@
#ifndef NVD_DEVICE
#define NVD_DEVICE
// workaround after GPU package Feb2021 update
// todo: make new neighbor code work with CUDA
#define LAL_USE_OLD_NEIGHBOR
#include <string>
#include <vector>
#include <iostream>

2
lib/gpu/geryon/ocl_device.h
View File

@ -29,7 +29,7 @@
#include <iostream>
#ifndef CL_TARGET_OPENCL_VERSION
#define CL_TARGET_OPENCL_VERSION 210
#define CL_TARGET_OPENCL_VERSION 300
#endif
#ifdef __APPLE__

2
lib/gpu/geryon/ocl_macros.h
View File

@ -5,7 +5,7 @@
#include <cassert>
#ifndef CL_TARGET_OPENCL_VERSION
#define CL_TARGET_OPENCL_VERSION 210
#define CL_TARGET_OPENCL_VERSION 300
#endif
#ifdef __APPLE__

18
lib/gpu/lal_device.cpp
View File

@ -18,6 +18,7 @@
#include <map>
#include <cmath>
#include <cstdlib>
#include <iostream>
#if (LAL_USE_OMP == 1)
#include <omp.h>
#endif
@ -777,6 +778,7 @@ void DeviceT::output_times(UCL_Timer &time_pair, Answer<numtyp,acctyp> &ans,
#ifdef USE_OPENCL
// Workaround for timing issue on Intel OpenCL
if (times[3] > 80e6) times[3]=0.0;
if (times[5] > 80e6) times[5]=0.0;
#endif
if (replica_me()==0)
@ -1025,6 +1027,22 @@ Device<PRECISION,ACC_PRECISION> global_device;
}
using namespace LAMMPS_AL;
bool lmp_has_gpu_device()
{
UCL_Device gpu;
return (gpu.num_platforms() > 0);
}
std::string lmp_gpu_device_info()
{
std::ostringstream out;
UCL_Device gpu;
if (gpu.num_platforms() > 0)
gpu.print_all(out);
return out.str();
}
int lmp_init_device(MPI_Comm world, MPI_Comm replica, const int ngpu,
const int first_gpu_id, const int gpu_mode,
const double particle_split, const int t_per_atom,

20
lib/gpu/lal_eam.cu
View File

@ -225,7 +225,7 @@ __kernel void k_energy(const __global numtyp4 *restrict x_,
const numtyp rdr, const numtyp rdrho,
const numtyp rhomax, const int nrho,
const int nr, const int t_per_atom) {
int tid, ii, offset, i, itype;
int tid, ii, offset, i, itype, tfrho;
atom_info(t_per_atom,ii,tid,offset);
int n_stride;
@ -242,6 +242,7 @@ __kernel void k_energy(const __global numtyp4 *restrict x_,
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
itype=ix.w;
tfrho=type2frho[itype];
for ( ; nbor<nbor_end; nbor+=n_stride) {
int j=dev_packed[nbor];
@ -270,7 +271,6 @@ __kernel void k_energy(const __global numtyp4 *restrict x_,
}
} // for nbor
} // if ii
const numtyp tfrho=type2frho[itype];
store_energy_fp(rho,energy,ii,inum,tid,t_per_atom,offset,
eflag,vflag,engv,rdrho,nrho,i,rhomax,tfrho);
}
@ -291,7 +291,7 @@ __kernel void k_energy_fast(const __global numtyp4 *restrict x_,
const numtyp rdrho, const numtyp rhomax,
const int nrho, const int nr,
const int t_per_atom) {
int tid, ii, offset, i, itype;
int tid, ii, offset, i, itype, tfrho;
atom_info(t_per_atom,ii,tid,offset);
#ifndef ONETYPE
@ -305,9 +305,9 @@ __kernel void k_energy_fast(const __global numtyp4 *restrict x_,
}
__syncthreads();
#else
const numtyp type2rhor_z2rx=
const int type2rhor_z2rx=
type2rhor_z2r_in[ONETYPE*MAX_SHARED_TYPES+ONETYPE].x;
const numtyp tfrho=type2frho_in[ONETYPE];
tfrho=type2frho_in[ONETYPE];
#endif
int n_stride;
@ -325,6 +325,7 @@ __kernel void k_energy_fast(const __global numtyp4 *restrict x_,
numtyp4 ix; fetch4(ix,i,pos_tex); //x_[i];
#ifndef ONETYPE
itype=ix.w;
tfrho=type2frho[itype];
#endif
for ( ; nbor<nbor_end; nbor+=n_stride) {
@ -347,7 +348,7 @@ __kernel void k_energy_fast(const __global numtyp4 *restrict x_,
p = MIN(p,(numtyp)1.0);
#ifndef ONETYPE
int jtype=fast_mul((int)MAX_SHARED_TYPES,jx.w);
int jtype = fast_mul((int)MAX_SHARED_TYPES,jx.w);
int mtype = jtype+itype;
int index = type2rhor_z2r[mtype].x*(nr+1)+m;
#else
@ -358,9 +359,6 @@ __kernel void k_energy_fast(const __global numtyp4 *restrict x_,
}
} // for nbor
} // if ii
#ifndef ONETYPE
const numtyp tfrho=type2frho[itype];
#endif
store_energy_fp(rho,energy,ii,inum,tid,t_per_atom,offset,
eflag,vflag,engv,rdrho,nrho,i,rhomax,tfrho);
}
@ -498,8 +496,8 @@ __kernel void k_eam_fast(const __global numtyp4 *x_,
__syncthreads();
#else
const int oi=ONETYPE*MAX_SHARED_TYPES+ONETYPE;
const numtyp type2rhor_z2rx=type2rhor_z2r_in[oi].x;
const numtyp type2rhor_z2ry=type2rhor_z2r_in[oi].y;
const int type2rhor_z2rx=type2rhor_z2r_in[oi].x;
const int type2rhor_z2ry=type2rhor_z2r_in[oi].y;
#endif
int n_stride;

2
lib/gpu/lal_neighbor_gpu.cu
View File

@ -115,7 +115,7 @@ __kernel void kernel_calc_cell_counts(const unsigned *restrict cell_id,
#define tagint int
#endif
#ifdef LAMMPS_BIGBIG
#define tagint long long int
#define tagint long
#endif
#ifdef LAMMPS_SMALLSMALL
#define tagint int

4
lib/pace/Install.py
View File

@ -6,7 +6,7 @@ used to automate the steps described in the README file in this dir
"""
from __future__ import print_function
import sys, os, subprocess, shutil
import sys, subprocess
from argparse import ArgumentParser
sys.path.append('..')
@ -102,4 +102,4 @@ if buildflag:
print("Removing pace build files and archive ...")
cmd = 'rm %s; make clean-build' % (download_filename)
subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)
subprocess.check_output(cmd, stderr=subprocess.STDOUT, shell=True)

6
potentials/CdTe.bop.table
View File

@ -19,11 +19,9 @@
0.00000000E+00 0.00000000E+00
0.50000000E+00 0.00000000E+00 0.10000000E+07
0.39687010E+00 0.88101950E+00 -0.27788970E+00
-0.10071280E+00 0.10000000E+01 0.10071280E+00
0.20952380E+00 0.60000000E+00 0.19047620E+00
0.39284960E-02 0.99992720E+00 -0.38556650E-02
-0.10071280E+00 0.10000000E+01 0.10071280E+00
0.20073220E+00 0.60000000E+00 0.19926780E+00
0.20952380E+00 0.60000000E+00 0.19047620E+00
0.39284960E-02 0.99992720E+00 -0.38556650E-02
0.11718170E+00 0.83481160E+00 0.48006700E-01
0.37236414E+06 0.37236414E+06 0.37236414E+06 0.37236414E+06 0.37236414E+06
@ -4831,3 +4829,5 @@
0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00

9664
potentials/CdTeZn.bop.table Normal file
View File

File diff suppressed because it is too large Load Diff

6
potentials/GaAs.bop.table
View File

@ -19,11 +19,9 @@
0.00000000E+00 0.00000000E+00
0.50000000E+00 0.00000000E+00 0.10E-04
9.21518535E-02 7.20160363E-01 1.87687783E-01
5.37304525E-02 7.20160376E-01 2.26109172E-01
2.34227480E-01 7.65772520E-01 0.00000000E+00
-3.18730636E-02 7.49789590E-01 2.82083474E-01
5.37304525E-02 7.20160376E-01 2.26109172E-01
2.56221496E-01 6.68582598E-01 7.51959058E-02
2.34227480E-01 7.65772520E-01 0.00000000E+00
-3.18730636E-02 7.49789590E-01 2.82083474E-01
-2.11787351E-02 7.66726837E-01 2.54451898E-01
0.67538417E+03 0.67538417E+03 0.67538417E+03 0.67538417E+03 0.67538417E+03
@ -4831,3 +4829,5 @@
0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00
0.00000000E+00

9
python/lammps/__init__.py
View File

@ -28,12 +28,17 @@ def get_version_number():
from importlib.metadata import version
try:
vstring = version('lammps')
except: pass
except:
# nothing to do, ignore
pass
else:
from pkg_resources import get_distribution
try:
vstring = get_distribution('lammps').version
except: pass
except:
# nothing to do, ignore
pass
if not vstring:
return 0

239
python/lammps/core.py
View File

@ -18,9 +18,6 @@ from __future__ import print_function
import os
import sys
import traceback
import types
import warnings
from ctypes import *
from os.path import dirname,abspath,join
from inspect import getsourcefile
@ -47,7 +44,7 @@ class ExceptionCheck:
def __enter__(self):
pass
def __exit__(self, type, value, traceback):
def __exit__(self, exc_type, exc_value, traceback):
if self.lmp.has_exceptions and self.lmp.lib.lammps_has_error(self.lmp.lmp):
raise self.lmp._lammps_exception
@ -284,6 +281,7 @@ class lammps(object):
self.lib.lammps_version.argtypes = [c_void_p]
self.lib.lammps_get_os_info.argtypes = [c_char_p, c_int]
self.lib.lammps_get_gpu_device_info.argtypes = [c_char_p, c_int]
self.lib.lammps_get_mpi_comm.argtypes = [c_void_p]
@ -310,6 +308,7 @@ class lammps(object):
# tested to work with mpi4py versions 2 and 3
self.has_mpi4py = mpi4py_version.split('.')[0] in ['2','3']
except:
# ignore failing import
pass
# if no ptr provided, create an instance of LAMMPS
@ -418,9 +417,16 @@ class lammps(object):
# shut-down LAMMPS instance
def __del__(self):
if self.lmp and self.opened:
self.lib.lammps_close(self.lmp)
self.opened = 0
self.close()
# -------------------------------------------------------------------------
# context manager implementation
def __enter__(self):
return self
def __exit__(self, ex_type, ex_value, ex_traceback):
self.close()
# -------------------------------------------------------------------------
@ -447,7 +453,8 @@ class lammps(object):
This is a wrapper around the :cpp:func:`lammps_close` function of the C-library interface.
"""
if self.opened: self.lib.lammps_close(self.lmp)
if self.lmp and self.opened:
self.lib.lammps_close(self.lmp)
self.lmp = None
self.opened = 0
@ -456,9 +463,7 @@ class lammps(object):
def finalize(self):
"""Shut down the MPI communication through the library interface by calling :cpp:func:`lammps_finalize`.
"""
if self.opened: self.lib.lammps_close(self.lmp)
self.lmp = None
self.opened = 0
self.close()
self.lib.lammps_finalize()
# -------------------------------------------------------------------------
@ -486,7 +491,7 @@ class lammps(object):
sb = create_string_buffer(512)
self.lib.lammps_get_os_info(sb,512)
return sb
return sb.value.decode()
# -------------------------------------------------------------------------
@ -779,6 +784,9 @@ class lammps(object):
target_type = float
elif dtype == LAMMPS_STRING:
self.lib.lammps_extract_global.restype = c_char_p
target_type = str
else:
target_type = None
ptr = self.lib.lammps_extract_global(self.lmp, name)
if ptr:
@ -878,71 +886,71 @@ class lammps(object):
# -------------------------------------------------------------------------
def extract_compute(self,id,style,type):
def extract_compute(self,cid,cstyle,ctype):
"""Retrieve data from a LAMMPS compute
This is a wrapper around the :cpp:func:`lammps_extract_compute`
function of the C-library interface.
This function returns ``None`` if either the compute id is not
recognized, or an invalid combination of :ref:`style <py_style_constants>`
and :ref:`type <py_type_constants>` constants is used. The
recognized, or an invalid combination of :ref:`cstyle <py_style_constants>`
and :ref:`ctype <py_type_constants>` constants is used. The
names and functionality of the constants are the same as for
the corresponding C-library function. For requests to return
a scalar or a size, the value is returned, otherwise a pointer.
:param id: compute ID
:type id: string
:param style: style of the data retrieve (global, atom, or local), see :ref:`py_style_constants`
:type style: int
:param type: type or size of the returned data (scalar, vector, or array), see :ref:`py_type_constants`
:type type: int
:param cid: compute ID
:type cid: string
:param cstyle: style of the data retrieve (global, atom, or local), see :ref:`py_style_constants`
:type cstyle: int
:param ctype: type or size of the returned data (scalar, vector, or array), see :ref:`py_type_constants`
:type ctype: int
:return: requested data as scalar, pointer to 1d or 2d double array, or None
:rtype: c_double, ctypes.POINTER(c_double), ctypes.POINTER(ctypes.POINTER(c_double)), or NoneType
"""
if id: id = id.encode()
if cid: cid = cid.encode()
else: return None
if type == LMP_TYPE_SCALAR:
if ctype == LMP_TYPE_SCALAR:
if style == LMP_STYLE_GLOBAL:
self.lib.lammps_extract_compute.restype = POINTER(c_double)
with ExceptionCheck(self):
ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type)
ptr = self.lib.lammps_extract_compute(self.lmp,cid,cstyle,ctype)
return ptr[0]
elif style == LMP_STYLE_ATOM:
return None
elif style == LMP_STYLE_LOCAL:
self.lib.lammps_extract_compute.restype = POINTER(c_int)
with ExceptionCheck(self):
ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type)
ptr = self.lib.lammps_extract_compute(self.lmp,cid,cstyle,ctype)
return ptr[0]
elif type == LMP_TYPE_VECTOR:
elif ctype == LMP_TYPE_VECTOR:
self.lib.lammps_extract_compute.restype = POINTER(c_double)
with ExceptionCheck(self):
ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type)
ptr = self.lib.lammps_extract_compute(self.lmp,cid,cstyle,ctype)
return ptr
elif type == LMP_TYPE_ARRAY:
elif ctype == LMP_TYPE_ARRAY:
self.lib.lammps_extract_compute.restype = POINTER(POINTER(c_double))
with ExceptionCheck(self):
ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type)
ptr = self.lib.lammps_extract_compute(self.lmp,cid,cstyle,ctype)
return ptr
elif type == LMP_SIZE_COLS:
if style == LMP_STYLE_GLOBAL \
or style == LMP_STYLE_ATOM \
or style == LMP_STYLE_LOCAL:
elif ctype == LMP_SIZE_COLS:
if cstyle == LMP_STYLE_GLOBAL \
or cstyle == LMP_STYLE_ATOM \
or cstyle == LMP_STYLE_LOCAL:
self.lib.lammps_extract_compute.restype = POINTER(c_int)
with ExceptionCheck(self):
ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type)
ptr = self.lib.lammps_extract_compute(self.lmp,cid,cstyle,ctype)
return ptr[0]
elif type == LMP_SIZE_VECTOR or type == LMP_SIZE_ROWS:
if style == LMP_STYLE_GLOBAL \
or style == LMP_STYLE_LOCAL:
elif ctype == LMP_SIZE_VECTOR or ctype == LMP_SIZE_ROWS:
if cstyle == LMP_STYLE_GLOBAL \
or cstyle == LMP_STYLE_LOCAL:
self.lib.lammps_extract_compute.restype = POINTER(c_int)
with ExceptionCheck(self):
ptr = self.lib.lammps_extract_compute(self.lmp,id,style,type)
ptr = self.lib.lammps_extract_compute(self.lmp,cid,cstyle,ctype)
return ptr[0]
return None
@ -952,25 +960,25 @@ class lammps(object):
# in case of global data, free memory for 1 double via lammps_free()
# double was allocated by library interface function
def extract_fix(self,id,style,type,nrow=0,ncol=0):
def extract_fix(self,fid,fstyle,ftype,nrow=0,ncol=0):
"""Retrieve data from a LAMMPS fix
This is a wrapper around the :cpp:func:`lammps_extract_fix`
function of the C-library interface.
This function returns ``None`` if either the fix id is not
recognized, or an invalid combination of :ref:`style <py_style_constants>`
and :ref:`type <py_type_constants>` constants is used. The
recognized, or an invalid combination of :ref:`fstyle <py_style_constants>`
and :ref:`ftype <py_type_constants>` constants is used. The
names and functionality of the constants are the same as for
the corresponding C-library function. For requests to return
a scalar or a size, the value is returned, also when accessing
global vectors or arrays, otherwise a pointer.
:param id: fix ID
:type id: string
:param style: style of the data retrieve (global, atom, or local), see :ref:`py_style_constants`
:type style: int
:param type: type or size of the returned data (scalar, vector, or array), see :ref:`py_type_constants`
:type type: int
:param fid: fix ID
:type fid: string
:param fstyle: style of the data retrieve (global, atom, or local), see :ref:`py_style_constants`
:type fstyle: int
:param ftype: type or size of the returned data (scalar, vector, or array), see :ref:`py_type_constants`
:type ftype: int
:param nrow: index of global vector element or row index of global array element
:type nrow: int
:param ncol: column index of global array element
@ -979,53 +987,53 @@ class lammps(object):
:rtype: c_double, ctypes.POINTER(c_double), ctypes.POINTER(ctypes.POINTER(c_double)), or NoneType
"""
if id: id = id.encode()
if fid: fid = fid.encode()
else: return None
if style == LMP_STYLE_GLOBAL:
if type in (LMP_TYPE_SCALAR, LMP_TYPE_VECTOR, LMP_TYPE_ARRAY):
if fstyle == LMP_STYLE_GLOBAL:
if ftype in (LMP_TYPE_SCALAR, LMP_TYPE_VECTOR, LMP_TYPE_ARRAY):
self.lib.lammps_extract_fix.restype = POINTER(c_double)
with ExceptionCheck(self):
ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,nrow,ncol)
ptr = self.lib.lammps_extract_fix(self.lmp,fid,fstyle,ftype,nrow,ncol)
result = ptr[0]
self.lib.lammps_free(ptr)
return result
elif type in (LMP_SIZE_VECTOR, LMP_SIZE_ROWS, LMP_SIZE_COLS):
elif ftype in (LMP_SIZE_VECTOR, LMP_SIZE_ROWS, LMP_SIZE_COLS):
self.lib.lammps_extract_fix.restype = POINTER(c_int)
with ExceptionCheck(self):
ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,nrow,ncol)
ptr = self.lib.lammps_extract_fix(self.lmp,fid,fstyle,ftype,nrow,ncol)
return ptr[0]
else:
return None
elif style == LMP_STYLE_ATOM:
if type == LMP_TYPE_VECTOR:
elif fstyle == LMP_STYLE_ATOM:
if ftype == LMP_TYPE_VECTOR:
self.lib.lammps_extract_fix.restype = POINTER(c_double)
elif type == LMP_TYPE_ARRAY:
elif ftype == LMP_TYPE_ARRAY:
self.lib.lammps_extract_fix.restype = POINTER(POINTER(c_double))
elif type == LMP_SIZE_COLS:
elif ftype == LMP_SIZE_COLS:
self.lib.lammps_extract_fix.restype = POINTER(c_int)
else:
return None
with ExceptionCheck(self):
ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,nrow,ncol)
if type == LMP_SIZE_COLS:
ptr = self.lib.lammps_extract_fix(self.lmp,fid,fstyle,ftype,nrow,ncol)
if ftype == LMP_SIZE_COLS:
return ptr[0]
else:
return ptr
elif style == LMP_STYLE_LOCAL:
if type == LMP_TYPE_VECTOR:
elif fstyle == LMP_STYLE_LOCAL:
if ftype == LMP_TYPE_VECTOR:
self.lib.lammps_extract_fix.restype = POINTER(c_double)
elif type == LMP_TYPE_ARRAY:
elif ftype == LMP_TYPE_ARRAY:
self.lib.lammps_extract_fix.restype = POINTER(POINTER(c_double))
elif type in (LMP_TYPE_SCALAR, LMP_SIZE_VECTOR, LMP_SIZE_ROWS, LMP_SIZE_COLS):
elif ftype in (LMP_TYPE_SCALAR, LMP_SIZE_VECTOR, LMP_SIZE_ROWS, LMP_SIZE_COLS):
self.lib.lammps_extract_fix.restype = POINTER(c_int)
else:
return None
with ExceptionCheck(self):
ptr = self.lib.lammps_extract_fix(self.lmp,id,style,type,nrow,ncol)
if type in (LMP_TYPE_VECTOR, LMP_TYPE_ARRAY):
ptr = self.lib.lammps_extract_fix(self.lmp,fid,fstyle,ftype,nrow,ncol)
if ftype in (LMP_TYPE_VECTOR, LMP_TYPE_ARRAY):
return ptr
else:
return ptr[0]
@ -1114,20 +1122,20 @@ class lammps(object):
# return vector of atom properties gathered across procs
# 3 variants to match src/library.cpp
# name = atom property recognized by LAMMPS in atom->extract()
# type = 0 for integer values, 1 for double values
# dtype = 0 for integer values, 1 for double values
# count = number of per-atom valus, 1 for type or charge, 3 for x or f
# returned data is a 1d vector - doc how it is ordered?
# NOTE: need to insure are converting to/from correct Python type
# e.g. for Python list or NumPy or ctypes
def gather_atoms(self,name,type,count):
def gather_atoms(self,name,dtype,count):
if name: name = name.encode()
natoms = self.get_natoms()
with ExceptionCheck(self):
if type == 0:
if dtype == 0:
data = ((count*natoms)*c_int)()
self.lib.lammps_gather_atoms(self.lmp,name,type,count,data)
elif type == 1:
elif dtype == 1:
data = ((count*natoms)*c_double)()
self.lib.lammps_gather_atoms(self.lmp,name,type,count,data)
else:
@ -1136,29 +1144,29 @@ class lammps(object):
# -------------------------------------------------------------------------
def gather_atoms_concat(self,name,type,count):
def gather_atoms_concat(self,name,dtype,count):
if name: name = name.encode()
natoms = self.get_natoms()
with ExceptionCheck(self):
if type == 0:
if dtype == 0:
data = ((count*natoms)*c_int)()
self.lib.lammps_gather_atoms_concat(self.lmp,name,type,count,data)
elif type == 1:
elif dtype == 1:
data = ((count*natoms)*c_double)()
self.lib.lammps_gather_atoms_concat(self.lmp,name,type,count,data)
else:
return None
return data
def gather_atoms_subset(self,name,type,count,ndata,ids):
def gather_atoms_subset(self,name,dtype,count,ndata,ids):
if name: name = name.encode()
with ExceptionCheck(self):
if type == 0:
if dtype == 0:
data = ((count*ndata)*c_int)()
self.lib.lammps_gather_atoms_subset(self.lmp,name,type,count,ndata,ids,data)
elif type == 1:
self.lib.lammps_gather_atoms_subset(self.lmp,name,dtype,count,ndata,ids,data)
elif dtype == 1:
data = ((count*ndata)*c_double)()
self.lib.lammps_gather_atoms_subset(self.lmp,name,type,count,ndata,ids,data)
self.lib.lammps_gather_atoms_subset(self.lmp,name,dtype,count,ndata,ids,data)
else:
return None
return data
@ -1174,17 +1182,17 @@ class lammps(object):
# NOTE: need to insure are converting to/from correct Python type
# e.g. for Python list or NumPy or ctypes
def scatter_atoms(self,name,type,count,data):
def scatter_atoms(self,name,dtype,count,data):
if name: name = name.encode()
with ExceptionCheck(self):
self.lib.lammps_scatter_atoms(self.lmp,name,type,count,data)
self.lib.lammps_scatter_atoms(self.lmp,name,dtype,count,data)
# -------------------------------------------------------------------------
def scatter_atoms_subset(self,name,type,count,ndata,ids,data):
def scatter_atoms_subset(self,name,dtype,count,ndata,ids,data):
if name: name = name.encode()
with ExceptionCheck(self):
self.lib.lammps_scatter_atoms_subset(self.lmp,name,type,count,ndata,ids,data)
self.lib.lammps_scatter_atoms_subset(self.lmp,name,dtype,count,ndata,ids,data)
# return vector of atom/compute/fix properties gathered across procs
# 3 variants to match src/library.cpp
@ -1194,43 +1202,43 @@ class lammps(object):
# returned data is a 1d vector - doc how it is ordered?
# NOTE: need to insure are converting to/from correct Python type
# e.g. for Python list or NumPy or ctypes
def gather(self,name,type,count):
def gather(self,name,dtype,count):
if name: name = name.encode()
natoms = self.get_natoms()
with ExceptionCheck(self):
if type == 0:
if dtype == 0:
data = ((count*natoms)*c_int)()
self.lib.lammps_gather(self.lmp,name,type,count,data)
elif type == 1:
self.lib.lammps_gather(self.lmp,name,dtype,count,data)
elif dtype == 1:
data = ((count*natoms)*c_double)()
self.lib.lammps_gather(self.lmp,name,type,count,data)
self.lib.lammps_gather(self.lmp,name,dtype,count,data)
else:
return None
return data
def gather_concat(self,name,type,count):
def gather_concat(self,name,dtype,count):
if name: name = name.encode()
natoms = self.get_natoms()
with ExceptionCheck(self):
if type == 0:
if dtype == 0:
data = ((count*natoms)*c_int)()
self.lib.lammps_gather_concat(self.lmp,name,type,count,data)
elif type == 1:
self.lib.lammps_gather_concat(self.lmp,name,dtype,count,data)
elif dtype == 1:
data = ((count*natoms)*c_double)()
self.lib.lammps_gather_concat(self.lmp,name,type,count,data)
self.lib.lammps_gather_concat(self.lmp,name,dtype,count,data)
else:
return None
return data
def gather_subset(self,name,type,count,ndata,ids):
def gather_subset(self,name,dtype,count,ndata,ids):
if name: name = name.encode()
with ExceptionCheck(self):
if type == 0:
if dtype == 0:
data = ((count*ndata)*c_int)()
self.lib.lammps_gather_subset(self.lmp,name,type,count,ndata,ids,data)
elif type == 1:
self.lib.lammps_gather_subset(self.lmp,name,dtype,count,ndata,ids,data)
elif dtype == 1:
data = ((count*ndata)*c_double)()
self.lib.lammps_gather_subset(self.lmp,name,type,count,ndata,ids,data)
self.lib.lammps_gather_subset(self.lmp,name,dtype,count,ndata,ids,data)
else:
return None
return data
@ -1244,15 +1252,15 @@ class lammps(object):
# NOTE: need to insure are converting to/from correct Python type
# e.g. for Python list or NumPy or ctypes
def scatter(self,name,type,count,data):
def scatter(self,name,dtype,count,data):
if name: name = name.encode()
with ExceptionCheck(self):
self.lib.lammps_scatter(self.lmp,name,type,count,data)
self.lib.lammps_scatter(self.lmp,name,dtype,count,data)
def scatter_subset(self,name,type,count,ndata,ids,data):
def scatter_subset(self,name,dtype,count,ndata,ids,data):
if name: name = name.encode()
with ExceptionCheck(self):
self.lib.lammps_scatter_subset(self.lmp,name,type,count,ndata,ids,data)
self.lib.lammps_scatter_subset(self.lmp,name,dtype,count,ndata,ids,data)
# -------------------------------------------------------------------------
@ -1545,6 +1553,37 @@ class lammps(object):
# -------------------------------------------------------------------------
@property
def has_gpu_device(self):
""" Availability of GPU package compatible device
This is a wrapper around the :cpp:func:`lammps_has_gpu_device`
function of the C library interface.
:return: True if a GPU package compatible device is present, otherwise False
:rtype: bool
"""
return self.lib.lammps_has_gpu_device() != 0
# -------------------------------------------------------------------------
def get_gpu_device_info(self):
"""Return a string with detailed information about any devices that are
usable by the GPU package.
This is a wrapper around the :cpp:func:`lammps_get_gpu_device_info`
function of the C-library interface.
:return: GPU device info string
:rtype: string
"""
sb = create_string_buffer(8192)
self.lib.lammps_get_gpu_device_info(sb,8192)
return sb.value.decode()
# -------------------------------------------------------------------------
@property
def installed_packages(self):
""" List of the names of enabled packages in the LAMMPS shared library

2
python/lammps/data.py
View File

@ -16,7 +16,7 @@
# Written by Richard Berger <richard.berger@temple.edu>
################################################################################
class NeighList:
class NeighList(object):
"""This is a wrapper class that exposes the contents of a neighbor list.
It can be used like a regular Python list. Each element is a tuple of:

2
python/lammps/formats.py
View File

@ -176,7 +176,7 @@ class AvgChunkFile:
current[data_column] = [value]
chunks_read += 1
assert (chunk == chunks_read)
assert chunk == chunks_read
else:
# do not support changing number of chunks
if not (num_chunks == int(parts[1])):

84
python/lammps/numpy_wrapper.py
View File

@ -142,7 +142,7 @@ class numpy_wrapper:
# -------------------------------------------------------------------------
def extract_compute(self, cid, style, type):
def extract_compute(self, cid, cstyle, ctype):
"""Retrieve data from a LAMMPS compute
This is a wrapper around the
@ -150,50 +150,50 @@ class numpy_wrapper:
It behaves the same as the original method, but returns NumPy arrays
instead of ``ctypes`` pointers.
:param id: compute ID
:type id: string
:param style: style of the data retrieve (global, atom, or local), see :ref:`py_style_constants`
:type style: int
:param type: type of the returned data (scalar, vector, or array), see :ref:`py_type_constants`
:type type: int
:param cid: compute ID
:type cid: string
:param cstyle: style of the data retrieve (global, atom, or local), see :ref:`py_style_constants`
:type cstyle: int
:param ctype: type of the returned data (scalar, vector, or array), see :ref:`py_type_constants`
:type ctype: int
:return: requested data either as float, as NumPy array with direct access to C data, or None
:rtype: float, numpy.array, or NoneType
"""
value = self.lmp.extract_compute(cid, style, type)
value = self.lmp.extract_compute(cid, cstyle, ctype)
if style in (LMP_STYLE_GLOBAL, LMP_STYLE_LOCAL):
if type == LMP_TYPE_VECTOR:
nrows = self.lmp.extract_compute(cid, style, LMP_SIZE_VECTOR)
if cstyle in (LMP_STYLE_GLOBAL, LMP_STYLE_LOCAL):
if ctype == LMP_TYPE_VECTOR:
nrows = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_VECTOR)
return self.darray(value, nrows)
elif type == LMP_TYPE_ARRAY:
nrows = self.lmp.extract_compute(cid, style, LMP_SIZE_ROWS)
ncols = self.lmp.extract_compute(cid, style, LMP_SIZE_COLS)
elif ctype == LMP_TYPE_ARRAY:
nrows = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_ROWS)
ncols = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_COLS)
return self.darray(value, nrows, ncols)
elif style == LMP_STYLE_ATOM:
if type == LMP_TYPE_VECTOR:
elif cstyle == LMP_STYLE_ATOM:
if ctype == LMP_TYPE_VECTOR:
nlocal = self.lmp.extract_global("nlocal")
return self.darray(value, nlocal)
elif type == LMP_TYPE_ARRAY:
elif ctype == LMP_TYPE_ARRAY:
nlocal = self.lmp.extract_global("nlocal")
ncols = self.lmp.extract_compute(cid, style, LMP_SIZE_COLS)
ncols = self.lmp.extract_compute(cid, cstyle, LMP_SIZE_COLS)
return self.darray(value, nlocal, ncols)
return value
# -------------------------------------------------------------------------
def extract_fix(self, fid, style, type, nrow=0, ncol=0):
def extract_fix(self, fid, fstyle, ftype, nrow=0, ncol=0):
"""Retrieve data from a LAMMPS fix
This is a wrapper around the :py:meth:`lammps.extract_fix() <lammps.lammps.extract_fix()>` method.
It behaves the same as the original method, but returns NumPy arrays
instead of ``ctypes`` pointers.
:param id: fix ID
:type id: string
:param style: style of the data retrieve (global, atom, or local), see :ref:`py_style_constants`
:type style: int
:param type: type or size of the returned data (scalar, vector, or array), see :ref:`py_type_constants`
:type type: int
:param fid: fix ID
:type fid: string
:param fstyle: style of the data retrieve (global, atom, or local), see :ref:`py_style_constants`
:type fstyle: int
:param ftype: type or size of the returned data (scalar, vector, or array), see :ref:`py_type_constants`
:type ftype: int
:param nrow: index of global vector element or row index of global array element
:type nrow: int
:param ncol: column index of global array element
@ -202,22 +202,22 @@ class numpy_wrapper:
:rtype: integer or double value, pointer to 1d or 2d double array or None
"""
value = self.lmp.extract_fix(fid, style, type, nrow, ncol)
if style == LMP_STYLE_ATOM:
if type == LMP_TYPE_VECTOR:
value = self.lmp.extract_fix(fid, fstyle, ftype, nrow, ncol)
if fstyle == LMP_STYLE_ATOM:
if ftype == LMP_TYPE_VECTOR:
nlocal = self.lmp.extract_global("nlocal")
return self.darray(value, nlocal)
elif type == LMP_TYPE_ARRAY:
elif ftype == LMP_TYPE_ARRAY:
nlocal = self.lmp.extract_global("nlocal")
ncols = self.lmp.extract_fix(fid, style, LMP_SIZE_COLS, 0, 0)
ncols = self.lmp.extract_fix(fid, fstyle, LMP_SIZE_COLS, 0, 0)
return self.darray(value, nlocal, ncols)
elif style == LMP_STYLE_LOCAL:
if type == LMP_TYPE_VECTOR:
nrows = self.lmp.extract_fix(fid, style, LMP_SIZE_ROWS, 0, 0)
elif fstyle == LMP_STYLE_LOCAL:
if ftype == LMP_TYPE_VECTOR:
nrows = self.lmp.extract_fix(fid, fstyle, LMP_SIZE_ROWS, 0, 0)
return self.darray(value, nrows)
elif type == LMP_TYPE_ARRAY:
nrows = self.lmp.extract_fix(fid, style, LMP_SIZE_ROWS, 0, 0)
ncols = self.lmp.extract_fix(fid, style, LMP_SIZE_COLS, 0, 0)
elif ftype == LMP_TYPE_ARRAY:
nrows = self.lmp.extract_fix(fid, fstyle, LMP_SIZE_ROWS, 0, 0)
ncols = self.lmp.extract_fix(fid, fstyle, LMP_SIZE_COLS, 0, 0)
return self.darray(value, nrows, ncols)
return value
@ -293,7 +293,11 @@ class numpy_wrapper:
ptr = cast(raw_ptr[0], POINTER(c_int_type * nelem * dim))
a = np.frombuffer(ptr.contents, dtype=np_int_type)
a.shape = (nelem, dim)
if dim > 1:
a.shape = (nelem, dim)
else:
a.shape = (nelem)
return a
# -------------------------------------------------------------------------
@ -306,7 +310,11 @@ class numpy_wrapper:
ptr = cast(raw_ptr[0], POINTER(c_double * nelem * dim))
a = np.frombuffer(ptr.contents)
a.shape = (nelem, dim)
if dim > 1:
a.shape = (nelem, dim)
else:
a.shape = (nelem)
return a
# -------------------------------------------------------------------------

10
python/lammps/pylammps.py
View File

@ -23,7 +23,6 @@ from __future__ import print_function
import os
import re
import select
import sys
from collections import namedtuple
from .core import lammps
@ -41,7 +40,7 @@ class OutputCapture(object):
os.dup2(self.stdout_pipe_write, self.stdout_fd)
return self
def __exit__(self, type, value, tracebac):
def __exit__(self, exc_type, exc_value, traceback):
os.dup2(self.stdout, self.stdout_fd)
os.close(self.stdout)
os.close(self.stdout_pipe_read)
@ -351,6 +350,7 @@ def get_thermo_data(output):
for i, col in enumerate(columns):
current_run[col].append(values[i])
except ValueError:
# cannot convert. must be a non-thermo output. ignore.
pass
return runs
@ -409,6 +409,12 @@ class PyLammps(object):
self._enable_cmd_history = False
self.runs = []
def __enter__(self):
return self
def __exit__(self, ex_type, ex_value, ex_traceback):
self.close()
def __del__(self):
if self.lmp: self.lmp.close()
self.lmp = None

8
src/.gitignore vendored
View File

@ -923,8 +923,6 @@
/pair_dipole_cut.h
/pair_dipole_sf.cpp
/pair_dipole_sf.h
/pair_dpd_mt.cpp
/pair_dpd_mt.h
/pair_dsmc.cpp
/pair_dsmc.h
/pair_e3b.cpp
@ -1400,8 +1398,10 @@
/pair_thole.h
/pair_buck_mdf.cpp
/pair_buck_mdf.h
/pair_dpd_conservative.cpp
/pair_dpd_conservative.h
/pair_dpd_ext.cpp
/pair_dpd_ext.h
/pair_dpd_ext_tstat.cpp
/pair_dpd_ext_tstat.h
/pair_dpd_fdt.cpp
/pair_dpd_fdt.h
/pair_dpd_fdt_energy.cpp

15
src/CLASS2/dihedral_class2.cpp
View File

@ -207,19 +207,8 @@ void DihedralClass2::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me = comm->me;
if (screen) {
error->warning(FLERR,fmt::format("Dihedral problem: {} {} {} {} {} {}",
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],
atom->tag[i3],atom->tag[i4]));
fmt::print(screen," 1st atom: {} {} {} {}\n",me,x[i1][0],x[i1][1],x[i1][2]);
fmt::print(screen," 2nd atom: {} {} {} {}\n",me,x[i2][0],x[i2][1],x[i2][2]);
fmt::print(screen," 3rd atom: {} {} {} {}\n",me,x[i3][0],x[i3][1],x[i3][2]);
fmt::print(screen," 4th atom: {} {} {} {}\n",me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;

26
src/CLASS2/improper_class2.cpp
View File

@ -18,18 +18,17 @@
#include "improper_class2.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "neighbor.h"
#include "update.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -154,19 +153,8 @@ void ImproperClass2::compute(int eflag, int vflag)
// angle error check
for (i = 0; i < 3; i++) {
if (costheta[i] == -1.0) {
int me = comm->me;
if (screen) {
error->warning(FLERR,fmt::format("Improper problem: {} {} {} {} {} {}",
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],
atom->tag[i3],atom->tag[i4]));
fmt::print(screen," 1st atom: {} {} {} {}\n",me,x[i1][0],x[i1][1],x[i1][2]);
fmt::print(screen," 2nd atom: {} {} {} {}\n",me,x[i2][0],x[i2][1],x[i2][2]);
fmt::print(screen," 3rd atom: {} {} {} {}\n",me,x[i3][0],x[i3][1],x[i3][2]);
fmt::print(screen," 4th atom: {} {} {} {}\n",me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (costheta[i] == -1.0)
problem(FLERR, i1, i2, i3, i4);
}
for (i = 0; i < 3; i++) {

2
src/COMPRESS/dump_atom_gz.cpp
View File

@ -73,7 +73,7 @@ void DumpAtomGZ::openfile()
++numfiles;
} else {
if (remove(nameslist[fileidx]) != 0) {
error->warning(FLERR, fmt::format("Could not delete {}", nameslist[fileidx]));
error->warning(FLERR, "Could not delete {}", nameslist[fileidx]);
}
delete[] nameslist[fileidx];
nameslist[fileidx] = utils::strdup(filecurrent);

2
src/COMPRESS/dump_atom_zstd.cpp
View File

@ -81,7 +81,7 @@ void DumpAtomZstd::openfile()
++numfiles;
} else {
if (remove(nameslist[fileidx]) != 0) {
error->warning(FLERR, fmt::format("Could not delete {}", nameslist[fileidx]));
error->warning(FLERR, "Could not delete {}", nameslist[fileidx]);
}
delete[] nameslist[fileidx];
nameslist[fileidx] = utils::strdup(filecurrent);

2
src/COMPRESS/dump_cfg_gz.cpp
View File

@ -75,7 +75,7 @@ void DumpCFGGZ::openfile()
++numfiles;
} else {
if (remove(nameslist[fileidx]) != 0) {
error->warning(FLERR, fmt::format("Could not delete {}", nameslist[fileidx]));
error->warning(FLERR, "Could not delete {}", nameslist[fileidx]);
}
delete[] nameslist[fileidx];
nameslist[fileidx] = utils::strdup(filecurrent);

2
src/COMPRESS/dump_cfg_zstd.cpp
View File

@ -83,7 +83,7 @@ void DumpCFGZstd::openfile()
++numfiles;
} else {
if (remove(nameslist[fileidx]) != 0) {
error->warning(FLERR, fmt::format("Could not delete {}", nameslist[fileidx]));
error->warning(FLERR, "Could not delete {}", nameslist[fileidx]);
}
delete[] nameslist[fileidx];
nameslist[fileidx] = utils::strdup(filecurrent);

2
src/COMPRESS/dump_custom_gz.cpp
View File

@ -73,7 +73,7 @@ void DumpCustomGZ::openfile()
++numfiles;
} else {
if (remove(nameslist[fileidx]) != 0) {
error->warning(FLERR, fmt::format("Could not delete {}", nameslist[fileidx]));
error->warning(FLERR, "Could not delete {}", nameslist[fileidx]);
}
delete[] nameslist[fileidx];
nameslist[fileidx] = utils::strdup(filecurrent);

2
src/COMPRESS/dump_custom_zstd.cpp
View File

@ -80,7 +80,7 @@ void DumpCustomZstd::openfile()
++numfiles;
} else {
if (remove(nameslist[fileidx]) != 0) {
error->warning(FLERR, fmt::format("Could not delete {}", nameslist[fileidx]));
error->warning(FLERR, "Could not delete {}", nameslist[fileidx]);
}
delete[] nameslist[fileidx];
nameslist[fileidx] = utils::strdup(filecurrent);

2
src/COMPRESS/dump_local_gz.cpp
View File

@ -73,7 +73,7 @@ void DumpLocalGZ::openfile()
++numfiles;
} else {
if (remove(nameslist[fileidx]) != 0) {
error->warning(FLERR, fmt::format("Could not delete {}", nameslist[fileidx]));
error->warning(FLERR, "Could not delete {}", nameslist[fileidx]);
}
delete[] nameslist[fileidx];
nameslist[fileidx] = utils::strdup(filecurrent);

2
src/COMPRESS/dump_local_zstd.cpp
View File

@ -80,7 +80,7 @@ void DumpLocalZstd::openfile()
++numfiles;
} else {
if (remove(nameslist[fileidx]) != 0) {
error->warning(FLERR, fmt::format("Could not delete {}", nameslist[fileidx]));
error->warning(FLERR, "Could not delete {}", nameslist[fileidx]);
}
delete[] nameslist[fileidx];
nameslist[fileidx] = utils::strdup(filecurrent);

5
src/GPU/pair_lj_charmm_coul_charmm_gpu.cpp
View File

@ -134,13 +134,12 @@ void PairLJCharmmCoulCharmmGPU::init_style()
error->all(FLERR,
"Cannot use newton pair with lj/charmm/coul/long/gpu pair style");
// Repeat cutsq calculation because done after call to init_style
// Repeated cutsq calculation in init_one() is required for GPU package
double cut;
for (int i = 1; i <= atom->ntypes; i++) {
for (int j = i; j <= atom->ntypes; j++) {
if (setflag[i][j] != 0 || (setflag[i][i] != 0 && setflag[j][j] != 0))
cut = init_one(i,j);
init_one(i,j);
}
}

41
src/GRANULAR/fix_pour.cpp
View File

@ -14,26 +14,26 @@
#include "fix_pour.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "atom_vec.h"
#include "force.h"
#include "update.h"
#include "comm.h"
#include "molecule.h"
#include "modify.h"
#include "fix_gravity.h"
#include "domain.h"
#include "error.h"
#include "fix_gravity.h"
#include "force.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "modify.h"
#include "molecule.h"
#include "random_park.h"
#include "region.h"
#include "region_block.h"
#include "region_cylinder.h"
#include "random_park.h"
#include "math_extra.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -394,10 +394,12 @@ void FixPour::pre_exchange()
if (next_reneighbor != update->ntimestep) return;
// clear ghost count and any ghost bonus data internal to AtomVec
// clear ghost count (and atom map) and any ghost bonus data
// internal to AtomVec
// same logic as beginning of Comm::exchange()
// do it now b/c inserting atoms will overwrite ghost atoms
if (atom->map_style != Atom::MAP_NONE) atom->map_clear();
atom->nghost = 0;
atom->avec->clear_bonus();
@ -701,7 +703,7 @@ void FixPour::pre_exchange()
int ninserted_mols = ninserted_atoms / natom;
ninserted += ninserted_mols;
if (ninserted_mols < nnew && me == 0)
error->warning(FLERR,"Less insertions than requested",0);
error->warning(FLERR,"Less insertions than requested");
// reset global natoms,nbonds,etc
// increment maxtag_all and maxmol_all if necessary
@ -721,10 +723,13 @@ void FixPour::pre_exchange()
}
if (maxtag_all >= MAXTAGINT)
error->all(FLERR,"New atom IDs exceed maximum allowed ID");
if (atom->map_style != Atom::MAP_NONE) {
atom->map_init();
atom->map_set();
}
}
// rebuild atom map
if (atom->map_style != Atom::MAP_NONE) {
if (success) atom->map_init();
atom->map_set();
}
// free local memory

18
src/GRANULAR/fix_wall_gran_region.cpp
View File

@ -99,21 +99,17 @@ void FixWallGranRegion::init()
// check if region properties changed between runs
// reset if restart info was inconsistent
if (strcmp(idregion,region->id) != 0 ||
strcmp(region_style,region->style) != 0 ||
nregion != region->nregion) {
char str[256];
snprintf(str,256,"Region properties for region %s changed between runs, "
"resetting its motion",idregion);
error->warning(FLERR,str);
if ((strcmp(idregion,region->id) != 0)
|| (strcmp(region_style,region->style) != 0)
|| (nregion != region->nregion)) {
error->warning(FLERR,"Region properties for region {} changed between "
"runs, resetting its motion",idregion);
region->reset_vel();
}
if (motion_resetflag) {
char str[256];
snprintf(str,256,"Region properties for region %s are inconsistent "
"with restart file, resetting its motion",idregion);
error->warning(FLERR,str);
error->warning(FLERR,"Region properties for region {} are inconsistent "
"with restart file, resetting its motion",idregion);
region->reset_vel();
}
}

34
src/GRANULAR/pair_gran_hooke_history.cpp
View File

@ -18,22 +18,22 @@
#include "pair_gran_hooke_history.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "comm.h"
#include "error.h"
#include "fix.h"
#include "fix_dummy.h"
#include "fix_neigh_history.h"
#include "comm.h"
#include "neighbor.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
@ -434,16 +434,10 @@ void PairGranHookeHistory::init_style()
// it replaces FixDummy, created in the constructor
// this is so its order in the fix list is preserved
if (history && fix_history == nullptr) {
char dnumstr[16];
sprintf(dnumstr,"%d",size_history);
char **fixarg = new char*[4];
fixarg[0] = (char *) "NEIGH_HISTORY_HH";
fixarg[1] = (char *) "all";
fixarg[2] = (char *) "NEIGH_HISTORY";
fixarg[3] = dnumstr;
modify->replace_fix("NEIGH_HISTORY_HH_DUMMY",4,fixarg,1);
delete [] fixarg;
if (history && (fix_history == nullptr)) {
auto cmd = fmt::format("NEIGH_HISTORY_HH all NEIGH_HISTORY {}",
size_history);
modify->replace_fix("NEIGH_HISTORY_HH_DUMMY",cmd,1);
int ifix = modify->find_fix("NEIGH_HISTORY_HH");
fix_history = (FixNeighHistory *) modify->fix[ifix];
fix_history->pair = this;

11
src/GRANULAR/pair_granular.cpp
View File

@ -1214,13 +1214,10 @@ double PairGranular::init_one(int i, int j)
(tangential_model[i][i] != tangential_model[j][j]) ||
(roll_model[i][i] != roll_model[j][j]) ||
(twist_model[i][i] != twist_model[j][j])) {
char str[512];
sprintf(str,"Granular pair style functional forms are different, "
"cannot mix coefficients for types %d and %d. \n"
"This combination must be set explicitly "
"via pair_coeff command",i,j);
error->one(FLERR,str);
error->all(FLERR,"Granular pair style functional forms are different, "
"cannot mix coefficients for types {} and {}. \n"
"This combination must be set explicitly via a "
"pair_coeff command",i,j);
}
if (normal_model[i][j] == HERTZ || normal_model[i][j] == HOOKE)

9
src/KOKKOS/bond_fene_kokkos.cpp
View File

@ -23,13 +23,14 @@
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory_kokkos.h"
#include "neighbor_kokkos.h"
#include <cmath>
using namespace LAMMPS_NS;
using MathConst::MY_CUBEROOT2;
/* ---------------------------------------------------------------------- */
@ -129,7 +130,7 @@ void BondFENEKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
k_warning_flag.template modify<DeviceType>();
k_warning_flag.template sync<LMPHostType>();
if (h_warning_flag())
error->warning(FLERR,"FENE bond too long",0);
error->warning(FLERR,"FENE bond too long");
k_error_flag.template modify<DeviceType>();
k_error_flag.template sync<LMPHostType>();
@ -200,7 +201,7 @@ void BondFENEKokkos<DeviceType>::operator()(TagBondFENECompute<NEWTON_BOND,EVFLA
// force from LJ term
F_FLOAT sr6 = 0.0;
if (rsq < TWO_1_3*d_sigma[type]*d_sigma[type]) {
if (rsq < MY_CUBEROOT2*d_sigma[type]*d_sigma[type]) {
const F_FLOAT sr2 = d_sigma[type]*d_sigma[type]/rsq;
sr6 = sr2*sr2*sr2;
fbond += 48.0*d_epsilon[type]*sr6*(sr6-0.5)/rsq;
@ -211,7 +212,7 @@ void BondFENEKokkos<DeviceType>::operator()(TagBondFENECompute<NEWTON_BOND,EVFLA
F_FLOAT ebond = 0.0;
if (eflag) {
ebond = -0.5 * d_k[type]*r0sq*log(rlogarg);
if (rsq < TWO_1_3*d_sigma[type]*d_sigma[type])
if (rsq < MY_CUBEROOT2*d_sigma[type]*d_sigma[type])
ebond += 4.0*d_epsilon[type]*sr6*(sr6-1.0) + d_epsilon[type];
}

2
src/KOKKOS/dihedral_charmm_kokkos.cpp
View File

@ -143,7 +143,7 @@ void DihedralCharmmKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
k_warning_flag.template modify<DeviceType>();
k_warning_flag.template sync<LMPHostType>();
if (h_warning_flag())
error->warning(FLERR,"Dihedral problem",0);
error->warning(FLERR,"Dihedral problem");
if (eflag_global) {
energy += evm.emol;

2
src/KOKKOS/dihedral_class2_kokkos.cpp
View File

@ -166,7 +166,7 @@ void DihedralClass2Kokkos<DeviceType>::compute(int eflag_in, int vflag_in)
k_warning_flag.template modify<DeviceType>();
k_warning_flag.template sync<LMPHostType>();
if (h_warning_flag())
error->warning(FLERR,"Dihedral problem",0);
error->warning(FLERR,"Dihedral problem");
if (eflag_global) energy += ev.evdwl;
if (vflag_global) {

2
src/KOKKOS/dihedral_harmonic_kokkos.cpp
View File

@ -129,7 +129,7 @@ void DihedralHarmonicKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
k_warning_flag.template modify<DeviceType>();
k_warning_flag.template sync<LMPHostType>();
if (h_warning_flag())
error->warning(FLERR,"Dihedral problem",0);
error->warning(FLERR,"Dihedral problem");
if (eflag_global) energy += ev.evdwl;
if (vflag_global) {

2
src/KOKKOS/dihedral_opls_kokkos.cpp
View File

@ -128,7 +128,7 @@ void DihedralOPLSKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
k_warning_flag.template modify<DeviceType>();
k_warning_flag.template sync<LMPHostType>();
if (h_warning_flag())
error->warning(FLERR,"Dihedral problem",0);
error->warning(FLERR,"Dihedral problem");
if (eflag_global) energy += ev.evdwl;
if (vflag_global) {

2
src/KOKKOS/fix_shake_kokkos.cpp
View File

@ -394,7 +394,7 @@ void FixShakeKokkos<DeviceType>::post_force(int vflag)
Kokkos::deep_copy(h_error_flag,d_error_flag);
if (h_error_flag() == 2)
error->warning(FLERR,"Shake determinant < 0.0",0);
error->warning(FLERR,"Shake determinant < 0.0");
else if (h_error_flag() == 3)
error->one(FLERR,"Shake determinant = 0.0");

2
src/KOKKOS/improper_class2_kokkos.cpp
View File

@ -142,7 +142,7 @@ void ImproperClass2Kokkos<DeviceType>::compute(int eflag_in, int vflag_in)
k_warning_flag.template modify<DeviceType>();
k_warning_flag.template sync<LMPHostType>();
if (h_warning_flag())
error->warning(FLERR,"Improper problem",0);
error->warning(FLERR,"Improper problem");
// Angle-Angle energy/force

2
src/KOKKOS/improper_harmonic_kokkos.cpp
View File

@ -130,7 +130,7 @@ void ImproperHarmonicKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
k_warning_flag.template modify<DeviceType>();
k_warning_flag.template sync<LMPHostType>();
if (h_warning_flag())
error->warning(FLERR,"Dihedral problem",0);
error->warning(FLERR,"Dihedral problem");
if (eflag_global) energy += ev.evdwl;
if (vflag_global) {

128
src/KOKKOS/neigh_bond_kokkos.cpp
View File

@ -288,24 +288,16 @@ void NeighBondKokkos<DeviceType>::bond_all()
}
} while (h_fail_flag());
if (nmissing && lostbond == Thermo::ERROR) {
char str[128];
sprintf(str,"Bond atoms missing on proc %d at step " BIGINT_FORMAT,
me,update->ntimestep);
error->one(FLERR,str);
}
if (nmissing && lostbond == Thermo::ERROR)
error->one(FLERR,"Bond atoms missing at step {}", update->ntimestep);
if (neighbor->cluster_check) bond_check();
if (lostbond == Thermo::IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Bond atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
if (all && me ==0)
error->warning(FLERR,"Bond atoms missing at step {}", update->ntimestep);
k_bondlist.modify<DeviceType>();
}
@ -383,24 +375,16 @@ void NeighBondKokkos<DeviceType>::bond_partial()
}
} while (h_fail_flag());
if (nmissing && lostbond == Thermo::ERROR) {
char str[128];
sprintf(str,"Bond atoms missing on proc %d at step " BIGINT_FORMAT,
me,update->ntimestep);
error->one(FLERR,str);
}
if (nmissing && lostbond == Thermo::ERROR)
error->one(FLERR,"Bond atoms missing at step {}", update->ntimestep);
if (neighbor->cluster_check) bond_check();
if (lostbond == Thermo::IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Bond atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
if (all && me ==0)
error->warning(FLERR,"Bond atoms missing at step {}", update->ntimestep);
k_bondlist.modify<DeviceType>();
}
@ -505,24 +489,16 @@ void NeighBondKokkos<DeviceType>::angle_all()
}
} while (h_fail_flag());
if (nmissing && lostbond == Thermo::ERROR) {
char str[128];
sprintf(str,"Angle atoms missing on proc %d at step " BIGINT_FORMAT,
me,update->ntimestep);
error->one(FLERR,str);
}
if (nmissing && lostbond == Thermo::ERROR)
error->one(FLERR,"Angle atoms missing at step {}", update->ntimestep);
if (neighbor->cluster_check) angle_check();
if (lostbond == Thermo::IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Angle atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
if (all && (me == 0))
error->warning(FLERR,"Angle atoms missing at step {}", update->ntimestep);
k_anglelist.modify<DeviceType>();
}
@ -607,24 +583,16 @@ void NeighBondKokkos<DeviceType>::angle_partial()
}
} while (h_fail_flag());
if (nmissing && lostbond == Thermo::ERROR) {
char str[128];
sprintf(str,"Angle atoms missing on proc %d at step " BIGINT_FORMAT,
me,update->ntimestep);
error->one(FLERR,str);
}
if (nmissing && lostbond == Thermo::ERROR)
error->one(FLERR,"Angle atoms missing at step {}", update->ntimestep);
if (neighbor->cluster_check) angle_check();
if (lostbond == Thermo::IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Angle atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
if (all && (me == 0))
error->warning(FLERR,"Angle atoms missing at step {}", update->ntimestep);
k_anglelist.modify<DeviceType>();
}
@ -749,24 +717,16 @@ void NeighBondKokkos<DeviceType>::dihedral_all()
}
} while (h_fail_flag());
if (nmissing && lostbond == Thermo::ERROR) {
char str[128];
sprintf(str,"Dihedral atoms missing on proc %d at step " BIGINT_FORMAT,
me,update->ntimestep);
error->one(FLERR,str);
}
if (nmissing && lostbond == Thermo::ERROR)
error->one(FLERR,"Dihedral atoms missing at step {}", update->ntimestep);
if (neighbor->cluster_check) dihedral_check(neighbor->ndihedrallist,v_dihedrallist);
if (lostbond == Thermo::IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Dihedral atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
if (all && (me == 0))
error->warning(FLERR,"Dihedral atoms missing at step {}", update->ntimestep);
k_dihedrallist.modify<DeviceType>();
}
@ -856,24 +816,16 @@ void NeighBondKokkos<DeviceType>::dihedral_partial()
}
} while (h_fail_flag());
if (nmissing && lostbond == Thermo::ERROR) {
char str[128];
sprintf(str,"Dihedral atoms missing on proc %d at step " BIGINT_FORMAT,
me,update->ntimestep);
error->one(FLERR,str);
}
if (nmissing && lostbond == Thermo::ERROR)
error->one(FLERR,"Dihedral atoms missing at step {}", update->ntimestep);
if (neighbor->cluster_check) dihedral_check(neighbor->ndihedrallist,v_dihedrallist);
if (lostbond == Thermo::IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Dihedral atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
if (all && (me == 0))
error->warning(FLERR,"Dihedral atoms missing at step {}", update->ntimestep);
k_dihedrallist.modify<DeviceType>();
}
@ -1020,24 +972,16 @@ void NeighBondKokkos<DeviceType>::improper_all()
}
} while (h_fail_flag());
if (nmissing && lostbond == Thermo::ERROR) {
char str[128];
sprintf(str,"Improper atoms missing on proc %d at step " BIGINT_FORMAT,
me,update->ntimestep);
error->one(FLERR,str);
}
if (nmissing && lostbond == Thermo::ERROR)
error->one(FLERR,"Improper atoms missing at step {}", update->ntimestep);
if (neighbor->cluster_check) dihedral_check(neighbor->nimproperlist,v_improperlist);
if (lostbond == Thermo::IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Improper atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
if (all && (me == 0))
error->warning(FLERR,"Improper atoms missing at step {}", update->ntimestep);
k_improperlist.modify<DeviceType>();
}
@ -1127,24 +1071,16 @@ void NeighBondKokkos<DeviceType>::improper_partial()
}
} while (h_fail_flag());
if (nmissing && lostbond == Thermo::ERROR) {
char str[128];
sprintf(str,"Improper atoms missing on proc %d at step " BIGINT_FORMAT,
me,update->ntimestep);
error->one(FLERR,str);
}
if (nmissing && lostbond == Thermo::ERROR)
error->one(FLERR,"Improper atoms missing at step {}", update->ntimestep);
if (neighbor->cluster_check) dihedral_check(neighbor->nimproperlist,v_improperlist);
if (lostbond == Thermo::IGNORE) return;
int all;
MPI_Allreduce(&nmissing,&all,1,MPI_INT,MPI_SUM,world);
if (all) {
char str[128];
sprintf(str,
"Improper atoms missing at step " BIGINT_FORMAT,update->ntimestep);
if (me == 0) error->warning(FLERR,str);
}
if (all && (me == 0))
error->warning(FLERR,"Improper atoms missing at step {}", update->ntimestep);
k_improperlist.modify<DeviceType>();
}

7
src/KSPACE/ewald_disp.cpp
View File

@ -214,11 +214,10 @@ void EwaldDisp::init()
error->all(FLERR,"Cannot (yet) use 'electron' units with dipoles");
if (qsqsum == 0.0 && bsbsum == 0.0 && M2 == 0.0)
error->all(FLERR,"Cannot use Ewald/disp solver "
"on system with no charge, dipole, or LJ particles");
error->all(FLERR,"Cannot use Ewald/disp solver on system without "
"charged, dipole, or LJ particles");
if (fabs(qsum) > SMALL && comm->me == 0)
error->warning(FLERR,fmt::format("System is not charge neutral, "
"net charge = {:.8g}",qsum));
error->warning(FLERR,"System is not charge neutral, net charge = {:.8g}",qsum);
if (!function[1] && !function[2]) dispersionflag = 0;
if (!function[3]) dipoleflag = 0;

4
src/KSPACE/msm.cpp
View File

@ -1086,8 +1086,8 @@ void MSM::set_grid_global()
*p_cutoff = cutoff;
if (me == 0)
error->warning(FLERR,fmt::format("Adjusting Coulombic cutoff for "
"MSM, new cutoff = {:.8}", cutoff));
error->warning(FLERR,"Adjusting Coulombic cutoff for "
"MSM, new cutoff = {:.8}", cutoff);
}
if (triclinic == 0) {

4
src/KSPACE/pppm_disp.cpp
View File

@ -1363,8 +1363,8 @@ void PPPMDisp::init_coeffs()
err = bmax/amax;
if (err > 1.0e-4 && comm->me == 0)
error->warning(FLERR,fmt::format("Estimated error in splitting of "
"dispersion coeffs is {}",err));
error->warning(FLERR,"Estimated error in splitting of "
"dispersion coeffs is {}",err);
// set B
B = new double[nsplit*n+nsplit];

7276
src/MANYBODY/pair_bop.cpp
View File

File diff suppressed because it is too large Load Diff

278
src/MANYBODY/pair_bop.h
View File

@ -28,8 +28,10 @@ PairStyle(bop,PairBOP);
#include "pair.h"
namespace LAMMPS_NS {
class TabularFunction;
class PairBOP : public Pair {
public:
PairBOP(class LAMMPS *);
virtual ~PairBOP();
@ -41,118 +43,51 @@ class PairBOP : public Pair {
double memory_usage();
private:
int maxneigh; // maximum size of neighbor list on this processor
int maxneigh3; // maximum size of neighbor list on this processor
int update_list; // check for changing maximum size of neighbor list
int maxbopn; // maximum size of bop neighbor list for allocation
int maxnall; // maximum size of bop neighbor list for allocation
int nr; // increments for the BOP pair potential
int ntheta; // increments for the angle function
int npower; // power of the angular function
int nBOt; // second BO increments
int bop_types; // number of elements in potential
int npairs; // number of element pairs
int bop_step;
int allocate_neigh;
int nb_pi,nb_sg;
int ago1;
int *BOP_index; // index for neighbor list position
int *BOP_total; // index for neighbor list position
int *BOP_index3; // index for neighbor list position
int *BOP_total3; // index for neighbor list position
int *neigh_index; // index for neighbor list position
int *neigh_index3; // index for neighbor list position
int neigh_total; // total number of neighbors stored
int neigh_total3; // total number of neighbors stored
int *cos_index; // index for neighbor cosine if not using on the fly
int *neigh_flag; // index for neighbor cosine if not using on the fly
int *neigh_flag3; // index for neighbor cosine if not using on the fly
int cos_total; // number of cosines stored if not using on the fly
int neigh_ct; // limit for large arrays
struct PairParameters {
double cutB, cutBsq, cutL, cutLsq;
TabularFunction *betaS;
TabularFunction *betaP;
TabularFunction *rep;
TabularFunction *cphi;
TabularFunction *bo;
PairParameters();
~PairParameters();
};
// Parameters variables
struct PairList1 {
double r, dis[3];
double betaS, dBetaS, betaP, dBetaP, rep, dRep;
PairList1(){};
};
int ncutoff,nfunc;
int a_flag;
double *pi_a,*pro_delta,*pi_delta;
double *pi_p,*pi_c,*sigma_r0,*pi_r0,*phi_r0;
double *sigma_rc,*pi_rc,*phi_rc,*r1,*sigma_beta0;
double *pi_beta0,*phi0,*sigma_n,*pi_n,*phi_m;
double *sigma_nc,*pi_nc,*phi_nc;
double *pro,*sigma_delta,*sigma_c,*sigma_a;
double *sigma_f,*sigma_k,*small3;
double small1,small2,small3g,small4,small5,small6,small7;
double which,alpha,alpha1,beta1,gamma1,alpha2,beta2,alpha3;
double beta3,rsmall,rbig,rcore;
char **words;
struct PairList2 {
double r, dis[3];
double rep, dRep;
PairList2(){};
};
double cutmax; // max cutoff for all elements
int otfly; // Defines whether to do on the fly
// calculations of angles and distances
// on the fly will slow down calculations
// but requires less memory on = 1, off=0
struct TripleList {
double G, dG, cosAng, dCosAngi[3], dCosAngj[3], dCosAngk[3];
TripleList(){};
};
// Neigh variables
struct B_SG {
double dAA[3];
double dBB[3];
double dCC[3];
double dDD[3];
double dEE1[3];
double dFF[3];
double dAAC[3];
double dSigB1[3];
double dSigB[3];
int temp;
int i;
int j;
};
double *rcut,*rcut3,*dr,*rdr,*dr3,*rdr3;
double *rcutsq,*rcutsq3;
double **disij,*rij;
double rcutall,rctroot;
// Triple variables
double *cosAng,***dcosAng,***dcAng;
// Double variables
double *betaS,*dBetaS,*betaP;
double *dBetaP,*repul,*dRepul;
// Sigma variables
int *itypeSigBk,nSigBk;
double sigB,sigB_0;
double sigB1;
// Pi variables
int *itypePiBk,nPiBk;
double piB,piB_0;
// Grids1 variables
double **pBetaS,**pBetaS1,**pBetaS2,**pBetaS3;
double **pBetaS4,**pBetaS5,**pBetaS6;
// Grids2 variables
double **pBetaP,**pBetaP1,**pBetaP2,**pBetaP3;
double **pBetaP4,**pBetaP5,**pBetaP6;
// Grids3 variables
double **pRepul,**pRepul1,**pRepul2,**pRepul3;
double **pRepul4,**pRepul5,**pRepul6;
double **pLong,**pLong1,**pLong2,**pLong3;
double **pLong4,**pLong5,**pLong6;
double ****gfunc,****gpara;
double ****gfunc1,****gfunc2,****gfunc3;
double ****gfunc4,****gfunc5,****gfunc6;
double dtheta,rdtheta;
// Grids4 variables
double **FsigBO,**FsigBO1,**FsigBO2,**FsigBO3;
double **FsigBO4,**FsigBO5,**FsigBO6;
double dBO,rdBO;
// End of BOP variables
double **rcmin,**rcmax,**rcmaxp;
struct B_PI{
struct B_PI {
double dAA[3];
double dBB[3];
double dPiB[3];
@ -160,68 +95,77 @@ class PairBOP : public Pair {
int i;
int j;
};
PairParameters *pairParameters;
TabularFunction *tripletParameters;
// Parameters variables
double small1, small2, small3g, small4, small5, small6, small7, *pi_p;
double *sigma_c, *sigma_a, *pi_c, *pi_a, *sigma_delta, *pi_delta;
double *sigma_f, *sigma_k, *small3;
double *pro_delta, *pro;
int bop_types; // number of elments in potential file
int npairs; // number of element pairs
int ntriples; // number of all triples
char **bop_elements; // names of elements in potential file
double bytes;
int otfly; // = 1 faster, more memory, = 0 slower, less memory
PairList1 *pairlist1;
PairList2 *pairlist2;
TripleList *triplelist;
B_SG *bt_sg;
B_PI *bt_pi;
struct B_SG{
double dAA[3];
double dBB[3];
double dCC[3];
double dDD[3];
double dEE[3];
double dEE1[3];
double dFF[3];
double dAAC[3];
double dBBC[3];
double dCCC[3];
double dDDC[3];
double dEEC[3];
double dFFC[3];
double dGGC[3];
double dUT[3];
double dSigB1[3];
double dSigB[3];
int temp;
int i;
int j;
};
B_SG *bt_sg;
int *BOP_index; // index for neighbor list position
int *BOP_total; // index for neighbor list position
int *BOP_index2; // index for neighbor list position
int *BOP_total2; // index for neighbor list position
int *neigh_index; // index for neighbor list position
int *neigh_index2; // index for neighbor list position
int atomlimit; // current size of atom based list
int neighlimit; // current size of neighbor based list
int neighlimit2; // current size of neighbor based list
int neineilimit; // current size of triple based list
int sglimit; // current size of bt_sg
int pilimit; // current size of bt_pi
int *cos_index; // index for neighbor cosine if not using on the fly
double cutmax;
#if defined(LMP_BOP_WRITE_TABLES)
void write_tables(int);
#endif
void setPbetaS();
void setPbetaP();
void setPrepul();
void setSign();
void gneigh();
double sigmaBo(int, int);
void angle(double, double *, double, double *, double &, double *, double *);
double SigmaBo(int, int);
double PiBo(int, int);
void memory_theta_create();
void memory_theta_destroy();
void memory_theta_grow();
void read_table(char *);
void allocate();
void allocate_tables();
void create_pi(int);
void create_sigma(int);
void destroy_pi();
void destroy_sigma();
void grow_pi(int,int);
void grow_sigma(int,int);
void memory_sg(int);
void memory_pi(int);
void initial_sg(int);
void initial_pi(int);
};
}
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
E: Illegal pair_style command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
Self-explanatory.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
Self-explanatory. Check the input script or data file.
E: Pair style BOP requires atom IDs
@ -229,34 +173,20 @@ This is a requirement to use the BOP potential.
E: Pair style BOP requires newton pair on
See the newton command. This is a restriction to use the BOP
potential.
This is a restriction to use the BOP potential.
E: Pair style bop requires comm ghost cutoff at least 3x larger than %g
E: Pair style bop requires a comm ghost cutoff of at least %lf
Use the communicate ghost command to set this. See the pair bop
doc page for more details.
Use the comm_modify cutoff to set this. See the pair bop doc page for
more details.
E: All pair coeffs are not set
All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.
E: Too many atom pairs for pair bop
The number of atomic pairs exceeds the expected number. Check your
atomic structure to ensure that it is realistic.
E: Too many atom triplets for pair bop
The number of three atom groups for angle determinations exceeds the
expected number. Check your atomic structure to ensure that it is
realistic.
Self-explanatory.
E: Cannot open BOP potential file %s
The specified BOP potential file cannot be opened. Check that the
path and name are correct.
Self-explanatory.
E: Incorrect table format check for element types
@ -264,6 +194,12 @@ Self-explanatory.
E: Unsupported BOP potential file format
Self-explanatory.
E: Pair style bop requires system dimension of at least %g
Self-explanatory.
UNDOCUMENTED
*/

4
src/MANYBODY/pair_lcbop.cpp
View File

@ -19,8 +19,6 @@
#include "pair_lcbop.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "force.h"
#include "comm.h"
@ -31,6 +29,8 @@
#include "memory.h"
#include "error.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;

382
src/MANYBODY/pair_polymorphic.cpp
View File

@ -17,6 +17,9 @@
This modifies from pair_tersoff.cpp by Aidan Thompson (SNL)
------------------------------------------------------------------------- */
// uncomment define to enable writing table files for debugging
// #define LMP_POLYMORPHIC_WRITE_TABLES 1
#include "pair_polymorphic.h"
#include "atom.h"
@ -29,6 +32,7 @@
#include "neigh_request.h"
#include "neighbor.h"
#include "potential_file_reader.h"
#include "tabular_function.h"
#include "tokenizer.h"
#include <cmath>
@ -40,6 +44,42 @@ using namespace MathExtra;
#define MAXLINE 1024
#define DELTA 4
/* ---------------------------------------------------------------------- */
PairPolymorphic::PairParameters::PairParameters()
{
cut = 0.0;
cutsq = 0.0;
xi = 0.0;
U = nullptr;
V = nullptr;
W = nullptr;
F = nullptr;
}
PairPolymorphic::PairParameters::~PairParameters()
{
delete U;
delete V;
delete W;
delete F;
}
/* ---------------------------------------------------------------------- */
PairPolymorphic::TripletParameters::TripletParameters()
{
P = nullptr;
G = nullptr;
}
PairPolymorphic::TripletParameters::~TripletParameters()
{
delete P;
delete G;
}
/* ---------------------------------------------------------------------- */
PairPolymorphic::PairPolymorphic(LAMMPS *lmp) : Pair(lmp)
@ -172,7 +212,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
if (eta == 1) {
iparam_ii = elem2param[itype][itype];
PairParameters & p = pairParameters[iparam_ii];
PairParameters &p = pairParameters[iparam_ii];
emb = (p.F)->get_vmax();
}
@ -192,13 +232,13 @@ void PairPolymorphic::compute(int eflag, int vflag)
r0 = sqrt(rsq);
iparam_ij = elem2param[itype][jtype];
PairParameters & p = pairParameters[iparam_ij];
PairParameters &p = pairParameters[iparam_ij];
// do not include the neighbor if get_vmax() <= epsilon because the function is near zero
if (eta == 1) {
if (emb > epsilon) {
iparam_jj = elem2param[jtype][jtype];
PairParameters & q = pairParameters[iparam_jj];
PairParameters &q = pairParameters[iparam_jj];
if (rsq < (q.W)->get_xmaxsq() && (q.W)->get_vmax() > epsilon) {
numneighW = numneighW + 1;
firstneighW[numneighW] = j;
@ -262,7 +302,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
if (emb > epsilon) {
iparam_ii = elem2param[itype][itype];
PairParameters & p = pairParameters[iparam_ii];
PairParameters &p = pairParameters[iparam_ii];
// accumulate bondorder zeta for each i-j interaction via loop over k
@ -273,7 +313,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
ktype = map[type[k]];
iparam_kk = elem2param[ktype][ktype];
PairParameters & q = pairParameters[iparam_kk];
PairParameters &q = pairParameters[iparam_kk];
(q.W)->value(drW[kk],wfac,1,fpair,0);
@ -300,7 +340,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
delr2[2] = -delzW[kk];
iparam_kk = elem2param[ktype][ktype];
PairParameters & q = pairParameters[iparam_kk];
PairParameters &q = pairParameters[iparam_kk];
(q.W)->value(drW[kk],wfac,0,fpair,1);
fpair = -prefactor*fpair/drW[kk];
@ -323,7 +363,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
jtype = map[type[j]];
iparam_ij = elem2param[itype][jtype];
PairParameters & p = pairParameters[iparam_ij];
PairParameters &p = pairParameters[iparam_ij];
delr1[0] = -delxV[jj];
delr1[1] = -delyV[jj];
@ -340,7 +380,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
if (j == k) continue;
ktype = map[type[k]];
iparam_ijk = elem3param[jtype][itype][ktype];
TripletParameters & trip = tripletParameters[iparam_ijk];
TripletParameters &trip = tripletParameters[iparam_ijk];
if ((trip.G)->get_vmax() <= epsilon) continue;
numneighW1 = numneighW1 + 1;
@ -355,7 +395,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
delr1[2]*delr2[2]) / (r1*r2);
iparam_ik = elem2param[itype][ktype];
PairParameters & q = pairParameters[iparam_ik];
PairParameters &q = pairParameters[iparam_ik];
(q.W)->value(r2,wfac,1,fpair,0);
(trip.P)->value(r1-(p.xi)*r2,pfac,1,fpair,0);
@ -389,7 +429,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
k = firstneighW[kk];
ktype = map[type[k]];
iparam_ijk = elem3param[jtype][itype][ktype];
TripletParameters & trip = tripletParameters[iparam_ijk];
TripletParameters &trip = tripletParameters[iparam_ijk];
delr2[0] = -delxW[kk];
delr2[1] = -delyW[kk];
@ -397,7 +437,7 @@ void PairPolymorphic::compute(int eflag, int vflag)
r2 = drW[kk];
iparam_ik = elem2param[itype][ktype];
PairParameters & q = pairParameters[iparam_ik];
PairParameters &q = pairParameters[iparam_ik];
attractive(&p,&q,&trip,prefactor,r1,r2,delr1,delr2,fi,fj,fk);
@ -568,7 +608,7 @@ void PairPolymorphic::read_file(char *file)
tripletParameters = new TripletParameters[ntriple];
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
PairParameters &p = pairParameters[i];
values = reader->next_values(2);
p.cut = values.next_double();
p.cutsq = p.cut*p.cut;
@ -582,6 +622,7 @@ void PairPolymorphic::read_file(char *file)
MPI_Bcast(&nr, 1, MPI_INT, 0, world);
MPI_Bcast(&ng, 1, MPI_INT, 0, world);
MPI_Bcast(&nx, 1, MPI_INT, 0, world);
MPI_Bcast(&eta, 1, MPI_INT, 0, world);
MPI_Bcast(&maxX, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&npair, 1, MPI_INT, 0, world);
@ -602,110 +643,92 @@ void PairPolymorphic::read_file(char *file)
// start reading tabular functions
double * singletable = new double[nr];
for (int i = 0; i < npair; i++) { // U
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
PairParameters &p = pairParameters[i];
if (comm->me == 0) reader->next_dvector(singletable, nr);
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.U = new tabularFunction(nr,0.0,p.cut);
(p.U)->set_values(nr,0.0,p.cut,singletable,epsilon);
p.U = new TabularFunction;
(p.U)->set_values(nr,0.0,p.cut,singletable);
}
for (int i = 0; i < npair; i++) { // V
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
PairParameters &p = pairParameters[i];
if (comm->me == 0) reader->next_dvector(singletable, nr);
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.V = new tabularFunction(nr,0.0,p.cut);
(p.V)->set_values(nr,0.0,p.cut,singletable,epsilon);
p.V = new TabularFunction;
(p.V)->set_values(nr,0.0,p.cut,singletable);
}
for (int i = 0; i < npair; i++) { // W
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
PairParameters &p = pairParameters[i];
if (comm->me == 0) reader->next_dvector(singletable, nr);
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.W = new tabularFunction(nr,0.0,p.cut);
(p.W)->set_values(nr,0.0,p.cut,singletable,epsilon);
p.W = new TabularFunction;
(p.W)->set_values(nr,0.0,p.cut,singletable);
}
cutmax = 0.0;
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
PairParameters &p = pairParameters[i];
if (p.cut > cutmax) cutmax = p.cut;
}
cutmaxsq = cutmax*cutmax;
if (eta != 3) {
for (int j = 0; j < nelements; j++) { // P
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
if (comm->me == 0) reader->next_dvector(singletable, nr);
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
for (int i = 0; i < nelements; i++) {
TripletParameters & p = tripletParameters[i*nelements*nelements+j*nelements+j];
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
TripletParameters &p = tripletParameters[i*nelements*nelements+j*nelements+j];
p.P = new TabularFunction;
(p.P)->set_values(nr,-cutmax,cutmax,singletable);
}
}
for (int j = 0; j < nelements-1; j++) { // P
for (int k = j+1; k < nelements; k++) {
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
for (int k = j+1; k < nelements; k++) {
if (comm->me == 0) reader->next_dvector(singletable, nr);
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
for (int i = 0; i < nelements; i++) {
TripletParameters &p = tripletParameters[i*nelements*nelements+j*nelements+k];
p.P = new TabularFunction;
(p.P)->set_values(nr,-cutmax,cutmax,singletable);
TripletParameters &q = tripletParameters[i*nelements*nelements+k*nelements+j];
q.P = new TabularFunction;
(q.P)->set_values(nr,-cutmax,cutmax,singletable);
}
}
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
for (int i = 0; i < nelements; i++) {
TripletParameters & p = tripletParameters[i*nelements*nelements+j*nelements+k];
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
TripletParameters & q = tripletParameters[i*nelements*nelements+k*nelements+j];
q.P = new tabularFunction(nr,-cutmax,cutmax);
(q.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
}
}
}
}
if (eta == 3) {
for (int i = 0; i < ntriple; i++) { // P
TripletParameters & p = tripletParameters[i];
if (comm->me == 0) {
reader->next_dvector(singletable, nr);
}
TripletParameters &p = tripletParameters[i];
if (comm->me == 0) reader->next_dvector(singletable, nr);
MPI_Bcast(singletable,nr,MPI_DOUBLE,0,world);
p.P = new tabularFunction(nr,-cutmax,cutmax);
(p.P)->set_values(nr,-cutmax,cutmax,singletable,epsilon);
p.P = new TabularFunction;
(p.P)->set_values(nr,-cutmax,cutmax,singletable);
}
}
delete[] singletable;
singletable = new double[ng];
for (int i = 0; i < ntriple; i++) { // G
TripletParameters & p = tripletParameters[i];
if (comm->me == 0) {
reader->next_dvector(singletable, ng);
}
TripletParameters &p = tripletParameters[i];
if (comm->me == 0) reader->next_dvector(singletable, ng);
MPI_Bcast(singletable,ng,MPI_DOUBLE,0,world);
p.G = new tabularFunction(ng,-1.0,1.0);
(p.G)->set_values(ng,-1.0,1.0,singletable,epsilon);
p.G = new TabularFunction;
(p.G)->set_values(ng,-1.0,1.0,singletable);
}
delete[] singletable;
singletable = new double[nx];
for (int i = 0; i < npair; i++) { // F
PairParameters & p = pairParameters[i];
if (comm->me == 0) {
reader->next_dvector(singletable, nx);
}
PairParameters &p = pairParameters[i];
if (comm->me == 0) reader->next_dvector(singletable, nx);
MPI_Bcast(singletable,nx,MPI_DOUBLE,0,world);
p.F = new tabularFunction(nx,0.0,maxX);
(p.F)->set_values(nx,0.0,maxX,singletable,epsilon);
p.F = new TabularFunction;
(p.F)->set_values(nx,0.0,maxX,singletable);
}
delete[] singletable;
if (comm->me == 0) {
delete reader;
}
if (comm->me == 0) delete reader;
// recalculate cutoffs of all params
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
PairParameters &p = pairParameters[i];
p.cut = (p.U)->get_xmax();
if (p.cut < (p.V)->get_xmax()) p.cut = (p.V)->get_xmax();
if (p.cut < (p.W)->get_xmax()) p.cut = (p.W)->get_xmax();
@ -715,7 +738,7 @@ void PairPolymorphic::read_file(char *file)
// set cutmax to max of all params
cutmax = 0.0;
for (int i = 0; i < npair; i++) {
PairParameters & p = pairParameters[i];
PairParameters &p = pairParameters[i];
if (cutmax < p.cut) cutmax = p.cut;
}
cutmaxsq = cutmax*cutmax;
@ -740,28 +763,28 @@ void PairPolymorphic::setup_params()
n++;
}
for (i = 0; i < nelements-1; i++) {
for (j = i+1; j < nelements; j++) {
elem2param[match[i]][match[j]] = n;
elem2param[match[j]][match[i]] = n;
n++;
}
for (j = i+1; j < nelements; j++) {
elem2param[match[i]][match[j]] = n;
elem2param[match[j]][match[i]] = n;
n++;
}
}
// map atom triplet to parameter index
n = 0;
for (i = 0; i < nelements; i++)
for (j = 0; j < nelements; j++)
for (k = 0; k < nelements; k++) {
elem3param[match[i]][match[j]][match[k]] = n;
n++;
}
for (j = 0; j < nelements; j++)
for (k = 0; k < nelements; k++) {
elem3param[match[i]][match[j]][match[k]] = n;
n++;
}
// for debugging, call write_tables() to check the tabular functions
// if (comm->me == 0) {
// write_tables(51);
// }
// error->all(FLERR,"Test potential tables");
// for debugging, call write_tables() to check the tabular functions
#if defined(LMP_POLYMORPHIC_WRITE_TABLES)
if (comm->me == 0) write_tables(51);
error->all(FLERR,"Test potential tables");
#endif
}
/* ----------------------------------------------------------------------
@ -769,10 +792,10 @@ void PairPolymorphic::setup_params()
------------------------------------------------------------------------- */
void PairPolymorphic::attractive(PairParameters *p, PairParameters *q,
TripletParameters *trip,
double prefactor, double rij, double rik,
double *delrij, double *delrik,
double *fi, double *fj, double *fk)
TripletParameters *trip,
double prefactor, double rij, double rik,
double *delrij, double *delrik,
double *fi, double *fj, double *fk)
{
double rij_hat[3],rik_hat[3];
double rijinv,rikinv;
@ -789,11 +812,11 @@ void PairPolymorphic::attractive(PairParameters *p, PairParameters *q,
/* ---------------------------------------------------------------------- */
void PairPolymorphic::ters_zetaterm_d(double prefactor,
double *rij_hat, double rij,
double *rik_hat, double rik,
double *dri, double *drj, double *drk,
PairParameters *p, PairParameters *q,
TripletParameters *trip)
double *rij_hat, double rij,
double *rik_hat, double rik,
double *dri, double *drj, double *drk,
PairParameters *p, PairParameters *q,
TripletParameters *trip)
{
double gijk,gijk_d,ex_delr,ex_delr_d,fc,dfc,cos_theta;
double dcosdri[3],dcosdrj[3],dcosdrk[3];
@ -831,8 +854,8 @@ void PairPolymorphic::ters_zetaterm_d(double prefactor,
/* ---------------------------------------------------------------------- */
void PairPolymorphic::costheta_d(double *rij_hat, double rij,
double *rik_hat, double rik,
double *dri, double *drj, double *drk)
double *rik_hat, double rik,
double *dri, double *drj, double *drk)
{
// first element is devative wrt Ri, second wrt Rj, third wrt Rk
@ -847,107 +870,96 @@ void PairPolymorphic::costheta_d(double *rij_hat, double rij,
}
/* ---------------------------------------------------------------------- */
#if defined(LMP_POLYMORPHIC_WRITE_TABLES)
void PairPolymorphic::write_tables(int npts)
{
char tag[6] = "";
if (comm->me != 0) sprintf(tag,"%d",comm->me);
FILE* fp = nullptr;
double xmin,xmax,x,uf,vf,wf,pf,gf,ff,ufp,vfp,wfp,pfp,gfp,ffp;
char line[MAXLINE];
std::string filename;
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,"_UVW");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem2param[i][j];
PairParameters & pair = pairParameters[iparam_ij];
xmin = (pair.U)->get_xmin();
xmax = (pair.U)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.U)->value(x, uf, 1, ufp, 1);
(pair.V)->value(x, vf, 1, vfp, 1);
(pair.W)->value(x, wf, 1, wfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f %12.4f %12.4f %12.4f %12.4f \n",x,uf,vf,wf,ufp,vfp,wfp);
for (int j = 0; j < nelements; j++) {
filename = fmt::format("{}{}_UVW{}",elements[i],
elements[j],comm->me);
fp = fopen(filename.c_str(), "w");
int iparam_ij = elem2param[i][j];
PairParameters &pair = pairParameters[iparam_ij];
xmin = (pair.U)->get_xmin();
xmax = (pair.U)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.U)->value(x, uf, 1, ufp, 1);
(pair.V)->value(x, vf, 1, vfp, 1);
(pair.W)->value(x, wf, 1, wfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f %12.4f %12.4f %12.4f %12.4f \n",
x,uf,vf,wf,ufp,vfp,wfp);
}
fclose(fp);
}
fclose(fp);
}
}
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,elements[k]);
strcat(line,"_P");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem3param[i][j][k];
TripletParameters & pair = tripletParameters[iparam_ij];
xmin = (pair.P)->get_xmin();
xmax = (pair.P)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int n = 0; n < npts; n++) {
x = xmin + (xmax-xmin) * n / (npts-1);
(pair.P)->value(x, pf, 1, pfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,pf,pfp);
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) {
filename = fmt::format("{}{}{}_P{}",elements[i],elements[j],
elements[k],comm->me);
fp = fopen(filename.c_str(), "w");
int iparam_ij = elem3param[i][j][k];
TripletParameters &pair = tripletParameters[iparam_ij];
xmin = (pair.P)->get_xmin();
xmax = (pair.P)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int n = 0; n < npts; n++) {
x = xmin + (xmax-xmin) * n / (npts-1);
(pair.P)->value(x, pf, 1, pfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,pf,pfp);
}
fclose(fp);
}
}
fclose(fp);
}
}
}
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,elements[k]);
strcat(line,"_G");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem3param[i][j][k];
TripletParameters & pair = tripletParameters[iparam_ij];
xmin = (pair.G)->get_xmin();
xmax = (pair.G)->get_xmax();
for (int n = 0; n < npts; n++) {
x = xmin + (xmax-xmin) * n / (npts-1);
(pair.G)->value(x, gf, 1, gfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,gf,gfp);
for (int j = 0; j < nelements; j++) {
for (int k = 0; k < nelements; k++) {
filename = fmt::format("{}{}{}_G{}",elements[i],elements[j],
elements[k],comm->me);
fp = fopen(filename.c_str(), "w");
int iparam_ij = elem3param[i][j][k];
TripletParameters &pair = tripletParameters[iparam_ij];
xmin = (pair.G)->get_xmin();
xmax = (pair.G)->get_xmax();
for (int n = 0; n < npts; n++) {
x = xmin + (xmax-xmin) * n / (npts-1);
(pair.G)->value(x, gf, 1, gfp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,gf,gfp);
}
fclose(fp);
}
}
fclose(fp);
}
}
}
for (int i = 0; i < nelements; i++) {
for (int j = 0; j < nelements; j++) {
strcpy(line,elements[i]);
strcat(line,elements[j]);
strcat(line,"_F");
strcat(line,tag);
fp = fopen(line, "w");
int iparam_ij = elem2param[i][j];
PairParameters & pair = pairParameters[iparam_ij];
xmin = (pair.F)->get_xmin();
xmax = (pair.F)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.F)->value(x, ff, 1, ffp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,ff,ffp);
for (int j = 0; j < nelements; j++) {
filename = fmt::format("{}{}_F{}",elements[i],
elements[j],comm->me);
fp = fopen(filename.c_str(), "w");
int iparam_ij = elem2param[i][j];
PairParameters &pair = pairParameters[iparam_ij];
xmin = (pair.F)->get_xmin();
xmax = (pair.F)->get_xmax();
double xl = xmax - xmin;
xmin = xmin - 0.5*xl;
xmax = xmax + 0.5*xl;
for (int k = 0; k < npts; k++) {
x = xmin + (xmax-xmin) * k / (npts-1);
(pair.F)->value(x, ff, 1, ffp, 1);
fprintf(fp,"%12.4f %12.4f %12.4f \n",x,ff,ffp);
}
fclose(fp);
}
fclose(fp);
}
}
}
#endif

264
src/MANYBODY/pair_polymorphic.h
View File

@ -21,13 +21,13 @@ PairStyle(polymorphic,PairPolymorphic);
#define LMP_PAIR_POLYMORPHIC_H
#include "pair.h"
#include <cmath>
namespace LAMMPS_NS {
// forward declaration
class TabularFunction;
class PairPolymorphic : public Pair {
public:
public:
PairPolymorphic(class LAMMPS *);
virtual ~PairPolymorphic();
@ -37,255 +37,25 @@ class PairPolymorphic : public Pair {
void init_style();
double init_one(int, int);
protected:
class tabularFunction {
public:
tabularFunction() {
size = 0;
xmin = 0.0;
xmax = 0.0;
xmaxsq = xmax*xmax;
vmax = 0.0;
xs = nullptr;
ys = nullptr;
ys1 = nullptr;
ys2 = nullptr;
ys3 = nullptr;
ys4 = nullptr;
ys5 = nullptr;
ys6 = nullptr;
}
tabularFunction(int n) {
size = n;
xmin = 0.0;
xmax = 0.0;
xmaxsq = xmax*xmax;
xs = new double[n];
ys = new double[n];
ys1 = new double[n];
ys2 = new double[n];
ys3 = new double[n];
ys4 = new double[n];
ys5 = new double[n];
ys6 = new double[n];
}
tabularFunction(int n, double x1, double x2) {
size = n;
xmin = x1;
xmax = x2;
xmaxsq = xmax*xmax;
xs = new double[n];
ys = new double[n];
ys1 = new double[n];
ys2 = new double[n];
ys3 = new double[n];
ys4 = new double[n];
ys5 = new double[n];
ys6 = new double[n];
}
virtual ~tabularFunction() {
delete [] xs;
delete [] ys;
delete [] ys1;
delete [] ys2;
delete [] ys3;
delete [] ys4;
delete [] ys5;
delete [] ys6;
}
void set_xrange(double x1, double x2) {
xmin = x1;
xmax = x2;
xmaxsq = xmax*xmax;
}
void set_values(int n, double x1, double x2, double * values, double epsilon)
{
int shrink = 1;
int ilo,ihi;
double vlo,vhi;
ilo = 0;
ihi = n-1;
for (int i = 0; i < n; i++) {
if (fabs(values[i]) <= epsilon) {
ilo = i;
} else {
break;
}
}
for (int i = n-1; i >= 0; i--) {
if (fabs(values[i]) <= epsilon) {
ihi = i;
} else {
break;
}
}
if (ihi < ilo) ihi = ilo;
vlo = values[ilo];
vhi = values[ilo];
for (int i = ilo; i <= ihi; i++) {
if (vlo > values[i]) vlo = values[i];
if (vhi < values[i]) vhi = values[i];
}
// shrink (remove near zero points) reduces cutoff radius, and therefore computational cost
// do not shrink when x2 < 1.1 (angular function) or x2 > 20.0 (non-radial function)
if (x2 < 1.1 || x2 > 20.0) {
shrink = 0;
}
// do not shrink when when list is abnormally small
if (ihi - ilo < 50) {
shrink = 0;
}
// shrink if it is a constant
if (vhi - vlo <= epsilon) {
// shrink = 1;
}
if (shrink == 0) {
ilo = 0;
ihi = n-1;
}
xmin = x1 + (x2-x1)/(n -1)*ilo;
xmax = xmin + (x2-x1)/(n -1)*(ihi-ilo);
xmaxsq = xmax*xmax;
n = ihi - ilo + 1;
resize(n);
for (int i = ilo; i <= ihi; i++) {
ys[i-ilo] = values[i];
}
initialize();
}
void value(double x, double &y, int ny, double &y1, int ny1)
{
double ps = (x - xmin) * rdx;
int ks = ps + 0.5;
if (ks > size-1) ks = size-1;
if (ks < 0 ) ks = 0;
ps = ps - ks;
if (ny) y = ((ys3[ks]*ps + ys2[ks])*ps + ys1[ks])*ps + ys[ks];
if (ny1) y1 = (ys6[ks]*ps + ys5[ks])*ps + ys4[ks];
}
void print_value()
{
printf("%d %f %f %f \n",size,xmin,xmax,rdx);
printf(" \n");
for (int i = 0; i < size; i++) {
printf("%f %f \n",xs[i],ys[i]);
}
}
double get_xmin() {
return xmin;
}
double get_xmax() {
return xmax;
}
double get_xmaxsq() {
return xmaxsq;
}
double get_rdx() {
return rdx;
}
double get_vmax() {
return vmax;
}
protected:
void resize(int n) {
if (n != size) {
size = n;
delete [] xs;
xs = new double[n];
delete [] ys;
ys = new double[n];
delete [] ys1;
ys1 = new double[n];
delete [] ys2;
ys2 = new double[n];
delete [] ys3;
ys3 = new double[n];
delete [] ys4;
ys4 = new double[n];
delete [] ys5;
ys5 = new double[n];
delete [] ys6;
ys6 = new double[n];
}
}
void initialize() {
int n = size;
rdx = (xmax-xmin)/(n-1.0);
vmax = 0.0;
for (int i = 0; i < n; i++) {
if (fabs(ys[i]) > vmax) vmax = fabs(ys[i]);
}
for (int i = 0; i < n; i++) {
xs[i] = xmin+i*rdx;
}
rdx = 1.0 / rdx;
ys1[0] = ys[1] - ys[0];
ys1[1] = 0.5 * (ys[2] - ys[0]);
ys1[n-2] = 0.5 * (ys[n-1] - ys[n-3]);
ys1[n-1] = ys[n-1] - ys[n-2];
for (int i = 2; i < n-2; i++) {
ys1[i]=((ys[i-2]-ys[i+2])+ 8.0*(ys[i+1]-ys[i-1]))/12.0;
}
for (int i = 0; i < n-1; i++) {
ys2[i]=3.0*(ys[i+1]-ys[i])-2.0*ys1[i]-ys1[i+1];
ys3[i]=ys1[i]+ys1[i+1]-2.0*(ys[i+1]-ys[i]);
}
ys2[n-1]=0.0;
ys3[n-1]=0.0;
for (int i = 0; i < n; i++) {
ys4[i]=ys1[i]*rdx;
ys5[i]=2.0*ys2[i]*rdx;
ys6[i]=3.0*ys3[i]*rdx;
}
}
int size;
double xmin,xmax,xmaxsq,rdx,vmax;
double *ys, *ys1, *ys2, *ys3, *ys4, *ys5, *ys6;
double *xs;
};
protected:
struct PairParameters {
double cut;
double cutsq;
double xi;
class tabularFunction *U;
class tabularFunction *V;
class tabularFunction *W;
class tabularFunction *F;
PairParameters() {
cut = 0.0;
cutsq = 0.0;
xi = 1.0;
U = nullptr;
V = nullptr;
W = nullptr;
F = nullptr;
};
~PairParameters() {
delete U;
delete V;
delete W;
delete F;
}
TabularFunction *U;
TabularFunction *V;
TabularFunction *W;
TabularFunction *F;
PairParameters();
~PairParameters();
};
struct TripletParameters {
class tabularFunction *P;
class tabularFunction *G;
TripletParameters() {
P = nullptr;
G = nullptr;
};
~TripletParameters() {
delete P;
delete G;
}
TabularFunction *P;
TabularFunction *G;
TripletParameters();
~TripletParameters();
};
double epsilon;
@ -311,8 +81,9 @@ class PairPolymorphic : public Pair {
virtual void read_file(char *);
void setup_params();
#if defined(LMP_POLYMORPHIC_WRITE_TABLES)
void write_tables(int);
#endif
void attractive(PairParameters *, PairParameters *, TripletParameters *,
double, double, double, double *, double *, double *,
double *, double *);
@ -323,7 +94,6 @@ class PairPolymorphic : public Pair {
void costheta_d(double *, double, double *, double,
double *, double *, double *);
};
}
#endif

14
src/MC/fix_charge_regulation.cpp
View File

@ -43,6 +43,7 @@
#include <cmath>
#include <cstring>
#include <memory>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -725,7 +726,9 @@ void FixChargeRegulation::forward_ions_multival() {
double energy_before = energy_stored;
double factor = 1;
double dummyp[3];
int mm[salt_charge_ratio + 1];// particle ID array for all ions to be inserted
// particle ID array for all ions to be inserted
auto mm = std::unique_ptr<int[]>(new int[salt_charge_ratio + 1]);
if (salt_charge[0] <= -salt_charge[1]) {
// insert one anion and (salt_charge_ratio) cations
@ -779,9 +782,12 @@ void FixChargeRegulation::backward_ions_multival() {
double energy_before = energy_stored;
double factor = 1;
double dummyp[3]; // dummy particle
int mm[salt_charge_ratio + 1]; // particle ID array for all deleted ions
double qq[salt_charge_ratio + 1]; // charge array for all deleted ions
int mask_tmp[salt_charge_ratio + 1]; // temporary mask array
// particle ID array for all deleted ions
auto mm = std::unique_ptr<int[]>(new int[salt_charge_ratio + 1]);
// charge array for all deleted ions
auto qq = std::unique_ptr<double[]>(new double[salt_charge_ratio + 1]);
// temporary mask array
auto mask_tmp = std::unique_ptr<int[]>(new int[salt_charge_ratio + 1]);
if (salt_charge[0] <= -salt_charge[1]) {
// delete one anion and (salt_charge_ratio) cations

4
src/MC/pair_dsmc.cpp
View File

@ -153,9 +153,9 @@ void PairDSMC::compute(int /*eflag*/, int /*vflag*/)
convert_double_to_equivalent_int(num_of_collisions_double);
if (num_of_collisions > number_of_A)
error->warning(FLERR,"Pair dsmc: num_of_collisions > number_of_A",0);
error->warning(FLERR,"Pair dsmc: num_of_collisions > number_of_A");
if (num_of_collisions > number_of_B)
error->warning(FLERR,"Pair dsmc: num_of_collisions > number_of_B",0);
error->warning(FLERR,"Pair dsmc: num_of_collisions > number_of_B");
// perform collisions on pairs of particles in icell

29
src/MISC/fix_deposit.cpp
View File

@ -288,8 +288,8 @@ void FixDeposit::init()
double separation = MAX(2.0*maxradinsert,maxradall+maxradinsert);
if (sqrt(nearsq) < separation && comm->me == 0)
error->warning(FLERR,fmt::format("Fix deposit near setting < possible "
"overlap separation {}",separation));
error->warning(FLERR,"Fix deposit near setting < possible "
"overlap separation {}",separation);
}
}
@ -316,10 +316,12 @@ void FixDeposit::pre_exchange()
if (next_reneighbor != update->ntimestep) return;
// clear ghost count and any ghost bonus data internal to AtomVec
// clear ghost count (and atom map) and any ghost bonus data
// internal to AtomVec
// same logic as beginning of Comm::exchange()
// do it now b/c inserting atoms will overwrite ghost atoms
if (atom->map_style != Atom::MAP_NONE) atom->map_clear();
atom->nghost = 0;
atom->avec->clear_bonus();
@ -579,14 +581,6 @@ void FixDeposit::pre_exchange()
else if (shakeflag)
fixshake->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
// old code: unsuccessful if no proc performed insertion of an atom
// don't think that check is necessary
// if get this far, should always be successful
// would be hard to undo partial insertion for a molecule
// better to check how many atoms could be inserted (w/out inserting)
// then sum to insure all are inserted, before doing actual insertion
// MPI_Allreduce(&flag,&success,1,MPI_INT,MPI_MAX,world);
success = 1;
break;
}
@ -594,7 +588,7 @@ void FixDeposit::pre_exchange()
// warn if not successful b/c too many attempts
if (!success && comm->me == 0)
error->warning(FLERR,"Particle deposition was unsuccessful",0);
error->warning(FLERR,"Particle deposition was unsuccessful");
// reset global natoms,nbonds,etc
// increment maxtag_all and maxmol_all if necessary
@ -622,10 +616,13 @@ void FixDeposit::pre_exchange()
maxmol_all++;
}
}
if (atom->map_style != Atom::MAP_NONE) {
atom->map_init();
atom->map_set();
}
}
// rebuild atom map
if (atom->map_style != Atom::MAP_NONE) {
if (success) atom->map_init();
atom->map_set();
}
// next timestep to insert

2
src/MISC/fix_ttm.cpp
View File

@ -435,7 +435,7 @@ void FixTTM::end_of_step()
num_inner_timesteps = static_cast<int>(update->dt/inner_dt) + 1;
inner_dt = update->dt/double(num_inner_timesteps);
if (num_inner_timesteps > 1000000)
error->warning(FLERR,"Too many inner timesteps in fix ttm",0);
error->warning(FLERR,"Too many inner timesteps in fix ttm");
}
for (int ith_inner_timestep = 0; ith_inner_timestep < num_inner_timesteps;

42
src/MOLECULE/bond_fene.cpp
View File

@ -14,26 +14,20 @@
#include "bond_fene.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "neighbor.h"
#include "comm.h"
#include "update.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondFENE::BondFENE(LAMMPS *lmp) : Bond(lmp)
{
TWO_1_3 = pow(2.0,(1.0/3.0));
}
using MathConst::MY_CUBEROOT2;
/* ---------------------------------------------------------------------- */
@ -86,8 +80,8 @@ void BondFENE::compute(int eflag, int vflag)
// if r > 2*r0 something serious is wrong, abort
if (rlogarg < 0.1) {
error->warning(FLERR,fmt::format("FENE bond too long: {} {} {} {}",
update->ntimestep,atom->tag[i1],atom->tag[i2],sqrt(rsq)));
error->warning(FLERR,"FENE bond too long: {} {} {} {}",
update->ntimestep,atom->tag[i1],atom->tag[i2],sqrt(rsq));
if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond");
rlogarg = 0.1;
}
@ -96,7 +90,7 @@ void BondFENE::compute(int eflag, int vflag)
// force from LJ term
if (rsq < TWO_1_3*sigma[type]*sigma[type]) {
if (rsq < MY_CUBEROOT2*sigma[type]*sigma[type]) {
sr2 = sigma[type]*sigma[type]/rsq;
sr6 = sr2*sr2*sr2;
fbond += 48.0*epsilon[type]*sr6*(sr6-0.5)/rsq;
@ -106,7 +100,7 @@ void BondFENE::compute(int eflag, int vflag)
if (eflag) {
ebond = -0.5 * k[type]*r0sq*log(rlogarg);
if (rsq < TWO_1_3*sigma[type]*sigma[type])
if (rsq < MY_CUBEROOT2*sigma[type]*sigma[type])
ebond += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
}
@ -252,17 +246,15 @@ double BondFENE::single(int type, double rsq, int /*i*/, int /*j*/,
// if r > 2*r0 something serious is wrong, abort
if (rlogarg < 0.1) {
char str[128];
sprintf(str,"FENE bond too long: " BIGINT_FORMAT " %g",
update->ntimestep,sqrt(rsq));
error->warning(FLERR,str,0);
error->warning(FLERR,"FENE bond too long: {} {:.8}",
update->ntimestep,sqrt(rsq));
if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond");
rlogarg = 0.1;
}
double eng = -0.5 * k[type]*r0sq*log(rlogarg);
fforce = -k[type]/rlogarg;
if (rsq < TWO_1_3*sigma[type]*sigma[type]) {
if (rsq < MY_CUBEROOT2*sigma[type]*sigma[type]) {
double sr2,sr6;
sr2 = sigma[type]*sigma[type]/rsq;
sr6 = sr2*sr2*sr2;

7
src/MOLECULE/bond_fene.h
View File

@ -26,7 +26,7 @@ namespace LAMMPS_NS {
class BondFENE : public Bond {
public:
BondFENE(class LAMMPS *);
BondFENE(class LAMMPS *lmp) : Bond(lmp) {}
virtual ~BondFENE();
virtual void compute(int, int);
virtual void coeff(int, char **);
@ -39,13 +39,12 @@ class BondFENE : public Bond {
virtual void *extract(const char *, int &);
protected:
double TWO_1_3;
double *k,*r0,*epsilon,*sigma;
double *k, *r0, *epsilon, *sigma;
virtual void allocate();
};
}
} // namespace LAMMPS_NS
#endif
#endif

38
src/MOLECULE/bond_fene_expand.cpp
View File

@ -14,24 +14,19 @@
#include "bond_fene_expand.h"
#include <cmath>
#include "atom.h"
#include "neighbor.h"
#include "comm.h"
#include "update.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondFENEExpand::BondFENEExpand(LAMMPS *lmp) : Bond(lmp)
{
TWO_1_3 = pow(2.0,(1.0/3.0));
}
using MathConst::MY_CUBEROOT2;
/* ---------------------------------------------------------------------- */
@ -89,11 +84,8 @@ void BondFENEExpand::compute(int eflag, int vflag)
// if r > 2*r0 something serious is wrong, abort
if (rlogarg < 0.1) {
char str[128];
sprintf(str,"FENE bond too long: " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " %g",
update->ntimestep,atom->tag[i1],atom->tag[i2],sqrt(rsq));
error->warning(FLERR,str,0);
error->warning(FLERR,"FENE bond too long: {} {} {} {:.8}",
update->ntimestep,atom->tag[i1],atom->tag[i2],sqrt(rsq));
if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond");
rlogarg = 0.1;
}
@ -102,7 +94,7 @@ void BondFENEExpand::compute(int eflag, int vflag)
// force from LJ term
if (rshiftsq < TWO_1_3*sigma[type]*sigma[type]) {
if (rshiftsq < MY_CUBEROOT2*sigma[type]*sigma[type]) {
sr2 = sigma[type]*sigma[type]/rshiftsq;
sr6 = sr2*sr2*sr2;
fbond += 48.0*epsilon[type]*sr6*(sr6-0.5)/rshift/r;
@ -112,7 +104,7 @@ void BondFENEExpand::compute(int eflag, int vflag)
if (eflag) {
ebond = -0.5 * k[type]*r0sq*log(rlogarg);
if (rshiftsq < TWO_1_3*sigma[type]*sigma[type])
if (rshiftsq < MY_CUBEROOT2*sigma[type]*sigma[type])
ebond += 4.0*epsilon[type]*sr6*(sr6-1.0) + epsilon[type];
}
@ -267,17 +259,15 @@ double BondFENEExpand::single(int type, double rsq, int /*i*/, int /*j*/,
// if r > 2*r0 something serious is wrong, abort
if (rlogarg < 0.1) {
char str[128];
sprintf(str,"FENE bond too long: " BIGINT_FORMAT " %g",
update->ntimestep,sqrt(rsq));
error->warning(FLERR,str,0);
error->warning(FLERR,"FENE bond too long: {} {:.8}",
update->ntimestep,sqrt(rsq));
if (rlogarg <= -3.0) error->one(FLERR,"Bad FENE bond");
rlogarg = 0.1;
}
double eng = -0.5 * k[type]*r0sq*log(rlogarg);
fforce = -k[type]*rshift/rlogarg/r;
if (rshiftsq < TWO_1_3*sigma[type]*sigma[type]) {
if (rshiftsq < MY_CUBEROOT2*sigma[type]*sigma[type]) {
double sr2,sr6;
sr2 = sigma[type]*sigma[type]/rshiftsq;
sr6 = sr2*sr2*sr2;

7
src/MOLECULE/bond_fene_expand.h
View File

@ -26,7 +26,7 @@ namespace LAMMPS_NS {
class BondFENEExpand : public Bond {
public:
BondFENEExpand(class LAMMPS *);
BondFENEExpand(class LAMMPS *lmp) : Bond(lmp) {}
virtual ~BondFENEExpand();
virtual void compute(int, int);
void coeff(int, char **);
@ -38,13 +38,12 @@ class BondFENEExpand : public Bond {
double single(int, double, int, int, double &);
protected:
double TWO_1_3;
double *k,*r0,*epsilon,*sigma,*shift;
double *k, *r0, *epsilon, *sigma, *shift;
void allocate();
};
}
} // namespace LAMMPS_NS
#endif
#endif

36
src/MOLECULE/bond_table.cpp
View File

@ -18,19 +18,17 @@
#include "bond_table.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "neighbor.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "error.h"
#include "tokenizer.h"
#include "neighbor.h"
#include "table_file_reader.h"
#include "tokenizer.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
@ -379,15 +377,15 @@ void BondTable::read_table(Table *tb, char *file, char *keyword)
}
if (ferror) {
error->warning(FLERR, fmt::format("{} of {} force values in table are inconsistent with -dE/dr.\n"
" Should only be flagged at inflection points",ferror,tb->ninput));
error->warning(FLERR, "{} of {} force values in table are inconsistent with -dE/dr.\n"
"WARNING: Should only be flagged at inflection points",ferror,tb->ninput);
}
// warn if data was read incompletely, e.g. columns were missing
if (cerror) {
error->warning(FLERR, fmt::format("{} of {} lines in table were incomplete or could not be"
" parsed completely",cerror,tb->ninput));
error->warning(FLERR, "{} of {} lines in table were incomplete or could not be"
" parsed completely",cerror,tb->ninput);
}
}
@ -597,18 +595,14 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f)
}
double fraction,a,b;
char estr[128];
const Table *tb = &tables[tabindex[type]];
const int itable = static_cast<int> ((x - tb->lo) * tb->invdelta);
if (itable < 0) {
sprintf(estr,"Bond length < table inner cutoff: "
"type %d length %g",type,x);
error->one(FLERR,estr);
} else if (itable >= tablength) {
sprintf(estr,"Bond length > table outer cutoff: "
"type %d length %g",type,x);
error->one(FLERR,estr);
}
if (itable < 0)
error->one(FLERR,"Bond length < table inner cutoff: "
"type {} length {:.8}",type,x);
else if (itable >= tablength)
error->one(FLERR,"Bond length > table outer cutoff: "
"type {} length {:.8}",type,x);
if (tabstyle == LINEAR) {
fraction = (x - tb->r[itable]) * tb->invdelta;

23
src/MOLECULE/dihedral_charmm.cpp
View File

@ -144,27 +144,8 @@ void DihedralCharmm::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me;
MPI_Comm_rank(world,&me);
if (screen) {
char str[128];
sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;

23
src/MOLECULE/dihedral_charmmfsw.cpp
View File

@ -147,27 +147,8 @@ void DihedralCharmmfsw::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me;
MPI_Comm_rank(world,&me);
if (screen) {
char str[128];
sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;

23
src/MOLECULE/dihedral_harmonic.cpp
View File

@ -130,27 +130,8 @@ void DihedralHarmonic::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me;
MPI_Comm_rank(world,&me);
if (screen) {
char str[128];
sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;

23
src/MOLECULE/dihedral_helix.cpp
View File

@ -161,27 +161,8 @@ void DihedralHelix::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me = comm->me;
if (screen) {
char str[128];
sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;

23
src/MOLECULE/dihedral_multi_harmonic.cpp
View File

@ -153,27 +153,8 @@ void DihedralMultiHarmonic::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me;
MPI_Comm_rank(world,&me);
if (screen) {
char str[128];
sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;

23
src/MOLECULE/dihedral_opls.cpp
View File

@ -159,27 +159,8 @@ void DihedralOPLS::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me;
MPI_Comm_rank(world,&me);
if (screen) {
char str[128];
sprintf(str,"Dihedral problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;

23
src/MOLECULE/improper_cvff.cpp
View File

@ -144,27 +144,8 @@ void ImproperCvff::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me;
MPI_Comm_rank(world,&me);
if (screen) {
char str[128];
sprintf(str,"Improper problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;

23
src/MOLECULE/improper_harmonic.cpp
View File

@ -118,27 +118,8 @@ void ImproperHarmonic::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me;
MPI_Comm_rank(world,&me);
if (screen) {
char str[128];
sprintf(str,"Improper problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;

23
src/MOLECULE/improper_umbrella.cpp
View File

@ -125,27 +125,8 @@ void ImproperUmbrella::compute(int eflag, int vflag)
// error check
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) {
int me;
MPI_Comm_rank(world,&me);
if (screen) {
char str[128];
sprintf(str,"Improper problem: %d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1][0],x[i1][1],x[i1][2]);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2][0],x[i2][1],x[i2][2]);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3][0],x[i3][1],x[i3][2]);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4][0],x[i4][1],x[i4][2]);
}
}
if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE))
problem(FLERR, i1, i2, i3, i4);
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;

45
src/PLUGIN/plugin.cpp
View File

@ -172,9 +172,8 @@ namespace LAMMPS_NS
auto pair_map = lmp->force->pair_map;
if (pair_map->find(plugin->name) != pair_map->end()) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,fmt::format("Overriding built-in pair "
"style {} from plugin",
plugin->name));
lmp->error->warning(FLERR,"Overriding built-in pair "
"style {} from plugin",plugin->name);
}
(*pair_map)[plugin->name] = (Force::PairCreator)plugin->creator.v1;
@ -182,9 +181,8 @@ namespace LAMMPS_NS
auto bond_map = lmp->force->bond_map;
if (bond_map->find(plugin->name) != bond_map->end()) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,fmt::format("Overriding built-in bond "
"style {} from plugin",
plugin->name));
lmp->error->warning(FLERR,"Overriding built-in bond "
"style {} from plugin",plugin->name);
}
(*bond_map)[plugin->name] = (Force::BondCreator)plugin->creator.v1;
@ -192,9 +190,8 @@ namespace LAMMPS_NS
auto angle_map = lmp->force->angle_map;
if (angle_map->find(plugin->name) != angle_map->end()) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,fmt::format("Overriding built-in angle "
"style {} from plugin",
plugin->name));
lmp->error->warning(FLERR,"Overriding built-in angle "
"style {} from plugin",plugin->name);
}
(*angle_map)[plugin->name] = (Force::AngleCreator)plugin->creator.v1;
@ -202,9 +199,8 @@ namespace LAMMPS_NS
auto dihedral_map = lmp->force->dihedral_map;
if (dihedral_map->find(plugin->name) != dihedral_map->end()) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,fmt::format("Overriding built-in dihedral "
"style {} from plugin",
plugin->name));
lmp->error->warning(FLERR,"Overriding built-in dihedral "
"style {} from plugin",plugin->name);
}
(*dihedral_map)[plugin->name] = (Force::DihedralCreator)plugin->creator.v1;
@ -212,9 +208,8 @@ namespace LAMMPS_NS
auto improper_map = lmp->force->improper_map;
if (improper_map->find(plugin->name) != improper_map->end()) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,fmt::format("Overriding built-in improper "
"style {} from plugin",
plugin->name));
lmp->error->warning(FLERR,"Overriding built-in improper "
"style {} from plugin",plugin->name);
}
(*improper_map)[plugin->name] = (Force::ImproperCreator)plugin->creator.v1;
@ -222,9 +217,8 @@ namespace LAMMPS_NS
auto compute_map = lmp->modify->compute_map;
if (compute_map->find(plugin->name) != compute_map->end()) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,fmt::format("Overriding built-in compute "
"style {} from plugin",
plugin->name));
lmp->error->warning(FLERR,"Overriding built-in compute "
"style {} from plugin",plugin->name);
}
(*compute_map)[plugin->name] = (Modify::ComputeCreator)plugin->creator.v2;
@ -232,9 +226,8 @@ namespace LAMMPS_NS
auto fix_map = lmp->modify->fix_map;
if (fix_map->find(plugin->name) != fix_map->end()) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,fmt::format("Overriding built-in fix "
"style {} from plugin",
plugin->name));
lmp->error->warning(FLERR,"Overriding built-in fix "
"style {} from plugin",plugin->name);
}
(*fix_map)[plugin->name] = (Modify::FixCreator)plugin->creator.v2;
@ -242,9 +235,8 @@ namespace LAMMPS_NS
auto region_map = lmp->domain->region_map;
if (region_map->find(plugin->name) != region_map->end()) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,fmt::format("Overriding built-in region "
"style {} from plugin",
plugin->name));
lmp->error->warning(FLERR,"Overriding built-in region "
"style {} from plugin",plugin->name);
}
(*region_map)[plugin->name] = (Domain::RegionCreator)plugin->creator.v2;
@ -252,9 +244,8 @@ namespace LAMMPS_NS
auto command_map = lmp->input->command_map;
if (command_map->find(plugin->name) != command_map->end()) {
if (lmp->comm->me == 0)
lmp->error->warning(FLERR,fmt::format("Overriding built-in command "
"style {} from plugin",
plugin->name));
lmp->error->warning(FLERR,"Overriding built-in command "
"style {} from plugin",plugin->name);
}
(*command_map)[plugin->name] = (Input::CommandCreator)plugin->creator.v1;

8
src/POEMS/fix_poems.cpp
View File

@ -942,11 +942,9 @@ void FixPOEMS::readfile(char *file)
if (me == 0) {
fp = fopen(file,"r");
if (fp == nullptr) {
char str[128];
snprintf(str,128,"Cannot open fix poems file %s",file);
error->one(FLERR,str);
}
if (fp == nullptr)
error->one(FLERR,"Cannot open fix poems file {}: {}",
file, utils::getsyserror());
}
nbody = 0;

7
src/QEQ/fix_qeq.cpp
View File

@ -388,10 +388,9 @@ int FixQEq::CG( double *b, double *x )
}
if ((comm->me == 0) && (loop >= maxiter))
error->warning(FLERR,fmt::format("Fix qeq CG convergence failed ({}) "
"after {} iterations at step {}",
sqrt(sig_new)/b_norm,loop,
update->ntimestep));
error->warning(FLERR,"Fix qeq CG convergence failed ({}) after {} "
"iterations at step {}",sqrt(sig_new)/b_norm,loop,
update->ntimestep);
return loop;
}

4
src/QEQ/fix_qeq_dynamic.cpp
View File

@ -156,8 +156,8 @@ void FixQEqDynamic::pre_force(int /*vflag*/)
}
if ((comm->me == 0) && (iloop >= maxiter))
error->warning(FLERR,fmt::format("Charges did not converge at step "
"{}: {}",update->ntimestep,enegchk));
error->warning(FLERR,"Charges did not converge at step {}: {}",
update->ntimestep,enegchk);
if (force->kspace) force->kspace->qsum_qsq();
}

4
src/QEQ/fix_qeq_fire.cpp
View File

@ -216,8 +216,8 @@ void FixQEqFire::pre_force(int /*vflag*/)
}
if ((comm->me == 0) && (iloop >= maxiter))
error->warning(FLERR,fmt::format("Charges did not converge at step "
"{}: {}",update->ntimestep,enegchk));
error->warning(FLERR,"Charges did not converge at step {}: {}",
update->ntimestep,enegchk);
if (force->kspace) force->kspace->qsum_qsq();
}

16
src/REPLICA/neb.cpp
View File

@ -60,6 +60,7 @@ NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
n1steps = n1steps_in;
n2steps = n2steps_in;
nevery = nevery_in;
verbose = false;
// replica info
@ -548,13 +549,11 @@ void NEB::open(char *file)
if (!compressed) fp = fopen(file,"r");
else {
#ifdef LAMMPS_GZIP
char gunzip[128];
snprintf(gunzip,128,"gzip -c -d %s",file);
auto gunzip = std::string("gzip -c -d ") + file;
#ifdef _WIN32
fp = _popen(gunzip,"rb");
fp = _popen(gunzip.c_str(),"rb");
#else
fp = popen(gunzip,"r");
fp = popen(gunzip.c_str(),"r");
#endif
#else
@ -562,11 +561,8 @@ void NEB::open(char *file)
#endif
}
if (fp == nullptr) {
char str[128];
snprintf(str,128,"Cannot open file %s",file);
error->one(FLERR,str);
}
if (fp == nullptr)
error->one(FLERR,"Cannot open file {}: {}",file,utils::getsyserror());
}
/* ----------------------------------------------------------------------

5
src/REPLICA/tad.cpp
View File

@ -44,7 +44,10 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
TAD::TAD(LAMMPS *lmp) : Command(lmp) {}
TAD::TAD(LAMMPS *lmp) : Command(lmp)
{
deltconf = deltstop = deltfirst = 0.0;
}
/* ---------------------------------------------------------------------- */

4
src/RIGID/fix_rigid.cpp
View File

@ -698,8 +698,8 @@ void FixRigid::init()
if (modify->fix[i]->rigid_flag) rflag = 1;
if (rflag && (modify->fmask[i] & POST_FORCE) &&
!modify->fix[i]->rigid_flag)
error->warning(FLERR,fmt::format("Fix {} alters forces after fix "
"rigid",modify->fix[i]->id));
error->warning(FLERR,"Fix {} alters forces after fix rigid",
modify->fix[i]->id);
}
}

4
src/RIGID/fix_rigid_small.cpp
View File

@ -549,8 +549,8 @@ void FixRigidSmall::init()
if (modify->fix[i]->rigid_flag) rflag = 1;
if (rflag && (modify->fmask[i] & POST_FORCE) &&
!modify->fix[i]->rigid_flag)
error->warning(FLERR,fmt::format("Fix {} alters forces after fix "
"rigid",modify->fix[i]->id));
error->warning(FLERR,"Fix {} alters forces after fix rigid",
modify->fix[i]->id);
}
}

2
src/RIGID/fix_shake.cpp
View File

@ -1717,7 +1717,7 @@ void FixShake::shake(int m)
double determ = b*b - 4.0*a*c;
if (determ < 0.0) {
error->warning(FLERR,"Shake determinant < 0.0",0);
error->warning(FLERR,"Shake determinant < 0.0");
determ = 0.0;
}

3
src/SNAP/pair_snap.cpp
View File

@ -505,7 +505,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
// initialize checklist for all required nelements
int elementflags[nelements];
int *elementflags = new int[nelements];
for (int jelem = 0; jelem < nelements; jelem++)
elementflags[jelem] = 0;
@ -602,6 +602,7 @@ void PairSNAP::read_files(char *coefffilename, char *paramfilename)
error->all(FLERR,"Element {} not found in SNAP coefficient "
"file", elements[jelem]);
}
delete[] elementflags;
// set flags for required keywords

15
src/SPIN/neb_spin.cpp
View File

@ -693,13 +693,11 @@ void NEBSpin::open(char *file)
if (!compressed) fp = fopen(file,"r");
else {
#ifdef LAMMPS_GZIP
char gunzip[128];
snprintf(gunzip,128,"gzip -c -d %s",file);
auto gunzip = std::string("gzip -c -d ") + file;
#ifdef _WIN32
fp = _popen(gunzip,"rb");
fp = _popen(gunzip.c_str(),"rb");
#else
fp = popen(gunzip,"r");
fp = popen(gunzip.c_str(),"r");
#endif
#else
@ -707,11 +705,8 @@ void NEBSpin::open(char *file)
#endif
}
if (fp == nullptr) {
char str[128];
snprintf(str,128,"Cannot open file %s",file);
error->one(FLERR,str);
}
if (fp == nullptr)
error->one(FLERR,"Cannot open file {}: {}",file,utils::getsyserror());
}
/* ----------------------------------------------------------------------

49
src/USER-BOCS/fix_bocs.cpp
View File

@ -630,7 +630,7 @@ int FixBocs::read_F_table( char *filename, int p_basis_type )
// NB: LAMMPS coding guidelines prefer cstdio so we are intentionally
// foregoing reading with getline
if (comm->me == 0) {
error->message(FLERR, fmt::format("INFO: About to read data file: {}", filename));
error->message(FLERR, "INFO: About to read data file: {}", filename);
}
// Data file lines hold two floating point numbers.
@ -646,7 +646,7 @@ int FixBocs::read_F_table( char *filename, int p_basis_type )
numEntries = inputLines.size();
if (comm->me == 0) {
error->message(FLERR, fmt::format("INFO: Read {} lines from file", numEntries));
error->message(FLERR, "INFO: Read {} lines from file", numEntries);
}
@ -671,34 +671,23 @@ int FixBocs::read_F_table( char *filename, int p_basis_type )
test_sscanf = sscanf(inputLines.at(i).c_str()," %f , %f ",&f1, &f2);
if (test_sscanf == 2)
{
//if (comm->me == 0) {
// error->message(FLERR, fmt::format("INFO: f1 = {}, f2 = {}", f1, f2));
//}
data[VOLUME][i] = (double)f1;
data[PRESSURE_CORRECTION][i] = (double)f2;
if (i == 1)
{
// second entry is used to compute the validation interval used below
stdVolumeInterval = data[VOLUME][i] - data[VOLUME][i-1];
//if (comm->me == 0) {
// error->message(FLERR, fmt::format("INFO: standard volume interval computed: {}", stdVolumeInterval));
//}
}
else if (i > 1)
{
// after second entry, all intervals are validated
currVolumeInterval = data[VOLUME][i] - data[VOLUME][i-1];
//if (comm->me == 0) {
// error->message(FLERR, fmt::format("INFO: current volume interval: {}", currVolumeInterval));
//}
if (fabs(currVolumeInterval - stdVolumeInterval) > volumeIntervalTolerance) {
if (comm->me == 0) {
message = fmt::format("Bad volume interval. Spline analysis requires uniform"
" volume distribution, found inconsistent volume"
" differential, line {} of file {}\n\tline: {}",
lineNum, filename, inputLines.at(i));
error->warning(FLERR, message);
}
if (comm->me == 0)
error->warning(FLERR,"Bad volume interval. Spline analysis requires uniform"
" volume distribution, found inconsistent volume"
" differential, line {} of file {}\nWARNING:\tline: {}",
lineNum, filename, inputLines.at(i));
badInput = true;
numBadVolumeIntervals++;
}
@ -707,12 +696,10 @@ int FixBocs::read_F_table( char *filename, int p_basis_type )
}
else
{
if (comm->me == 0) {
message = fmt::format("Bad input format: did not find 2 comma separated numeric"
" values in line {} of file {}\n\tline: {}",
lineNum, filename, inputLines.at(i));
error->warning(FLERR, message);
}
if (comm->me == 0)
error->warning(FLERR,"Bad input format: did not find 2 comma separated numeric"
" values in line {} of file {}\nWARNING:\tline: {}",
lineNum, filename, inputLines.at(i));
badInput = true;
}
if (badInput)
@ -722,7 +709,7 @@ int FixBocs::read_F_table( char *filename, int p_basis_type )
}
if (numBadVolumeIntervals > 0 && comm->me == 0) {
error->message(FLERR, fmt::format("INFO: total number bad volume intervals = {}", numBadVolumeIntervals));
error->message(FLERR, "INFO: total number bad volume intervals = {}", numBadVolumeIntervals);
}
}
else {
@ -730,7 +717,7 @@ int FixBocs::read_F_table( char *filename, int p_basis_type )
}
if (badInput && comm->me == 0) {
error->warning(FLERR,fmt::format("Bad volume / pressure-correction data: {}\nSee details above", filename));
error->warning(FLERR,"Bad volume / pressure-correction data: {}\nSee details above", filename);
}
if (p_basis_type == BASIS_LINEAR_SPLINE)
@ -753,9 +740,6 @@ int FixBocs::read_F_table( char *filename, int p_basis_type )
}
int FixBocs::build_linear_splines(double **data) {
//if (comm->me == 0) {
//error->message(FLERR, fmt::format("INFO: entering build_linear_splines, spline_length = {}", spline_length));
//}
splines = (double **) calloc(NUM_LINEAR_SPLINE_COLUMNS,sizeof(double *));
splines[VOLUME] = (double *) calloc(spline_length,sizeof(double));
splines[PRESSURE_CORRECTION] = (double *) calloc(spline_length,sizeof(double));
@ -767,7 +751,7 @@ int FixBocs::build_linear_splines(double **data) {
}
if (comm->me == 0) {
error->message(FLERR, fmt::format("INFO: leaving build_linear_splines, spline_length = {}", spline_length));
error->message(FLERR, "INFO: leaving build_linear_splines, spline_length = {}", spline_length);
}
return spline_length;
@ -775,9 +759,6 @@ int FixBocs::build_linear_splines(double **data) {
int FixBocs::build_cubic_splines(double **data)
{
//if (comm->me == 0) {
//error->message(FLERR, fmt::format("INFO: entering build_cubic_splines, spline_length = {}", spline_length));
//}
int n = spline_length;
double *a, *b, *d, *h, *alpha, *c, *l, *mu, *z;
// 2020-07-17 ag:
@ -868,7 +849,7 @@ int FixBocs::build_cubic_splines(double **data)
memory->destroy(z);
if (comm->me == 0) {
error->message(FLERR, fmt::format("INFO: leaving build_cubic_splines, numSplines = {}", numSplines));
error->message(FLERR, "INFO: leaving build_cubic_splines, numSplines = {}", numSplines);
}
// Tell the caller how many splines we created

4
src/USER-CGDNA/bond_oxdna_fene.cpp
View File

@ -222,7 +222,7 @@ void BondOxdnaFene::compute(int eflag, int vflag)
sprintf(str,"FENE bond too long: " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " %g",
update->ntimestep,atom->tag[a],atom->tag[b],r);
error->warning(FLERR,str,0);
error->warning(FLERR,str);
rlogarg = 0.1;
}
@ -401,7 +401,7 @@ double BondOxdnaFene::single(int type, double rsq, int /*i*/, int /*j*/,
char str[128];
sprintf(str,"FENE bond too long: " BIGINT_FORMAT " %g",
update->ntimestep,sqrt(rsq));
error->warning(FLERR,str,0);
error->warning(FLERR,str);
rlogarg = 0.1;
}

13
src/USER-INTEL/bond_fene_intel.cpp
View File

@ -23,6 +23,7 @@
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
@ -33,6 +34,7 @@
#include "omp_compat.h"
using namespace LAMMPS_NS;
using MathConst::MY_CUBEROOT2;
typedef struct { int a,b,t; } int3_t;
@ -184,11 +186,8 @@ void BondFENEIntel::eval(const int vflag,
// if r > 2*r0 something serious is wrong, abort
if (rlogarg < (flt_t)0.1) {
char str[128];
sprintf(str,"FENE bond too long: " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " %g",
update->ntimestep,atom->tag[i1],atom->tag[i2],sqrt(rsq));
error->warning(FLERR,str,0);
error->warning(FLERR,"FENE bond too long: {} {} {} {:.8}",
update->ntimestep,atom->tag[i1],atom->tag[i2],sqrt(rsq));
if (rlogarg <= (flt_t)-3.0) error->one(FLERR,"Bad FENE bond");
rlogarg = (flt_t)0.1;
}
@ -198,7 +197,7 @@ void BondFENEIntel::eval(const int vflag,
// force from LJ term
flt_t sr2,sr6;
if (rsq < (flt_t)TWO_1_3*sigmasq) {
if (rsq < (flt_t)MY_CUBEROOT2*sigmasq) {
sr2 = sigmasq * irsq;
sr6 = sr2 * sr2 * sr2;
fbond += (flt_t)48.0 * epsilon * sr6 * (sr6 - (flt_t)0.5) * irsq;
@ -209,7 +208,7 @@ void BondFENEIntel::eval(const int vflag,
flt_t ebond;
if (EFLAG) {
ebond = (flt_t)-0.5 * k / ir0sq * log(rlogarg);
if (rsq < (flt_t)TWO_1_3 * sigmasq)
if (rsq < (flt_t)MY_CUBEROOT2 * sigmasq)
ebond += (flt_t)4.0 * epsilon * sr6 * (sr6 - (flt_t)1.0) + epsilon;
}

23
src/USER-INTEL/dihedral_charmm_intel.cpp
View File

@ -246,27 +246,8 @@ void DihedralCharmmIntel::eval(const int vflag,
// error check
#ifndef LMP_SIMD_COMPILER_TEST
if (c > PTOLERANCE || c < MTOLERANCE) {
int me = comm->me;
if (screen) {
char str[128];
sprintf(str,"Dihedral problem: %d/%d " BIGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT " "
TAGINT_FORMAT " " TAGINT_FORMAT,
me,tid,update->ntimestep,
atom->tag[i1],atom->tag[i2],atom->tag[i3],atom->tag[i4]);
error->warning(FLERR,str,0);
fprintf(screen," 1st atom: %d %g %g %g\n",
me,x[i1].x,x[i1].y,x[i1].z);
fprintf(screen," 2nd atom: %d %g %g %g\n",
me,x[i2].x,x[i2].y,x[i2].z);
fprintf(screen," 3rd atom: %d %g %g %g\n",
me,x[i3].x,x[i3].y,x[i3].z);
fprintf(screen," 4th atom: %d %g %g %g\n",
me,x[i4].x,x[i4].y,x[i4].z);
}
}
if (c > PTOLERANCE || c < MTOLERANCE)
problem(FLERR, i1, i2, i3, i4);
#endif
if (c > (flt_t)1.0) c = (flt_t)1.0;

Some files were not shown because too many files have changed in this diff Show More