diff --git a/cmake/Modules/CodeCoverage.cmake b/cmake/Modules/CodeCoverage.cmake
index d018db43d9..3b323e37ff 100644
--- a/cmake/Modules/CodeCoverage.cmake
+++ b/cmake/Modules/CodeCoverage.cmake
@@ -15,14 +15,35 @@ if(ENABLE_COVERAGE)
gen_coverage_xml
COMMAND ${GCOVR_BINARY} -s -x -r ${ABSOLUTE_LAMMPS_SOURCE_DIR} --object-directory=${CMAKE_BINARY_DIR} -o coverage.xml
WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
- COMMENT "Generating XML Coverage Report..."
+ COMMENT "Generating XML coverage report..."
)
+ set(COVERAGE_HTML_DIR ${CMAKE_BINARY_DIR}/coverage_html)
+
+ add_custom_target(coverage_html_folder
+ COMMAND ${CMAKE_COMMAND} -E make_directory ${COVERAGE_HTML_DIR})
+
add_custom_target(
gen_coverage_html
- COMMAND ${GCOVR_BINARY} -s --html --html-details -r ${ABSOLUTE_LAMMPS_SOURCE_DIR} --object-directory=${CMAKE_BINARY_DIR} -o coverage.html
+ COMMAND ${GCOVR_BINARY} -s --html --html-details -r ${ABSOLUTE_LAMMPS_SOURCE_DIR} --object-directory=${CMAKE_BINARY_DIR} -o ${COVERAGE_HTML_DIR}/index.html
WORKING_DIRECTORY ${CMAKE_BINARY_DIR}
- COMMENT "Generating HTML Coverage Report..."
+ COMMENT "Generating HTML coverage report..."
)
+ add_dependencies(gen_coverage_html coverage_html_folder)
+
+ add_custom_target(clean_coverage_html
+ ${CMAKE_COMMAND} -E remove_directory ${COVERAGE_HTML_DIR}
+ COMMENT "Deleting HTML coverage report..."
+ )
+
+ add_custom_target(reset_coverage
+ ${CMAKE_COMMAND} -E remove -f */*.gcda */*/*.gcda */*/*/*.gcda
+ */*/*/*/*.gcda */*/*/*/*/*.gcda */*/*/*/*/*/*.gcda
+ */*/*/*/*/*/*/*.gcda */*/*/*/*/*/*/*/*.gcda
+ */*/*/*/*/*/*/*/*/*.gcda */*/*/*/*/*/*/*/*/*/*.gcda
+ WORKIND_DIRECTORY ${CMAKE_BINARY_DIR}
+ COMMENT "Deleting coverage data files..."
+ )
+ add_dependencies(reset_coverage clean_coverage_html)
endif()
endif()
diff --git a/doc/lammps.1 b/doc/lammps.1
index 80d07cbf0a..d9af339e0a 100644
--- a/doc/lammps.1
+++ b/doc/lammps.1
@@ -1,4 +1,4 @@
-.TH LAMMPS "5 May 2020" "2020-05-5"
+.TH LAMMPS "2 June 2020" "2020-06-02"
.SH NAME
.B LAMMPS
\- Molecular Dynamics Simulator.
diff --git a/doc/src/Build_development.rst b/doc/src/Build_development.rst
index d477d2d949..94b78a4a9c 100644
--- a/doc/src/Build_development.rst
+++ b/doc/src/Build_development.rst
@@ -244,22 +244,35 @@ and working.
of mis-compiled code (or undesired large of precision due to
reordering of operations).
-------------
+Collect and visualize code coverage metrics
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-You can also collect code coverage metrics while running the tests by
-enabling coverage support during building.
+You can also collect code coverage metrics while running LAMMPS or the
+tests by enabling code coverage support during the CMake configuration:
.. code-block:: bash
- -D ENABLE_COVERAGE=value # enable coverage measurements, value = no (default) or yes
+ -D ENABLE_COVERAGE=on # enable coverage measurements (off by default)
-This will also add the following targets to generate coverage reports
-after running the LAMMPS executable or the unit tests:
+This will instrument all object files to write information about which
+lines of code were accessed during execution in files next to the
+corresponding object files. These can be post-processed to visually
+show the degree of coverage and which code paths are accessed and which
+are not taken. When working on unit tests (see above), this can be
+extremely helpful to determine which parts of the code are not executed
+and thus what kind of tests are still missing. The coverage data is
+cumulative, i.e. new data is added with each new run.
+
+Enabling code coverage will also add the following build targets to
+generate coverage reports after running the LAMMPS executable or the
+unit tests:
.. code-block:: bash
- make gen_coverage_html # generate coverage report in HTML format
- make gen_coverage_xml # generate coverage report in XML format
+ make gen_coverage_html # generate coverage report in HTML format
+ make gen_coverage_xml # generate coverage report in XML format
+ make clean_coverage_html # delete folder with HTML format coverage report
+ make reset_coverage # delete all collected coverage data and HTML output
These reports require `GCOVR `_ to be installed. The easiest way
to do this to install it via pip:
@@ -267,3 +280,29 @@ to do this to install it via pip:
.. code-block:: bash
pip install git+https://github.com/gcovr/gcovr.git
+
+After post-processing with ``gen_coverage_html`` the results are in
+a folder ``coverage_html`` and can be viewed with a web browser.
+The images below illustrate how the data is presented.
+
+.. list-table::
+
+ * - .. figure:: JPG/coverage-overview-top.png
+ :target: JPG/coverage-overview-top.png
+
+ Top of the overview page
+
+ - .. figure:: JPG/coverage-overview-manybody.png
+ :target: JPG/coverage-overview-manybody.png
+
+ Styles with good coverage
+
+ - .. figure:: JPG/coverage-file-top.png
+ :target: JPG/coverage-file-top.png
+
+ Top of individual source page
+
+ - .. figure:: JPG/coverage-file-branches.png
+ :target: JPG/coverage-file-branches.png
+
+ Source page with branches
diff --git a/doc/src/JPG/coverage-file-branches.png b/doc/src/JPG/coverage-file-branches.png
new file mode 100644
index 0000000000..8bf3b6eefd
Binary files /dev/null and b/doc/src/JPG/coverage-file-branches.png differ
diff --git a/doc/src/JPG/coverage-file-top.png b/doc/src/JPG/coverage-file-top.png
new file mode 100644
index 0000000000..cc834be27f
Binary files /dev/null and b/doc/src/JPG/coverage-file-top.png differ
diff --git a/doc/src/JPG/coverage-overview-manybody.png b/doc/src/JPG/coverage-overview-manybody.png
new file mode 100644
index 0000000000..85fad27f74
Binary files /dev/null and b/doc/src/JPG/coverage-overview-manybody.png differ
diff --git a/doc/src/JPG/coverage-overview-top.png b/doc/src/JPG/coverage-overview-top.png
new file mode 100644
index 0000000000..e57486a293
Binary files /dev/null and b/doc/src/JPG/coverage-overview-top.png differ
diff --git a/doc/src/Packages_details.rst b/doc/src/Packages_details.rst
index 8251c5301e..2353df777f 100644
--- a/doc/src/Packages_details.rst
+++ b/doc/src/Packages_details.rst
@@ -2064,7 +2064,7 @@ molecules, and chiral-sensitive reactions.
* examples/USER/reaction
* `2017 LAMMPS Workshop `_
* `2019 LAMMPS Workshop `_
-* disarmmd.org
+* reacter.org
----------
diff --git a/doc/src/fix_bond_react.rst b/doc/src/fix_bond_react.rst
index fc260de324..d5f917bfdb 100755
--- a/doc/src/fix_bond_react.rst
+++ b/doc/src/fix_bond_react.rst
@@ -300,7 +300,8 @@ either 'none' or 'charges.' Further details are provided in the
discussion of the 'update_edges' keyword. The fifth optional section
begins with the keyword 'Constraints' and lists additional criteria
that must be satisfied in order for the reaction to occur. Currently,
-there are four types of constraints available, as discussed below.
+there are four types of constraints available, as discussed below:
+'distance', 'angle', 'dihedral', and 'arrhenius'.
A sample map file is given below:
@@ -353,8 +354,9 @@ has syntax as follows:
distance *ID1* *ID2* *rmin* *rmax*
where 'distance' is the required keyword, *ID1* and *ID2* are
-pre-reaction atom IDs, and these two atoms must be separated by a
-distance between *rmin* and *rmax* for the reaction to occur.
+pre-reaction atom IDs (or molecule-fragment IDs, see below), and these
+two atoms must be separated by a distance between *rmin* and *rmax*
+for the reaction to occur.
The constraint of type 'angle' has the following syntax:
@@ -363,11 +365,11 @@ The constraint of type 'angle' has the following syntax:
angle *ID1* *ID2* *ID3* *amin* *amax*
where 'angle' is the required keyword, *ID1*\ , *ID2* and *ID3* are
-pre-reaction atom IDs, and these three atoms must form an angle
-between *amin* and *amax* for the reaction to occur (where *ID2* is
-the central atom). Angles must be specified in degrees. This
-constraint can be used to enforce a certain orientation between
-reacting molecules.
+pre-reaction atom IDs (or molecule-fragment IDs, see below), and these
+three atoms must form an angle between *amin* and *amax* for the
+reaction to occur (where *ID2* is the central atom). Angles must be
+specified in degrees. This constraint can be used to enforce a certain
+orientation between reacting molecules.
The constraint of type 'dihedral' has the following syntax:
@@ -376,15 +378,23 @@ The constraint of type 'dihedral' has the following syntax:
dihedral *ID1* *ID2* *ID3* *ID4* *amin* *amax* *amin2* *amax2*
where 'dihedral' is the required keyword, and *ID1*\ , *ID2*\ , *ID3*
-and *ID4* are pre-reaction atom IDs. Dihedral angles are calculated in
-the interval (-180,180]. Refer to the :doc:`dihedral style `
-documentation for further details on convention. If *amin* is less
-than *amax*, these four atoms must form a dihedral angle greater than
-*amin* **and** less than *amax* for the reaction to occur. If *amin*
-is greater than *amax*, these four atoms must form a dihedral angle
-greater than *amin* **or** less than *amax* for the reaction to occur.
-Angles must be specified in degrees. Optionally, a second range of
-permissible angles *amin2*-*amax2* can be specified.
+and *ID4* are pre-reaction atom IDs (or molecule-fragment IDs, see
+below). Dihedral angles are calculated in the interval (-180,180].
+Refer to the :doc:`dihedral style ` documentation for
+further details on convention. If *amin* is less than *amax*, these
+four atoms must form a dihedral angle greater than *amin* **and** less
+than *amax* for the reaction to occur. If *amin* is greater than
+*amax*, these four atoms must form a dihedral angle greater than
+*amin* **or** less than *amax* for the reaction to occur. Angles must
+be specified in degrees. Optionally, a second range of permissible
+angles *amin2*-*amax2* can be specified.
+
+For the 'distance', 'angle', and 'dihedral' constraints (explained
+above), atom IDs can be replaced by pre-reaction molecule-fragment
+IDs. The molecule-fragment ID must begin with a letter. The location
+of the ID is the geometric center of all atom positions in the
+fragment. The molecule fragment must have been defined in the
+:doc:`molecule ` command for the pre-reaction template.
The constraint of type 'arrhenius' imposes an additional reaction
probability according to the temperature-dependent Arrhenius equation:
diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt
index 4e1ca2c958..ffe18f5cb1 100644
--- a/doc/utils/sphinx-config/false_positives.txt
+++ b/doc/utils/sphinx-config/false_positives.txt
@@ -631,7 +631,6 @@ dipolar
dir
Direc
dirname
-disarmmd
discoverable
discretization
discretized
@@ -2462,6 +2461,7 @@ rdc
rdf
RDideal
rdx
+reacter
README
realtime
reamin
diff --git a/lib/kokkos/Makefile.kokkos b/lib/kokkos/Makefile.kokkos
index dc53de0a12..a9478c6e09 100644
--- a/lib/kokkos/Makefile.kokkos
+++ b/lib/kokkos/Makefile.kokkos
@@ -4,7 +4,11 @@
ifndef KOKKOS_PATH
KOKKOS_PATH=../../lib/kokkos
endif
+
CXXFLAGS=$(CCFLAGS)
+ifeq ($(mode),shared)
+CXXFLAGS += $(SHFLAGS)
+endif
KOKKOS_VERSION_MAJOR = 3
KOKKOS_VERSION_MINOR = 1
diff --git a/potentials/ffield.comb b/potentials/ffield.comb
index 0ff3c16835..9edc478afd 100644
--- a/potentials/ffield.comb
+++ b/potentials/ffield.comb
@@ -19,13 +19,13 @@
# 27 entries for a system containing three elements A, B and C
# These entries are in LAMMPS "metal" units
#
-Hf Hf Hf 1 0 1 0 1.011011 0.046511 0.959614 0.959614 55.9421 55.9421 3.85 0.20 2.069563 2.069563 707.53 707.53 0 0 0.008 0 0 0 1 1 1 1 -4 4 0.26152 -0.25918 -4 4 0.26152 -0.25918 0 3.139520d0 0 0.009410d0 0 0.679131 -3.92875 -3.92875 4.83958 4.83958 12 0
+Hf Hf Hf 1 0 1 0 1.011011 0.046511 0.959614 0.959614 55.9421 55.9421 3.85 0.20 2.069563 2.069563 707.53 707.53 0 0 0.008 0 0 0 1 1 1 1 -4 4 0.26152 -0.25918 -4 4 0.26152 -0.25918 0 3.139520 0 0.009410 0 0.679131 -3.92875 -3.92875 4.83958 4.83958 12 0
Ti Ti Ti 1 0 1 0 1.255016 0.089078 1.226342 1.226342 99.3916 99.3916 3.43 0.05 2.082408 2.082408 546.386 546.386 0 0 0.0084 0 0 0 1 1 1 1 -4 4 2.508854 -2.511416 -4 4 2.508854 -2.511416 0 2.46820415900968 0 0.151351003255176 0 0.873685 0.392632 0.392632 1.78349 1.78349 12 0
O O O 1 6.6 1 -0.229 1 2 2.68 2.68 260.893 260.893 2.8 0.2 5.36 5.36 3326.69 3326.69 0 0 0 0 0 0 1 1 1 1 -1.8349 5.5046 0.00148 -0.00112 -1.8349 5.5046 0.00148 -0.00112 5.63441383 7.689598017 4.51426991 1.330079082 0 2.243072 -3.922011 -3.922011 0.971086 0.971086 12 0
Cu Cu Cu 1 0 1 0 1 0.140835 1.681711 1.681711 146.987 146.987 3.0 0.05 2.794608 2.794608 952.693 952.693 0.077 0.0095 0 0 0 0 1 1 1 1 -6 2 0.1677645 -0.161007 -6 2 0.1677645 -0.161007 0 5.946437 0 0 0 0.454784 0.72571 0.72571 0.274649 0.274649 12 0
Si Si Si 3 100390 16.218 -0.59826 0.78734 1.0999E-06 1.7322 1.7322 471.18 471.18 2.90 0.10 2.4799 2.4799 1830.8 1830.8 0 0 0 0 0 0 1 1 1 1 -4 4 1.651725 -1.658949 -4 4 1.651725 -1.658949 0 3.625144859 0 0.087067714 0 0.772871 -0.499378 -0.499378 2.999911 2.999911 12 0
-Zr Zr Zr 1 0 1 0 1 0 0.929 0.929 39.9454 39.9454 3.8 0.31 1.857 1.857 382.6 382.6 0 0 0 0 0 0 1 1 1 1 -4 4 1.64 -1.5 -4 4 1.64 -1.5 0 3.139520d0 0 0.009410d0 0 0.679131 -3.92875 -3.92875 4.83958 4.83958 12 0
-U U U 1 0 1 0 4.346966 0.77617 0.832 0.832 162.6 162.6 3.9 0.15 1.835 1.835 795.6 795.6 0 0 0 0 0 0 1 1 1 1 -4 4 2 -2 -4 4 2 -2 0 3.139520d0 0 0.009410d0 0 0.679131 -3.92875 -3.92875 4.83958 4.83958 12 0
+Zr Zr Zr 1 0 1 0 1 0 0.929 0.929 39.9454 39.9454 3.8 0.31 1.857 1.857 382.6 382.6 0 0 0 0 0 0 1 1 1 1 -4 4 1.64 -1.5 -4 4 1.64 -1.5 0 3.139520 0 0.009410 0 0.679131 -3.92875 -3.92875 4.83958 4.83958 12 0
+U U U 1 0 1 0 4.346966 0.77617 0.832 0.832 162.6 162.6 3.9 0.15 1.835 1.835 795.6 795.6 0 0 0 0 0 0 1 1 1 1 -4 4 2 -2 -4 4 2 -2 0 3.139520 0 0.009410 0 0.679131 -3.92875 -3.92875 4.83958 4.83958 12 0
#
Si O O 3 100390 16.218 -0.59826 0.78734 1.0999E-06 1.7322 2.68 471.18 260.893 2.80 0.25 2.4799 5.36 1830.8 3326.69 0 0 0 109.47 0.3122 0 1 1 1 1 -4 4 1.651725 -1.658949 -1.8349 5.5046 0.00148 -0.00112 0 3.625144859 0 0.087067714 0 0.772871 -0.499378 -3.922011 2.999911 0.971086 12 0
O Si Si 1 6.6 1 -0.229 1 2 2.68 1.7322 260.893 471.18 2.80 0.25 5.36 2.4799 3326.69 1830.8 0 0 0 143.73 2.6 0 1 1 1 1 -1.8349 5.5046 0.00148 -0.00112 -4 4 1.651725 -1.658949 5.63441383 7.689598017 4.51426991 1.330079082 0 2.243072 -3.922011 -0.499378 0.971086 2.999911 12 0
diff --git a/src/GPU/pair_eam_alloy_gpu.cpp b/src/GPU/pair_eam_alloy_gpu.cpp
index 266a79bc3c..fcbc3c5d3c 100644
--- a/src/GPU/pair_eam_alloy_gpu.cpp
+++ b/src/GPU/pair_eam_alloy_gpu.cpp
@@ -31,11 +31,11 @@
#include "domain.h"
#include "utils.h"
#include "suffix.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
using namespace LAMMPS_NS;
-#define MAXLINE 1024
-
// External functions from cuda library for atom decomposition
int eam_alloy_gpu_init(const int ntypes, double host_cutforcesq,
@@ -369,94 +369,112 @@ void PairEAMAlloyGPU::read_file(char *filename)
{
Setfl *file = setfl;
- // open potential file
+ // read potential file
+ if(comm->me == 0) {
+ PotentialFileReader reader(lmp, filename, "EAM");
- int me = comm->me;
- FILE *fptr;
- char line[MAXLINE];
+ try {
+ char * line = nullptr;
- if (me == 0) {
- fptr = fopen(filename,"r");
- if (fptr == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open EAM potential file %s",filename);
- error->one(FLERR,str);
+ reader.skip_line();
+ reader.skip_line();
+ reader.skip_line();
+
+ // extract element names from nelements line
+ line = reader.next_line(1);
+ ValueTokenizer values(line);
+ file->nelements = values.next_int();
+
+ if (values.count() != file->nelements + 1)
+ error->one(FLERR,"Incorrect element names in EAM potential file");
+
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) {
+ const std::string word = values.next_string();
+ const int n = word.length() + 1;
+ file->elements[i] = new char[n];
+ strcpy(file->elements[i], word.c_str());
+ }
+
+ //
+
+ line = reader.next_line(5);
+ values = ValueTokenizer(line);
+ file->nrho = values.next_int();
+ file->drho = values.next_double();
+ file->nr = values.next_int();
+ file->dr = values.next_double();
+ file->cut = values.next_double();
+
+ if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
+ error->one(FLERR,"Invalid EAM potential file");
+
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
+
+ for (int i = 0; i < file->nelements; i++) {
+ line = reader.next_line(2);
+ values = ValueTokenizer(line);
+ values.next_int(); // ignore
+ file->mass[i] = values.next_double();
+
+ reader.next_dvector(file->nrho, &file->frho[i][1]);
+ reader.next_dvector(file->nr, &file->rhor[i][1]);
+ }
+
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ reader.next_dvector(file->nr, &file->z2r[i][j][1]);
+ }
+ }
+ } catch (TokenizerException & e) {
+ error->one(FLERR, e.what());
}
}
- // read and broadcast header
- // extract element names from nelements line
+ // broadcast potential information
+ MPI_Bcast(&file->nelements, 1, MPI_INT, 0, world);
- int n;
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- n = strlen(line) + 1;
+ MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->nr, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world);
+
+ // allocate memory on other procs
+ if (comm->me != 0) {
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) file->elements[i] = nullptr;
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
}
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
- sscanf(line,"%d",&file->nelements);
- int nwords = atom->count_words(line);
- if (nwords != file->nelements + 1)
- error->all(FLERR,"Incorrect element names in EAM potential file");
-
- char **words = new char*[file->nelements+1];
- nwords = 0;
- strtok(line," \t\n\r\f");
- while ( (words[nwords++] = strtok(NULL," \t\n\r\f")) ) continue;
-
- file->elements = new char*[file->nelements];
+ // broadcast file->elements string array
for (int i = 0; i < file->nelements; i++) {
- n = strlen(words[i]) + 1;
- file->elements[i] = new char[n];
- strcpy(file->elements[i],words[i]);
- }
- delete [] words;
-
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- sscanf(line,"%d %lg %d %lg %lg",
- &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);
+ int n;
+ if (comm->me == 0) n = strlen(file->elements[i]) + 1;
+ MPI_Bcast(&n, 1, MPI_INT, 0, world);
+ if (comm->me != 0) file->elements[i] = new char[n];
+ MPI_Bcast(file->elements[i], n, MPI_CHAR, 0, world);
}
- MPI_Bcast(&file->nrho,1,MPI_INT,0,world);
- MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->nr,1,MPI_INT,0,world);
- MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world);
+ // broadcast file->mass, frho, rhor
+ for (int i = 0; i < file->nelements; i++) {
+ MPI_Bcast(&file->mass[i], 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->frho[i][1], file->nrho, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->rhor[i][1], file->nr, MPI_DOUBLE, 0, world);
+ }
- file->mass = new double[file->nelements];
- memory->create(file->frho,file->nelements,file->nrho+1,"pair:frho");
- memory->create(file->rhor,file->nelements,file->nr+1,"pair:rhor");
- memory->create(file->z2r,file->nelements,file->nelements,file->nr+1,
- "pair:z2r");
-
- int i,j,tmp;
- for (i = 0; i < file->nelements; i++) {
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- sscanf(line,"%d %lg",&tmp,&file->mass[i]);
+ // broadcast file->z2r
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world);
}
- MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world);
-
- if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]);
- MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world);
- if (me == 0) grab(fptr,file->nr,&file->rhor[i][1]);
- MPI_Bcast(&file->rhor[i][1],file->nr,MPI_DOUBLE,0,world);
}
-
- for (i = 0; i < file->nelements; i++)
- for (j = 0; j <= i; j++) {
- if (me == 0) grab(fptr,file->nr,&file->z2r[i][j][1]);
- MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
- }
-
- // close the potential file
-
- if (me == 0) fclose(fptr);
}
/* ----------------------------------------------------------------------
diff --git a/src/GPU/pair_eam_fs_gpu.cpp b/src/GPU/pair_eam_fs_gpu.cpp
index 2bd5a53324..e674594498 100644
--- a/src/GPU/pair_eam_fs_gpu.cpp
+++ b/src/GPU/pair_eam_fs_gpu.cpp
@@ -30,11 +30,12 @@
#include "gpu_extra.h"
#include "domain.h"
#include "suffix.h"
+#include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
using namespace LAMMPS_NS;
-#define MAXLINE 1024
-
// External functions from cuda library for atom decomposition
int eam_fs_gpu_init(const int ntypes, double host_cutforcesq,
@@ -368,99 +369,118 @@ void PairEAMFSGPU::read_file(char *filename)
{
Fs *file = fs;
- // open potential file
+ // read potential file
+ if(comm->me == 0) {
+ PotentialFileReader reader(lmp, filename, "EAM");
- int me = comm->me;
- FILE *fptr;
- char line[MAXLINE];
+ try {
+ char * line = nullptr;
- if (me == 0) {
- fptr = force->open_potential(filename);
- if (fptr == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open EAM potential file %s",filename);
- error->one(FLERR,str);
+ reader.skip_line();
+ reader.skip_line();
+ reader.skip_line();
+
+ // extract element names from nelements line
+ line = reader.next_line(1);
+ ValueTokenizer values(line);
+ file->nelements = values.next_int();
+
+ if (values.count() != file->nelements + 1)
+ error->one(FLERR,"Incorrect element names in EAM potential file");
+
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) {
+ const std::string word = values.next_string();
+ const int n = word.length() + 1;
+ file->elements[i] = new char[n];
+ strcpy(file->elements[i], word.c_str());
+ }
+
+ //
+
+ line = reader.next_line(5);
+ values = ValueTokenizer(line);
+ file->nrho = values.next_int();
+ file->drho = values.next_double();
+ file->nr = values.next_int();
+ file->dr = values.next_double();
+ file->cut = values.next_double();
+
+ if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
+ error->one(FLERR,"Invalid EAM potential file");
+
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
+
+ for (int i = 0; i < file->nelements; i++) {
+ line = reader.next_line(2);
+ values = ValueTokenizer(line);
+ values.next_int(); // ignore
+ file->mass[i] = values.next_double();
+
+ reader.next_dvector(file->nrho, &file->frho[i][1]);
+
+ for (int j = 0; j < file->nelements; j++) {
+ reader.next_dvector(file->nr, &file->rhor[i][j][1]);
+ }
+ }
+
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ reader.next_dvector(file->nr, &file->z2r[i][j][1]);
+ }
+ }
+ } catch (TokenizerException & e) {
+ error->one(FLERR, e.what());
}
}
- // read and broadcast header
- // extract element names from nelements line
+ // broadcast potential information
+ MPI_Bcast(&file->nelements, 1, MPI_INT, 0, world);
- int n;
- if (me == 0) {
- fgets(line,MAXLINE,fptr);
- fgets(line,MAXLINE,fptr);
- fgets(line,MAXLINE,fptr);
- fgets(line,MAXLINE,fptr);
- n = strlen(line) + 1;
+ MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->nr, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world);
+
+ // allocate memory on other procs
+ if (comm->me != 0) {
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) file->elements[i] = nullptr;
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
}
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
- sscanf(line,"%d",&file->nelements);
- int nwords = atom->count_words(line);
- if (nwords != file->nelements + 1)
- error->all(FLERR,"Incorrect element names in EAM potential file");
-
- char **words = new char*[file->nelements+1];
- nwords = 0;
- strtok(line," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
-
- file->elements = new char*[file->nelements];
+ // broadcast file->elements string array
for (int i = 0; i < file->nelements; i++) {
- n = strlen(words[i]) + 1;
- file->elements[i] = new char[n];
- strcpy(file->elements[i],words[i]);
- }
- delete [] words;
-
- if (me == 0) {
- fgets(line,MAXLINE,fptr);
- sscanf(line,"%d %lg %d %lg %lg",
- &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);
+ int n;
+ if (comm->me == 0) n = strlen(file->elements[i]) + 1;
+ MPI_Bcast(&n, 1, MPI_INT, 0, world);
+ if (comm->me != 0) file->elements[i] = new char[n];
+ MPI_Bcast(file->elements[i], n, MPI_CHAR, 0, world);
}
- MPI_Bcast(&file->nrho,1,MPI_INT,0,world);
- MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->nr,1,MPI_INT,0,world);
- MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world);
+ // broadcast file->mass, frho, rhor
+ for (int i = 0; i < file->nelements; i++) {
+ MPI_Bcast(&file->mass[i], 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->frho[i][1], file->nrho, MPI_DOUBLE, 0, world);
- file->mass = new double[file->nelements];
- memory->create(file->frho,file->nelements,file->nrho+1,
- "pair:frho");
- memory->create(file->rhor,file->nelements,file->nelements,
- file->nr+1,"pair:rhor");
- memory->create(file->z2r,file->nelements,file->nelements,
- file->nr+1,"pair:z2r");
-
- int i,j,tmp;
- for (i = 0; i < file->nelements; i++) {
- if (me == 0) {
- fgets(line,MAXLINE,fptr);
- sscanf(line,"%d %lg",&tmp,&file->mass[i]);
- }
- MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world);
-
- if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]);
- MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world);
-
- for (j = 0; j < file->nelements; j++) {
- if (me == 0) grab(fptr,file->nr,&file->rhor[i][j][1]);
- MPI_Bcast(&file->rhor[i][j][1],file->nr,MPI_DOUBLE,0,world);
+ for (int j = 0; j < file->nelements; j++) {
+ MPI_Bcast(&file->rhor[i][j][1], file->nr, MPI_DOUBLE, 0, world);
}
}
- for (i = 0; i < file->nelements; i++)
- for (j = 0; j <= i; j++) {
- if (me == 0) grab(fptr,file->nr,&file->z2r[i][j][1]);
- MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
+ // broadcast file->z2r
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world);
}
-
- // close the potential file
-
- if (me == 0) fclose(fptr);
+ }
}
/* ----------------------------------------------------------------------
diff --git a/src/KIM/kim_init.cpp b/src/KIM/kim_init.cpp
index 1d3830c275..5328769d22 100644
--- a/src/KIM/kim_init.cpp
+++ b/src/KIM/kim_init.cpp
@@ -275,7 +275,7 @@ void KimInit::determine_model_type_and_units(char * model_name,
std::string mesg("Incompatible units for KIM Simulator Model, "
"required units = ");
mesg += *model_units;
- error->all(FLERR,mesg.c_str());
+ error->all(FLERR,mesg);
}
}
}
@@ -328,8 +328,7 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units, KIM
mesg += "\n";
mesg += "#\n";
- if (screen) fputs(mesg.c_str(),screen);
- if (logfile) fputs(mesg.c_str(),logfile);
+ utils::logmesg(lmp,mesg);
}
fix_store->setptr("simulator_model", (void *) simulatorModel);
@@ -346,7 +345,7 @@ void KimInit::do_init(char *model_name, char *user_units, char *model_units, KIM
std::string cmd("units ");
cmd += model_units;
- input->one(cmd.c_str());
+ input->one(cmd);
if (model_type == SM) {
int sim_fields, sim_lines;
@@ -519,7 +518,7 @@ void KimInit::do_variables(char *user_units, char *model_units)
"unit = " + units[i] + "; "
"from = " + from + "; "
"to = " + to + ".";
- error->all(FLERR,err.c_str());
+ error->all(FLERR,err);
}
variable->internal_set(v_unit,conversion_factor);
if (comm->me == 0) {
diff --git a/src/KIM/kim_interactions.cpp b/src/KIM/kim_interactions.cpp
index d0aa0e3c4e..2167cd6069 100644
--- a/src/KIM/kim_interactions.cpp
+++ b/src/KIM/kim_interactions.cpp
@@ -192,7 +192,7 @@ void KimInteractions::do_setup(int narg, char **arg)
std::string msg("Species '");
msg += strword;
msg += "' is not supported by this KIM Simulator Model";
- error->all(FLERR,msg.c_str());
+ error->all(FLERR,msg);
}
strword = strtok(NULL," \t");
}
@@ -276,8 +276,8 @@ void KimInteractions::do_setup(int narg, char **arg)
cmd2 += " ";
}
- input->one(cmd1.c_str());
- input->one(cmd2.c_str());
+ input->one(cmd1);
+ input->one(cmd2);
}
// End output to log file
diff --git a/src/KIM/kim_param.cpp b/src/KIM/kim_param.cpp
index f5a154f2d9..ac995c8367 100644
--- a/src/KIM/kim_param.cpp
+++ b/src/KIM/kim_param.cpp
@@ -163,7 +163,7 @@ void KimParam::command(int narg, char **arg)
if (!kim_param_get && !kim_param_set) {
std::string msg("Incorrect arguments in kim_param command.\n");
msg += "'kim_param get/set' is mandatory.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
// Check if we called a kim_init command
@@ -225,7 +225,7 @@ void KimParam::command(int narg, char **arg)
msg += "To set the new parameter values, pair style must be assigned.\n";
msg += "Must use 'kim_interactions' or";
msg += "'pair_style kim ' before 'kim_param set'";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
} else {
KIM_LengthUnit lengthUnit;
KIM_EnergyUnit energyUnit;
@@ -296,7 +296,7 @@ void KimParam::command(int narg, char **arg)
msg += "This Model does not have the requested '";
msg += paramname;
msg += "' parameter.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
// Get the index_range for the requested parameter
@@ -313,7 +313,7 @@ void KimParam::command(int narg, char **arg)
msg += "Expected integer parameter(s) instead of '";
msg += argtostr;
msg += "' in index_range.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
std::string::size_type npos = argtostr.find(':');
@@ -330,7 +330,7 @@ void KimParam::command(int narg, char **arg)
msg += "' parameter with extent of '";
msg += SNUM(extent);
msg += "' .";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
} else {
std::stringstream str(argtostr);
@@ -342,7 +342,7 @@ void KimParam::command(int narg, char **arg)
msg += "' with the extent of '";
msg += SNUM(extent);
msg += "' .";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
nubound = nlbound;
}
@@ -350,7 +350,7 @@ void KimParam::command(int narg, char **arg)
std::string msg("Wrong number of arguments in ");
msg += "kim_param get command.\n";
msg += "Index range after parameter name is mandatory.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
int const nvars = nubound - nlbound + 1;
@@ -363,7 +363,7 @@ void KimParam::command(int narg, char **arg)
std::string msg("Wrong number of arguments in ");
msg += "kim_param get command.\n";
msg += "The LAMMPS variable name is mandatory.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
// indicator flag for list request
@@ -400,7 +400,7 @@ void KimParam::command(int narg, char **arg)
msg += "' variable names or '";
msg += varname;
msg += " split' is mandatory.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
} else {
std::string msg("Wrong number of arguments in ");
@@ -410,7 +410,7 @@ void KimParam::command(int narg, char **arg)
msg += "' variable names or '";
msg += varname;
msg += " split/list' is mandatory.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
} else {
varsname = new char *[1];
@@ -524,7 +524,7 @@ void KimParam::command(int narg, char **arg)
set_cmd += " ";
set_cmd += arg[i];
}
- input->one(set_cmd.c_str());
+ input->one(set_cmd);
}
} else
error->all(FLERR, "This model has No mutable parameters.");
diff --git a/src/KIM/kim_property.cpp b/src/KIM/kim_property.cpp
index 4e63a25ffb..1020f380b7 100644
--- a/src/KIM/kim_property.cpp
+++ b/src/KIM/kim_property.cpp
@@ -100,13 +100,13 @@ void kimProperty::command(int narg, char **arg)
std::string msg("Error incorrect arguments in kim_property command.\n");
msg += "`kim_property create/destroy/modify/remove/dump` ";
msg += "is mandatory.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
if (comm->me == 0) {
std::string msg;
msg = "#=== kim-property ===========================================\n";
- input->write_echo(msg.c_str());
+ input->write_echo(msg);
}
// Get the kim_str ptr to the data associated with a kim_property_str
diff --git a/src/KIM/pair_kim.cpp b/src/KIM/pair_kim.cpp
index 8ff0ddf455..76c18bb38a 100644
--- a/src/KIM/pair_kim.cpp
+++ b/src/KIM/pair_kim.cpp
@@ -429,7 +429,7 @@ void PairKIM::coeff(int narg, char **arg)
} else {
std::string msg("create_kim_particle_codes: symbol not found: ");
msg += lmps_unique_elements[i];
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
}
// Set the new values for PM parameters
@@ -441,7 +441,7 @@ void PairKIM::coeff(int narg, char **arg)
if (!numberOfParameters) {
std::string msg("Incorrect args for pair coefficients \n");
msg += "This model has No mutable parameters.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
int kimerror;
@@ -477,7 +477,7 @@ void PairKIM::coeff(int narg, char **arg)
msg += "This Model does not have the requested '";
msg += paramname;
msg += "' parameter.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
// Get the index_range for the requested parameter
@@ -493,7 +493,7 @@ void PairKIM::coeff(int narg, char **arg)
msg += "Expected integer parameter(s) instead of '";
msg += argtostr;
msg += "' in index_range.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
std::string::size_type npos = argtostr.find(':');
@@ -510,7 +510,7 @@ void PairKIM::coeff(int narg, char **arg)
msg += "' parameter with extent of '";
msg += SNUM(extent);
msg += "' .";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
} else {
std::stringstream str(argtostr);
@@ -522,7 +522,7 @@ void PairKIM::coeff(int narg, char **arg)
msg += "' parameter with extent of '";
msg += SNUM(extent);
msg += "' .";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
nubound = nlbound;
}
@@ -530,7 +530,7 @@ void PairKIM::coeff(int narg, char **arg)
std::string msg =
"Wrong number of arguments for pair coefficients.\n";
msg += "Index range after parameter name is mandatory.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
// Parameter values
@@ -561,7 +561,7 @@ void PairKIM::coeff(int narg, char **arg)
msg += "' values are requested for '";
msg += paramname;
msg += "' parameter.";
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
}
@@ -1122,7 +1122,7 @@ void PairKIM::set_kim_model_has_flags()
KIM_SUPPORT_STATUS_required)) {
std::string msg("KIM Model requires unsupported compute argument: ");
msg += KIM_ComputeArgumentName_ToString(computeArgumentName);
- error->all(FLERR, msg.c_str());
+ error->all(FLERR, msg);
}
}
diff --git a/src/KOKKOS/atom_kokkos.cpp b/src/KOKKOS/atom_kokkos.cpp
index 49fdbb95bb..e826a7c392 100644
--- a/src/KOKKOS/atom_kokkos.cpp
+++ b/src/KOKKOS/atom_kokkos.cpp
@@ -340,7 +340,8 @@ void AtomKokkos::sync_modify(ExecutionSpace execution_space,
modified(execution_space,datamask_modify);
}
-AtomVec *AtomKokkos::new_avec(const char *style, int trysuffix, int &sflag)
+AtomVec *AtomKokkos::new_avec(const std::string &style,
+ int trysuffix, int &sflag)
{
AtomVec* avec = Atom::new_avec(style,trysuffix,sflag);
if (!avec->kokkosable)
diff --git a/src/KOKKOS/atom_kokkos.h b/src/KOKKOS/atom_kokkos.h
index 0ae032032a..82ea55beb1 100644
--- a/src/KOKKOS/atom_kokkos.h
+++ b/src/KOKKOS/atom_kokkos.h
@@ -74,7 +74,7 @@ class AtomKokkos : public Atom {
virtual void deallocate_topology();
void sync_modify(ExecutionSpace, unsigned int, unsigned int);
private:
- class AtomVec *new_avec(const char *, int, int &);
+ class AtomVec *new_avec(const std::string &, int, int &);
};
template
diff --git a/src/KOKKOS/pair_eam_alloy_kokkos.cpp b/src/KOKKOS/pair_eam_alloy_kokkos.cpp
index 1963722d86..7f598f966d 100644
--- a/src/KOKKOS/pair_eam_alloy_kokkos.cpp
+++ b/src/KOKKOS/pair_eam_alloy_kokkos.cpp
@@ -29,11 +29,12 @@
#include "memory_kokkos.h"
#include "error.h"
#include "atom_masks.h"
+#include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
using namespace LAMMPS_NS;
-#define MAXLINE 1024
-
// Cannot use virtual inheritance on the GPU, so must duplicate code
/* ---------------------------------------------------------------------- */
@@ -982,94 +983,112 @@ void PairEAMAlloyKokkos::read_file(char *filename)
{
Setfl *file = setfl;
- // open potential file
+ // read potential file
+ if(comm->me == 0) {
+ PotentialFileReader reader(lmp, filename, "EAM");
- int me = comm->me;
- FILE *fptr;
- char line[MAXLINE];
+ try {
+ char * line = nullptr;
- if (me == 0) {
- fptr = force->open_potential(filename);
- if (fptr == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open EAM potential file %s",filename);
- error->one(FLERR,str);
+ reader.skip_line();
+ reader.skip_line();
+ reader.skip_line();
+
+ // extract element names from nelements line
+ line = reader.next_line(1);
+ ValueTokenizer values(line);
+ file->nelements = values.next_int();
+
+ if (values.count() != file->nelements + 1)
+ error->one(FLERR,"Incorrect element names in EAM potential file");
+
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) {
+ const std::string word = values.next_string();
+ const int n = word.length() + 1;
+ file->elements[i] = new char[n];
+ strcpy(file->elements[i], word.c_str());
+ }
+
+ //
+
+ line = reader.next_line(5);
+ values = ValueTokenizer(line);
+ file->nrho = values.next_int();
+ file->drho = values.next_double();
+ file->nr = values.next_int();
+ file->dr = values.next_double();
+ file->cut = values.next_double();
+
+ if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
+ error->one(FLERR,"Invalid EAM potential file");
+
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
+
+ for (int i = 0; i < file->nelements; i++) {
+ line = reader.next_line(2);
+ values = ValueTokenizer(line);
+ values.next_int(); // ignore
+ file->mass[i] = values.next_double();
+
+ reader.next_dvector(file->nrho, &file->frho[i][1]);
+ reader.next_dvector(file->nr, &file->rhor[i][1]);
+ }
+
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ reader.next_dvector(file->nr, &file->z2r[i][j][1]);
+ }
+ }
+ } catch (TokenizerException & e) {
+ error->one(FLERR, e.what());
}
}
- // read and broadcast header
- // extract element names from nelements line
+ // broadcast potential information
+ MPI_Bcast(&file->nelements, 1, MPI_INT, 0, world);
- int n;
- if (me == 0) {
- fgets(line,MAXLINE,fptr);
- fgets(line,MAXLINE,fptr);
- fgets(line,MAXLINE,fptr);
- fgets(line,MAXLINE,fptr);
- n = strlen(line) + 1;
+ MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->nr, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world);
+
+ // allocate memory on other procs
+ if (comm->me != 0) {
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) file->elements[i] = nullptr;
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
}
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
- sscanf(line,"%d",&file->nelements);
- int nwords = atom->count_words(line);
- if (nwords != file->nelements + 1)
- error->all(FLERR,"Incorrect element names in EAM potential file");
-
- char **words = new char*[file->nelements+1];
- nwords = 0;
- strtok(line," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
-
- file->elements = new char*[file->nelements];
+ // broadcast file->elements string array
for (int i = 0; i < file->nelements; i++) {
- n = strlen(words[i]) + 1;
- file->elements[i] = new char[n];
- strcpy(file->elements[i],words[i]);
- }
- delete [] words;
-
- if (me == 0) {
- fgets(line,MAXLINE,fptr);
- sscanf(line,"%d %lg %d %lg %lg",
- &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);
+ int n;
+ if (comm->me == 0) n = strlen(file->elements[i]) + 1;
+ MPI_Bcast(&n, 1, MPI_INT, 0, world);
+ if (comm->me != 0) file->elements[i] = new char[n];
+ MPI_Bcast(file->elements[i], n, MPI_CHAR, 0, world);
}
- MPI_Bcast(&file->nrho,1,MPI_INT,0,world);
- MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->nr,1,MPI_INT,0,world);
- MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world);
+ // broadcast file->mass, frho, rhor
+ for (int i = 0; i < file->nelements; i++) {
+ MPI_Bcast(&file->mass[i], 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->frho[i][1], file->nrho, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->rhor[i][1], file->nr, MPI_DOUBLE, 0, world);
+ }
- file->mass = new double[file->nelements];
- memory->create(file->frho,file->nelements,file->nrho+1,"pair:frho");
- memory->create(file->rhor,file->nelements,file->nr+1,"pair:rhor");
- memory->create(file->z2r,file->nelements,file->nelements,file->nr+1,
- "pair:z2r");
-
- int i,j,tmp;
- for (i = 0; i < file->nelements; i++) {
- if (me == 0) {
- fgets(line,MAXLINE,fptr);
- sscanf(line,"%d %lg",&tmp,&file->mass[i]);
+ // broadcast file->z2r
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world);
}
- MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world);
-
- if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]);
- MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world);
- if (me == 0) grab(fptr,file->nr,&file->rhor[i][1]);
- MPI_Bcast(&file->rhor[i][1],file->nr,MPI_DOUBLE,0,world);
}
-
- for (i = 0; i < file->nelements; i++)
- for (j = 0; j <= i; j++) {
- if (me == 0) grab(fptr,file->nr,&file->z2r[i][j][1]);
- MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
- }
-
- // close the potential file
-
- if (me == 0) fclose(fptr);
}
/* ----------------------------------------------------------------------
diff --git a/src/KOKKOS/pair_eam_fs_kokkos.cpp b/src/KOKKOS/pair_eam_fs_kokkos.cpp
index a33923887b..e532a4ec38 100644
--- a/src/KOKKOS/pair_eam_fs_kokkos.cpp
+++ b/src/KOKKOS/pair_eam_fs_kokkos.cpp
@@ -29,11 +29,12 @@
#include "memory_kokkos.h"
#include "error.h"
#include "atom_masks.h"
+#include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
using namespace LAMMPS_NS;
-#define MAXLINE 1024
-
// Cannot use virtual inheritance on the GPU, so must duplicate code
/* ---------------------------------------------------------------------- */
@@ -982,99 +983,118 @@ void PairEAMFSKokkos::read_file(char *filename)
{
Fs *file = fs;
- // open potential file
+ // read potential file
+ if(comm->me == 0) {
+ PotentialFileReader reader(lmp, filename, "EAM");
- int me = comm->me;
- FILE *fptr;
- char line[MAXLINE];
+ try {
+ char * line = nullptr;
- if (me == 0) {
- fptr = force->open_potential(filename);
- if (fptr == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open EAM potential file %s",filename);
- error->one(FLERR,str);
+ reader.skip_line();
+ reader.skip_line();
+ reader.skip_line();
+
+ // extract element names from nelements line
+ line = reader.next_line(1);
+ ValueTokenizer values(line);
+ file->nelements = values.next_int();
+
+ if (values.count() != file->nelements + 1)
+ error->one(FLERR,"Incorrect element names in EAM potential file");
+
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) {
+ const std::string word = values.next_string();
+ const int n = word.length() + 1;
+ file->elements[i] = new char[n];
+ strcpy(file->elements[i], word.c_str());
+ }
+
+ //
+
+ line = reader.next_line(5);
+ values = ValueTokenizer(line);
+ file->nrho = values.next_int();
+ file->drho = values.next_double();
+ file->nr = values.next_int();
+ file->dr = values.next_double();
+ file->cut = values.next_double();
+
+ if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
+ error->one(FLERR,"Invalid EAM potential file");
+
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
+
+ for (int i = 0; i < file->nelements; i++) {
+ line = reader.next_line(2);
+ values = ValueTokenizer(line);
+ values.next_int(); // ignore
+ file->mass[i] = values.next_double();
+
+ reader.next_dvector(file->nrho, &file->frho[i][1]);
+
+ for (int j = 0; j < file->nelements; j++) {
+ reader.next_dvector(file->nr, &file->rhor[i][j][1]);
+ }
+ }
+
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ reader.next_dvector(file->nr, &file->z2r[i][j][1]);
+ }
+ }
+ } catch (TokenizerException & e) {
+ error->one(FLERR, e.what());
}
}
- // read and broadcast header
- // extract element names from nelements line
+ // broadcast potential information
+ MPI_Bcast(&file->nelements, 1, MPI_INT, 0, world);
- int n;
- if (me == 0) {
- fgets(line,MAXLINE,fptr);
- fgets(line,MAXLINE,fptr);
- fgets(line,MAXLINE,fptr);
- fgets(line,MAXLINE,fptr);
- n = strlen(line) + 1;
+ MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->nr, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world);
+
+ // allocate memory on other procs
+ if (comm->me != 0) {
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) file->elements[i] = nullptr;
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
}
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
- sscanf(line,"%d",&file->nelements);
- int nwords = atom->count_words(line);
- if (nwords != file->nelements + 1)
- error->all(FLERR,"Incorrect element names in EAM potential file");
-
- char **words = new char*[file->nelements+1];
- nwords = 0;
- strtok(line," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
-
- file->elements = new char*[file->nelements];
+ // broadcast file->elements string array
for (int i = 0; i < file->nelements; i++) {
- n = strlen(words[i]) + 1;
- file->elements[i] = new char[n];
- strcpy(file->elements[i],words[i]);
- }
- delete [] words;
-
- if (me == 0) {
- fgets(line,MAXLINE,fptr);
- sscanf(line,"%d %lg %d %lg %lg",
- &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);
+ int n;
+ if (comm->me == 0) n = strlen(file->elements[i]) + 1;
+ MPI_Bcast(&n, 1, MPI_INT, 0, world);
+ if (comm->me != 0) file->elements[i] = new char[n];
+ MPI_Bcast(file->elements[i], n, MPI_CHAR, 0, world);
}
- MPI_Bcast(&file->nrho,1,MPI_INT,0,world);
- MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->nr,1,MPI_INT,0,world);
- MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world);
+ // broadcast file->mass, frho, rhor
+ for (int i = 0; i < file->nelements; i++) {
+ MPI_Bcast(&file->mass[i], 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->frho[i][1], file->nrho, MPI_DOUBLE, 0, world);
- file->mass = new double[file->nelements];
- memory->create(file->frho,file->nelements,file->nrho+1,
- "pair:frho");
- memory->create(file->rhor,file->nelements,file->nelements,
- file->nr+1,"pair:rhor");
- memory->create(file->z2r,file->nelements,file->nelements,
- file->nr+1,"pair:z2r");
-
- int i,j,tmp;
- for (i = 0; i < file->nelements; i++) {
- if (me == 0) {
- fgets(line,MAXLINE,fptr);
- sscanf(line,"%d %lg",&tmp,&file->mass[i]);
- }
- MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world);
-
- if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]);
- MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world);
-
- for (j = 0; j < file->nelements; j++) {
- if (me == 0) grab(fptr,file->nr,&file->rhor[i][j][1]);
- MPI_Bcast(&file->rhor[i][j][1],file->nr,MPI_DOUBLE,0,world);
+ for (int j = 0; j < file->nelements; j++) {
+ MPI_Bcast(&file->rhor[i][j][1], file->nr, MPI_DOUBLE, 0, world);
}
}
- for (i = 0; i < file->nelements; i++)
- for (j = 0; j <= i; j++) {
- if (me == 0) grab(fptr,file->nr,&file->z2r[i][j][1]);
- MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
+ // broadcast file->z2r
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world);
}
-
- // close the potential file
-
- if (me == 0) fclose(fptr);
+ }
}
/* ----------------------------------------------------------------------
diff --git a/src/KOKKOS/pair_exp6_rx_kokkos.cpp b/src/KOKKOS/pair_exp6_rx_kokkos.cpp
index a64df08127..8c98e32d5b 100644
--- a/src/KOKKOS/pair_exp6_rx_kokkos.cpp
+++ b/src/KOKKOS/pair_exp6_rx_kokkos.cpp
@@ -32,6 +32,7 @@
#include "neigh_request.h"
#include "atom_kokkos.h"
#include "kokkos.h"
+#include "utils.h"
#ifdef _OPENMP
#include
@@ -1752,7 +1753,7 @@ void PairExp6rxKokkos::read_file(char *file)
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
- nwords = atom->count_words(line);
+ nwords = utils::count_words(line);
if (nwords == 0) continue;
// concatenate additional lines until have params_per_line words
@@ -1771,7 +1772,7 @@ void PairExp6rxKokkos::read_file(char *file)
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
if ((ptr = strchr(line,'#'))) *ptr = '\0';
- nwords = atom->count_words(line);
+ nwords = utils::count_words(line);
}
if (nwords != params_per_line)
diff --git a/src/KOKKOS/pair_kokkos.h b/src/KOKKOS/pair_kokkos.h
index c28c858113..2934ddf621 100644
--- a/src/KOKKOS/pair_kokkos.h
+++ b/src/KOKKOS/pair_kokkos.h
@@ -51,7 +51,8 @@ struct DoCoul<1> {
//Specialisation for Neighborlist types Half, HalfThread, Full
template
-struct PairComputeFunctor {
+class PairComputeFunctor {
+ public:
typedef typename PairStyle::device_type device_type ;
typedef ArrayTypes AT;
diff --git a/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp b/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp
index 0371ac42c5..36d8126eec 100644
--- a/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp
+++ b/src/KSPACE/pair_lj_charmmfsw_coul_long.cpp
@@ -739,37 +739,18 @@ void PairLJCharmmfswCoulLong::init_style()
int irequest;
+ int respa = 0;
if (update->whichflag == 1 && strstr(update->integrate_style,"respa")) {
- int respa = 0;
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
+ }
- if (respa == 0) irequest = neighbor->request(this,instance_me);
- else if (respa == 1) {
- irequest = neighbor->request(this,instance_me);
- neighbor->requests[irequest]->id = 1;
- neighbor->requests[irequest]->half = 0;
- neighbor->requests[irequest]->respainner = 1;
- irequest = neighbor->request(this,instance_me);
- neighbor->requests[irequest]->id = 3;
- neighbor->requests[irequest]->half = 0;
- neighbor->requests[irequest]->respaouter = 1;
- } else {
- irequest = neighbor->request(this,instance_me);
- neighbor->requests[irequest]->id = 1;
- neighbor->requests[irequest]->half = 0;
- neighbor->requests[irequest]->respainner = 1;
- irequest = neighbor->request(this,instance_me);
- neighbor->requests[irequest]->id = 2;
- neighbor->requests[irequest]->half = 0;
- neighbor->requests[irequest]->respamiddle = 1;
- irequest = neighbor->request(this,instance_me);
- neighbor->requests[irequest]->id = 3;
- neighbor->requests[irequest]->half = 0;
- neighbor->requests[irequest]->respaouter = 1;
- }
-
- } else irequest = neighbor->request(this,instance_me);
+ irequest = neighbor->request(this,instance_me);
+ if (respa >= 1) {
+ neighbor->requests[irequest]->respaouter = 1;
+ neighbor->requests[irequest]->respainner = 1;
+ }
+ if (respa == 2) neighbor->requests[irequest]->respamiddle = 1;
// require cut_lj_inner < cut_lj
diff --git a/src/MANYBODY/pair_adp.cpp b/src/MANYBODY/pair_adp.cpp
index e730466fcb..2a64839e98 100644
--- a/src/MANYBODY/pair_adp.cpp
+++ b/src/MANYBODY/pair_adp.cpp
@@ -29,11 +29,11 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
using namespace LAMMPS_NS;
-#define MAXLINE 1024
-
/* ---------------------------------------------------------------------- */
PairADP::PairADP(LAMMPS *lmp) : Pair(lmp)
@@ -96,7 +96,7 @@ PairADP::~PairADP()
if (setfl) {
for (int i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i];
delete [] setfl->elements;
- delete [] setfl->mass;
+ memory->destroy(setfl->mass);
memory->destroy(setfl->frho);
memory->destroy(setfl->rhor);
memory->destroy(setfl->z2r);
@@ -453,7 +453,7 @@ void PairADP::coeff(int narg, char **arg)
if (setfl) {
for (i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i];
delete [] setfl->elements;
- delete [] setfl->mass;
+ memory->destroy(setfl->mass);
memory->destroy(setfl->frho);
memory->destroy(setfl->rhor);
memory->destroy(setfl->z2r);
@@ -541,110 +541,130 @@ void PairADP::read_file(char *filename)
{
Setfl *file = setfl;
- // open potential file
+ // read potential file
+ if(comm->me == 0) {
+ PotentialFileReader reader(lmp, filename, "ADP");
- int me = comm->me;
- FILE *fp;
- char line[MAXLINE];
+ try {
+ char * line = nullptr;
- if (me == 0) {
- fp = force->open_potential(filename);
- if (fp == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open ADP potential file %s",filename);
- error->one(FLERR,str);
+ reader.skip_line();
+ reader.skip_line();
+ reader.skip_line();
+
+ // extract element names from nelements line
+ line = reader.next_line(1);
+ ValueTokenizer values(line);
+ file->nelements = values.next_int();
+
+ if (values.count() != file->nelements + 1)
+ error->one(FLERR,"Incorrect element names in ADP potential file");
+
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) {
+ const std::string word = values.next_string();
+ const int n = word.length() + 1;
+ file->elements[i] = new char[n];
+ strcpy(file->elements[i], word.c_str());
+ }
+
+ //
+
+ line = reader.next_line(5);
+ values = ValueTokenizer(line);
+ file->nrho = values.next_int();
+ file->drho = values.next_double();
+ file->nr = values.next_int();
+ file->dr = values.next_double();
+ file->cut = values.next_double();
+
+ if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
+ error->one(FLERR,"Invalid EAM potential file");
+
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
+ memory->create(file->u2r, file->nelements, file->nelements, file->nr + 1, "pair:u2r");
+ memory->create(file->w2r, file->nelements, file->nelements, file->nr + 1, "pair:w2r");
+
+ for (int i = 0; i < file->nelements; i++) {
+ line = reader.next_line(2);
+ values = ValueTokenizer(line);
+ values.next_int(); // ignore
+ file->mass[i] = values.next_double();
+
+ reader.next_dvector(file->nrho, &file->frho[i][1]);
+ reader.next_dvector(file->nr, &file->rhor[i][1]);
+ }
+
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ reader.next_dvector(file->nr, &file->z2r[i][j][1]);
+ }
+ }
+
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ reader.next_dvector(file->nr, &file->u2r[i][j][1]);
+ }
+ }
+
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ reader.next_dvector(file->nr, &file->w2r[i][j][1]);
+ }
+ }
+ } catch (TokenizerException & e) {
+ error->one(FLERR, e.what());
}
}
- // read and broadcast header
- // extract element names from nelements line
+ // broadcast potential information
+ MPI_Bcast(&file->nelements, 1, MPI_INT, 0, world);
- int n;
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fp,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fp,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fp,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fp,filename,error);
- n = strlen(line) + 1;
+ MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->nr, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world);
+
+ // allocate memory on other procs
+ if (comm->me != 0) {
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) file->elements[i] = nullptr;
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
+ memory->create(file->u2r, file->nelements, file->nelements, file->nr + 1, "pair:u2r");
+ memory->create(file->w2r, file->nelements, file->nelements, file->nr + 1, "pair:w2r");
}
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
- sscanf(line,"%d",&file->nelements);
- int nwords = atom->count_words(line);
- if (nwords != file->nelements + 1)
- error->all(FLERR,"Incorrect element names in ADP potential file");
-
- char **words = new char*[file->nelements+1];
- nwords = 0;
- strtok(line," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
-
- file->elements = new char*[file->nelements];
+ // broadcast file->elements string array
for (int i = 0; i < file->nelements; i++) {
- n = strlen(words[i]) + 1;
- file->elements[i] = new char[n];
- strcpy(file->elements[i],words[i]);
- }
- delete [] words;
-
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fp,filename,error);
- sscanf(line,"%d %lg %d %lg %lg",
- &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);
+ int n;
+ if (comm->me == 0) n = strlen(file->elements[i]) + 1;
+ MPI_Bcast(&n, 1, MPI_INT, 0, world);
+ if (comm->me != 0) file->elements[i] = new char[n];
+ MPI_Bcast(file->elements[i], n, MPI_CHAR, 0, world);
}
- MPI_Bcast(&file->nrho,1,MPI_INT,0,world);
- MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->nr,1,MPI_INT,0,world);
- MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world);
-
- file->mass = new double[file->nelements];
- memory->create(file->frho,file->nelements,file->nrho+1,"pair:frho");
- memory->create(file->rhor,file->nelements,file->nr+1,"pair:rhor");
- memory->create(file->z2r,file->nelements,file->nelements,file->nr+1,
- "pair:z2r");
- memory->create(file->u2r,file->nelements,file->nelements,file->nr+1,
- "pair:u2r");
- memory->create(file->w2r,file->nelements,file->nelements,file->nr+1,
- "pair:w2r");
-
- int i,j,tmp;
- for (i = 0; i < file->nelements; i++) {
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fp,filename,error);
- sscanf(line,"%d %lg",&tmp,&file->mass[i]);
- }
+ // broadcast file->mass, frho, rhor
+ for (int i = 0; i < file->nelements; i++) {
MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world);
-
- if (me == 0) grab(fp,filename,file->nrho,&file->frho[i][1]);
MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world);
- if (me == 0) grab(fp,filename,file->nr,&file->rhor[i][1]);
MPI_Bcast(&file->rhor[i][1],file->nr,MPI_DOUBLE,0,world);
}
- for (i = 0; i < file->nelements; i++)
- for (j = 0; j <= i; j++) {
- if (me == 0) grab(fp,filename,file->nr,&file->z2r[i][j][1]);
+ // broadcast file->z2r, u2r, w2r
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
- }
-
- for (i = 0; i < file->nelements; i++)
- for (j = 0; j <= i; j++) {
- if (me == 0) grab(fp,filename,file->nr,&file->u2r[i][j][1]);
MPI_Bcast(&file->u2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
- }
-
- for (i = 0; i < file->nelements; i++)
- for (j = 0; j <= i; j++) {
- if (me == 0) grab(fp,filename,file->nr,&file->w2r[i][j][1]);
MPI_Bcast(&file->w2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
}
-
- // close the potential file
-
- if (me == 0) fclose(fp);
+ }
}
/* ----------------------------------------------------------------------
@@ -912,26 +932,6 @@ void PairADP::interpolate(int n, double delta, double *f, double **spline)
}
}
-/* ----------------------------------------------------------------------
- grab n values from file fp and put them in list
- values can be several to a line
- only called by proc 0
-------------------------------------------------------------------------- */
-
-void PairADP::grab(FILE *fp, char *filename, int n, double *list)
-{
- char *ptr;
- char line[MAXLINE];
-
- int i = 0;
- while (i < n) {
- utils::sfgets(FLERR,line,MAXLINE,fp,filename,error);
- ptr = strtok(line," \t\n\r\f");
- if (ptr) list[i++] = atof(ptr);
- while ((ptr = strtok(NULL," \t\n\r\f"))) list[i++] = atof(ptr);
- }
-}
-
/* ---------------------------------------------------------------------- */
int PairADP::pack_forward_comm(int n, int *list, double *buf,
diff --git a/src/MANYBODY/pair_adp.h b/src/MANYBODY/pair_adp.h
index 360e910c2a..65596113ba 100644
--- a/src/MANYBODY/pair_adp.h
+++ b/src/MANYBODY/pair_adp.h
@@ -82,7 +82,6 @@ class PairADP : public Pair {
void allocate();
void array2spline();
void interpolate(int, double, double *, double **);
- void grab(FILE *, char *, int, double *);
void read_file(char *);
void file2array();
diff --git a/src/MANYBODY/pair_bop.cpp b/src/MANYBODY/pair_bop.cpp
index 10ab1d7080..14b0d35a4c 100644
--- a/src/MANYBODY/pair_bop.cpp
+++ b/src/MANYBODY/pair_bop.cpp
@@ -1035,145 +1035,6 @@ void PairBOP::gneigh()
/* ---------------------------------------------------------------------- */
-void PairBOP::theta()
-{
- int i,j,ii,jj,kk;
- int itype,jtype,i12;
- int temp_ij,temp_ik,temp_ijk;
- int n,nlocal,nall,ks;
- int nlisti;
- int *ilist;
- int *iilist;
- int **firstneigh;
- double rj2,rk2,rsq,ps;
- double rj1k1,rj2k2;
- double **x = atom->x;
- int *type = atom->type;
-
- nlocal = atom->nlocal;
- nall = nlocal+atom->nghost;
- ilist = list->ilist;
- firstneigh = list->firstneigh;
- if(update_list!=0)
- memory_theta_grow();
- else
- memory_theta_create();
- for (ii = 0; ii < nall; ii++) {
- if(ii=npairs) {
- error->one(FLERR,"Too many atom pairs for pair bop");
- }
- disij[0][temp_ij]=x[j][0]-x[i][0];
- disij[1][temp_ij]=x[j][1]-x[i][1];
- disij[2][temp_ij]=x[j][2]-x[i][2];
- rsq=disij[0][temp_ij]*disij[0][temp_ij]
- +disij[1][temp_ij]*disij[1][temp_ij]
- +disij[2][temp_ij]*disij[2][temp_ij];
- rij[temp_ij]=sqrt(rsq);
- if(rij[temp_ij]<=rcut[i12])
- neigh_flag[temp_ij]=1;
- else
- neigh_flag[temp_ij]=0;
- if(rij[temp_ij]<=rcut3[i12])
- neigh_flag3[temp_ij]=1;
- else
- neigh_flag3[temp_ij]=0;
- ps=rij[temp_ij]*rdr[i12]+1.0;
- ks=(int)ps;
- if(nr-11.0)
- ps=1.0;
- betaS[temp_ij]=((pBetaS3[i12][ks-1]*ps+pBetaS2[i12][ks-1])*ps+pBetaS1[i12][ks-1])*ps+pBetaS[i12][ks-1];
- dBetaS[temp_ij]=(pBetaS6[i12][ks-1]*ps+pBetaS5[i12][ks-1])*ps
- +pBetaS4[i12][ks-1];
- betaP[temp_ij]=((pBetaP3[i12][ks-1]*ps+pBetaP2[i12][ks-1])*ps
- +pBetaP1[i12][ks-1])*ps+pBetaP[i12][ks-1];
- dBetaP[temp_ij]=(pBetaP6[i12][ks-1]*ps+pBetaP5[i12][ks-1])*ps
- +pBetaP4[i12][ks-1];
- repul[temp_ij]=((pRepul3[i12][ks-1]*ps+pRepul2[i12][ks-1])*ps
- +pRepul1[i12][ks-1])*ps+pRepul[i12][ks-1];
- dRepul[temp_ij]=(pRepul6[i12][ks-1]*ps+pRepul5[i12][ks-1])*ps
- +pRepul4[i12][ks-1];
- }
- }
- for (ii = 0; ii < nall; ii++) {
- n=0;
- if(ii=cos_total) {
- error->one(FLERR,"Too many atom triplets for pair bop");
- }
- temp_ik=BOP_index[i]+kk;
- temp_ijk=cos_index[i]+n;
- if(temp_ijk>=cos_total) {
- error->one(FLERR,"Too many atom triplets for pair bop");
- }
- rk2=rij[temp_ik]*rij[temp_ik];
- rj1k1=rij[temp_ij]*rij[temp_ik];
- rj2k2=rj1k1*rj1k1;
- if(temp_ijk>=cos_total) {
- error->one(FLERR,"Too many atom triplets for pair bop");
- }
- cosAng[temp_ijk]=(disij[0][temp_ij]*disij[0][temp_ik]+disij[1][temp_ij]
- *disij[1][temp_ik]+disij[2][temp_ij]*disij[2][temp_ik])/rj1k1;
- dcAng[temp_ijk][0][0]=(disij[0][temp_ik]*rj1k1-cosAng[temp_ijk]
- *disij[0][temp_ij]*rk2)/(rj2k2);
- dcAng[temp_ijk][1][0]=(disij[1][temp_ik]*rj1k1-cosAng[temp_ijk]
- *disij[1][temp_ij]*rk2)/(rj2k2);
- dcAng[temp_ijk][2][0]=(disij[2][temp_ik]*rj1k1-cosAng[temp_ijk]
- *disij[2][temp_ij]*rk2)/(rj2k2);
- dcAng[temp_ijk][0][1]=(disij[0][temp_ij]*rj1k1-cosAng[temp_ijk]
- *disij[0][temp_ik]*rj2)/(rj2k2);
- dcAng[temp_ijk][1][1]=(disij[1][temp_ij]*rj1k1-cosAng[temp_ijk]
- *disij[1][temp_ik]*rj2)/(rj2k2);
- dcAng[temp_ijk][2][1]=(disij[2][temp_ij]*rj1k1-cosAng[temp_ijk]
- *disij[2][temp_ik]*rj2)/(rj2k2);
- n++;
- }
- }
- }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairBOP::theta_mod()
-{
- if(update_list!=0)
- memory_theta_grow();
- else
- memory_theta_create();
-}
-
-/* ---------------------------------------------------------------------- */
-
/* The formulation differs slightly to avoid negative square roots
in the calculation of Sigma^(1/2) of (a) Eq. 6 and (b) Eq. 11 */
@@ -5322,11 +5183,13 @@ void _noopt PairBOP::read_table(char *filename)
if(rcut3[i]>rcutall)
rcutall=rcut3[i];
rcutsq[i]=rcut[i]*rcut[i];
- rcutsq3[i]=rcut3[i]*rcut3[i];
dr[i]=rcut[i]/((double)nr-1.0);
rdr[i]=1.0/dr[i];
- dr3[i]=rcut3[i]/((double)nr-1.0);
- rdr3[i]=1.0/dr3[i];
+ if (nws==3) {
+ rcutsq3[i]=rcut3[i]*rcut3[i];
+ dr3[i]=rcut3[i]/((double)nr-1.0);
+ rdr3[i]=1.0/dr3[i];
+ }
}
rctroot=rcutall;
dtheta=2.0/((double)ntheta-1.0);
@@ -5412,9 +5275,11 @@ void _noopt PairBOP::read_table(char *filename)
pRepul4[i][k]=pRepul1[i][k]/dr[i];
pRepul5[i][k]=2.0*pRepul2[i][k]/dr[i];
pRepul6[i][k]=3.0*pRepul3[i][k]/dr[i];
- pLong4[i][k]=pLong1[i][k]/dr3[i];
- pLong5[i][k]=2.0*pLong2[i][k]/dr3[i];
- pLong6[i][k]=3.0*pLong3[i][k]/dr3[i];
+ if (nws==3) {
+ pLong4[i][k]=pLong1[i][k]/dr3[i];
+ pLong5[i][k]=2.0*pLong2[i][k]/dr3[i];
+ pLong6[i][k]=3.0*pLong3[i][k]/dr3[i];
+ }
}
for(k=0;ksfree(params);
- params = NULL;
+ params = nullptr;
nparams = 0;
maxparam = 0;
// open file on proc 0
-
- FILE *fp;
if (comm->me == 0) {
- fp = force->open_potential(file);
- if (fp == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open COMB potential file %s",file);
- error->one(FLERR,str);
- }
- }
+ PotentialFileReader reader(lmp, file, "COMB");
+ char * line;
- // read each line out of file, skipping blank lines or leading '#'
- // store line of params if all 3 element tags are in element list
+ while((line = reader.next_line(NPARAMS_PER_LINE))) {
+ try {
+ ValueTokenizer values(line," \t\n\r\f");
- int n,nwords,ielement,jelement,kelement;
- char line[MAXLINE],*ptr;
- int eof = 0;
+ std::string iname = values.next_string();
+ std::string jname = values.next_string();
+ std::string kname = values.next_string();
- while (1) {
- if (comm->me == 0) {
- ptr = fgets(line,MAXLINE,fp);
- if (ptr == NULL) {
- eof = 1;
- fclose(fp);
- } else n = strlen(line) + 1;
- }
- MPI_Bcast(&eof,1,MPI_INT,0,world);
- if (eof) break;
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
+ // ielement,jelement,kelement = 1st args
+ // if all 3 args are in element list, then parse this line
+ // else skip to next line
+ int ielement, jelement, kelement;
- // strip comment, skip line if blank
+ for (ielement = 0; ielement < nelements; ielement++)
+ if (iname == elements[ielement]) break;
+ if (ielement == nelements) continue;
+ for (jelement = 0; jelement < nelements; jelement++)
+ if (jname == elements[jelement]) break;
+ if (jelement == nelements) continue;
+ for (kelement = 0; kelement < nelements; kelement++)
+ if (kname == elements[kelement]) break;
+ if (kelement == nelements) continue;
- if ((ptr = strchr(line,'#'))) *ptr = '\0';
- nwords = atom->count_words(line);
- if (nwords == 0) continue;
+ // load up parameter settings and error check their values
- // concatenate additional lines until have params_per_line words
+ if (nparams == maxparam) {
+ maxparam += DELTA;
+ params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
+ "pair:params");
+ }
- while (nwords < params_per_line) {
- n = strlen(line);
- if (comm->me == 0) {
- ptr = fgets(&line[n],MAXLINE-n,fp);
- if (ptr == NULL) {
- eof = 1;
- fclose(fp);
- } else n = strlen(line) + 1;
+ params[nparams].ielement = ielement;
+ params[nparams].jelement = jelement;
+ params[nparams].kelement = kelement;
+ params[nparams].powerm = values.next_double();
+ params[nparams].c = values.next_double();
+ params[nparams].d = values.next_double();
+ params[nparams].h = values.next_double();
+ params[nparams].powern = values.next_double();
+ params[nparams].beta = values.next_double();
+ params[nparams].lam21 = values.next_double();
+ params[nparams].lam22 = values.next_double();
+ params[nparams].bigb1 = values.next_double();
+ params[nparams].bigb2 = values.next_double();
+ params[nparams].bigr = values.next_double();
+ params[nparams].bigd = values.next_double();
+ params[nparams].lam11 = values.next_double();
+ params[nparams].lam12 = values.next_double();
+ params[nparams].biga1 = values.next_double();
+ params[nparams].biga2 = values.next_double();
+ params[nparams].plp1 = values.next_double();
+ params[nparams].plp3 = values.next_double();
+ params[nparams].plp6 = values.next_double();
+ params[nparams].a123 = values.next_double();
+ params[nparams].aconf = values.next_double();
+ params[nparams].addrep = values.next_double();
+ params[nparams].romigb = values.next_double();
+ params[nparams].romigc = values.next_double();
+ params[nparams].romigd = values.next_double();
+ params[nparams].romiga = values.next_double();
+ params[nparams].QL1 = values.next_double();
+ params[nparams].QU1 = values.next_double();
+ params[nparams].DL1 = values.next_double();
+ params[nparams].DU1 = values.next_double();
+ params[nparams].QL2 = values.next_double();
+ params[nparams].QU2 = values.next_double();
+ params[nparams].DL2 = values.next_double();
+ params[nparams].DU2 = values.next_double();
+ params[nparams].chi = values.next_double();
+ params[nparams].dj = values.next_double();
+ params[nparams].dk = values.next_double();
+ params[nparams].dl = values.next_double();
+ params[nparams].dm = values.next_double();
+ params[nparams].esm1 = values.next_double();
+ params[nparams].cmn1 = values.next_double();
+ params[nparams].cml1 = values.next_double();
+ params[nparams].cmn2 = values.next_double();
+ params[nparams].cml2 = values.next_double();
+ params[nparams].coulcut = values.next_double();
+ params[nparams].hfocor = values.next_double();
+ } catch (TokenizerException & e) {
+ error->one(FLERR, e.what());
}
- MPI_Bcast(&eof,1,MPI_INT,0,world);
- if (eof) break;
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
- if ((ptr = strchr(line,'#'))) *ptr = '\0';
- nwords = atom->count_words(line);
+
+ params[nparams].powermint = int(params[nparams].powerm);
+
+ // parameter sanity checks
+
+ if (params[nparams].lam11 < 0.0 || params[nparams].lam12 < 0.0 ||
+ params[nparams].c < 0.0 || params[nparams].d < 0.0 ||
+ params[nparams].powern < 0.0 || params[nparams].beta < 0.0 ||
+ params[nparams].lam21 < 0.0 || params[nparams].lam22 < 0.0 ||
+ params[nparams].bigb1< 0.0 || params[nparams].bigb2< 0.0 ||
+ params[nparams].biga1< 0.0 || params[nparams].biga2< 0.0 ||
+ params[nparams].bigr < 0.0 || params[nparams].bigd < 0.0 ||
+ params[nparams].bigd > params[nparams].bigr ||
+ params[nparams].powerm - params[nparams].powermint != 0.0 ||
+ (params[nparams].powermint != 3 && params[nparams].powermint != 1) ||
+ params[nparams].plp1 < 0.0 || params[nparams].plp3 < 0.0 ||
+ params[nparams].plp6 < 0.0 ||
+ params[nparams].a123 > 360.0 || params[nparams].aconf < 0.0 ||
+ params[nparams].addrep < 0.0 || params[nparams].romigb < 0.0 ||
+ params[nparams].romigc < 0.0 || params[nparams].romigd < 0.0 ||
+ params[nparams].romiga < 0.0 ||
+ params[nparams].QL1 > 0.0 || params[nparams].QU1 < 0.0 ||
+ params[nparams].DL1 < 0.0 || params[nparams].DU1 > 0.0 ||
+ params[nparams].QL2 > 0.0 || params[nparams].QU2 < 0.0 ||
+ params[nparams].DL2 < 0.0 || params[nparams].DU2 > 0.0 ||
+ params[nparams].chi < 0.0 ||
+ // params[nparams].dj < 0.0 || params[nparams].dk < 0.0 ||
+ // params[nparams].dl < 0.0 || params[nparams].dm < 0.0 ||
+ params[nparams].esm1 < 0.0)
+ error->one(FLERR,"Illegal COMB parameter");
+
+ if (params[nparams].lam11 < params[nparams].lam21 ||
+ params[nparams].lam12 < params[nparams].lam22 ||
+ params[nparams].biga1< params[nparams].bigb1 ||
+ params[nparams].biga2< params[nparams].bigb2)
+ error->one(FLERR,"Illegal COMB parameter");
+
+ nparams++;
}
-
- if (nwords != params_per_line)
- error->all(FLERR,"Incorrect format in COMB potential file");
-
- // words = ptrs to all words in line
-
- nwords = 0;
- words[nwords++] = strtok(line," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
-
- // ielement,jelement,kelement = 1st args
- // if all 3 args are in element list, then parse this line
- // else skip to next line
-
- for (ielement = 0; ielement < nelements; ielement++)
- if (strcmp(words[0],elements[ielement]) == 0) break;
- if (ielement == nelements) continue;
- for (jelement = 0; jelement < nelements; jelement++)
- if (strcmp(words[1],elements[jelement]) == 0) break;
- if (jelement == nelements) continue;
- for (kelement = 0; kelement < nelements; kelement++)
- if (strcmp(words[2],elements[kelement]) == 0) break;
- if (kelement == nelements) continue;
-
- // load up parameter settings and error check their values
-
- if (nparams == maxparam) {
- maxparam += DELTA;
- params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
- "pair:params");
- }
-
- params[nparams].ielement = ielement;
- params[nparams].jelement = jelement;
- params[nparams].kelement = kelement;
- params[nparams].powerm = atof(words[3]);
- params[nparams].c = atof(words[4]);
- params[nparams].d = atof(words[5]);
- params[nparams].h = atof(words[6]);
- params[nparams].powern = atof(words[7]);
- params[nparams].beta = atof(words[8]);
- params[nparams].lam21 = atof(words[9]);
- params[nparams].lam22 = atof(words[10]);
- params[nparams].bigb1 = atof(words[11]);
- params[nparams].bigb2 = atof(words[12]);
- params[nparams].bigr = atof(words[13]);
- params[nparams].bigd = atof(words[14]);
- params[nparams].lam11 = atof(words[15]);
- params[nparams].lam12 = atof(words[16]);
- params[nparams].biga1 = atof(words[17]);
- params[nparams].biga2 = atof(words[18]);
- params[nparams].plp1 = atof(words[19]);
- params[nparams].plp3 = atof(words[20]);
- params[nparams].plp6 = atof(words[21]);
- params[nparams].a123 = atof(words[22]);
- params[nparams].aconf= atof(words[23]);
- params[nparams].addrep = atof(words[24]);
- params[nparams].romigb = atof(words[25]);
- params[nparams].romigc = atof(words[26]);
- params[nparams].romigd = atof(words[27]);
- params[nparams].romiga = atof(words[28]);
- params[nparams].QL1 = atof(words[29]);
- params[nparams].QU1 = atof(words[30]);
- params[nparams].DL1 = atof(words[31]);
- params[nparams].DU1 = atof(words[32]);
- params[nparams].QL2 = atof(words[33]);
- params[nparams].QU2 = atof(words[34]);
- params[nparams].DL2 = atof(words[35]);
- params[nparams].DU2 = atof(words[36]);
- params[nparams].chi = atof(words[37]);
- params[nparams].dj = atof(words[38]);
- params[nparams].dk = atof(words[39]);
- params[nparams].dl = atof(words[40]);
- params[nparams].dm = atof(words[41]);
- params[nparams].esm1 = atof(words[42]);
- params[nparams].cmn1 = atof(words[43]);
- params[nparams].cml1 = atof(words[44]);
- params[nparams].cmn2 = atof(words[45]);
- params[nparams].cml2 = atof(words[46]);
- params[nparams].coulcut = atof(words[47]);
- params[nparams].hfocor = atof(words[48]);
-
- params[nparams].powermint = int(params[nparams].powerm);
-
- // parameter sanity checks
-
- if (params[nparams].lam11 < 0.0 || params[nparams].lam12 < 0.0 ||
- params[nparams].c < 0.0 || params[nparams].d < 0.0 ||
- params[nparams].powern < 0.0 || params[nparams].beta < 0.0 ||
- params[nparams].lam21 < 0.0 || params[nparams].lam22 < 0.0 ||
- params[nparams].bigb1< 0.0 || params[nparams].bigb2< 0.0 ||
- params[nparams].biga1< 0.0 || params[nparams].biga2< 0.0 ||
- params[nparams].bigr < 0.0 || params[nparams].bigd < 0.0 ||
- params[nparams].bigd > params[nparams].bigr ||
- params[nparams].powerm - params[nparams].powermint != 0.0 ||
- (params[nparams].powermint != 3 && params[nparams].powermint != 1) ||
- params[nparams].plp1 < 0.0 || params[nparams].plp3 < 0.0 ||
- params[nparams].plp6 < 0.0 ||
- params[nparams].a123 > 360.0 || params[nparams].aconf < 0.0 ||
- params[nparams].addrep < 0.0 || params[nparams].romigb < 0.0 ||
- params[nparams].romigc < 0.0 || params[nparams].romigd < 0.0 ||
- params[nparams].romiga < 0.0 ||
- params[nparams].QL1 > 0.0 || params[nparams].QU1 < 0.0 ||
- params[nparams].DL1 < 0.0 || params[nparams].DU1 > 0.0 ||
- params[nparams].QL2 > 0.0 || params[nparams].QU2 < 0.0 ||
- params[nparams].DL2 < 0.0 || params[nparams].DU2 > 0.0 ||
- params[nparams].chi < 0.0 ||
-// params[nparams].dj < 0.0 || params[nparams].dk < 0.0 ||
-// params[nparams].dl < 0.0 || params[nparams].dm < 0.0 ||
- params[nparams].esm1 < 0.0)
- error->all(FLERR,"Illegal COMB parameter");
-
- if (params[nparams].lam11 < params[nparams].lam21 ||
- params[nparams].lam12 < params[nparams].lam22 ||
- params[nparams].biga1< params[nparams].bigb1 ||
- params[nparams].biga2< params[nparams].bigb2)
- error->all(FLERR,"Illegal COMB parameter");
-
- nparams++;
}
- delete [] words;
+ MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
+ MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
+
+ if(comm->me != 0) {
+ params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
+ }
+
+ MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/MANYBODY/pair_comb.h b/src/MANYBODY/pair_comb.h
index 19788dfab1..7a3d279033 100644
--- a/src/MANYBODY/pair_comb.h
+++ b/src/MANYBODY/pair_comb.h
@@ -38,6 +38,8 @@ class PairComb : public Pair {
virtual double yasu_char(double *, int &);
double enegtot;
+ static const int NPARAMS_PER_LINE = 49;
+
protected:
struct Param {
double lam11,lam12,lam21,lam22;
diff --git a/src/MANYBODY/pair_comb3.cpp b/src/MANYBODY/pair_comb3.cpp
index a79d919d64..8a010df55d 100644
--- a/src/MANYBODY/pair_comb3.cpp
+++ b/src/MANYBODY/pair_comb3.cpp
@@ -34,11 +34,12 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
using namespace LAMMPS_NS;
using namespace MathConst;
-#define MAXLINE 1024
#define DELTA 4
#define PGDELTA 1
#define MAXNEIGH 24
@@ -312,220 +313,211 @@ double PairComb3::init_one(int i, int j)
void PairComb3::read_lib()
{
const unsigned int MAXLIB = 1024;
- int i,j,k,l,nwords,m;
+ int i,j,k,l,m;
int ii,jj,kk,ll,mm,iii;
char s[MAXLIB];
- char **words = new char*[80];
- // open libraray file on proc 0
+ // open library file on proc 0
if (comm->me == 0) {
const char filename[] = "lib.comb3";
FILE *fp = force->open_potential(filename);
if (fp == NULL) error->one(FLERR,"Cannot open COMB3 lib.comb3 file");
- // read and store at the same time
- utils::sfgets(FLERR,s,MAXLIB,fp,filename,error);
- utils::sfgets(FLERR,s,MAXLIB,fp,filename,error);
- nwords = 0;
- words[nwords++] = strtok(s," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
- ccutoff[0] = atof(words[0]);
- ccutoff[1] = atof(words[1]);
- ccutoff[2] = atof(words[2]);
- ccutoff[3] = atof(words[3]);
- ccutoff[4] = atof(words[4]);
- ccutoff[5] = atof(words[5]);
+ try {
+ // read and store at the same time
+ utils::sfgets(FLERR, s, MAXLIB, fp, filename, error);
- utils::sfgets(FLERR,s,MAXLIB,fp,filename,error);
- nwords = 0;
- words[nwords++] = strtok(s," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
- ch_a[0] = atof(words[0]);
- ch_a[1] = atof(words[1]);
- ch_a[2] = atof(words[2]);
- ch_a[3] = atof(words[3]);
- ch_a[4] = atof(words[4]);
- ch_a[5] = atof(words[5]);
- ch_a[6] = atof(words[6]);
+ utils::sfgets(FLERR, s, MAXLIB, fp, filename, error);
+ ValueTokenizer values(s, " \t\n\r\f");
- utils::sfgets(FLERR,s,MAXLIB,fp,filename,error);
- nwords = 0;
- words[nwords++] = strtok(s," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
- nsplpcn = atoi(words[0]);
- nsplrad = atoi(words[1]);
- nspltor = atoi(words[2]);
+ ccutoff[0] = values.next_double();
+ ccutoff[1] = values.next_double();
+ ccutoff[2] = values.next_double();
+ ccutoff[3] = values.next_double();
+ ccutoff[4] = values.next_double();
+ ccutoff[5] = values.next_double();
- utils::sfgets(FLERR,s,MAXLIB,fp,filename,error);
- nwords = 0;
- words[nwords++] = strtok(s," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
- maxx = atoi(words[0]);
- maxy = atoi(words[1]);
- maxz = atoi(words[2]);
+ utils::sfgets(FLERR, s, MAXLIB, fp, filename, error);
+ values = ValueTokenizer(s, " \t\n\r\f");
- utils::sfgets(FLERR,s,MAXLIB,fp,filename,error);
- nwords = 0;
- words[nwords++] = strtok(s," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f")))continue;
- maxxc = atoi(words[0]);
- maxyc = atoi(words[1]);
- maxconj = atoi(words[2]);
+ ch_a[0] = values.next_double();
+ ch_a[1] = values.next_double();
+ ch_a[2] = values.next_double();
+ ch_a[3] = values.next_double();
+ ch_a[4] = values.next_double();
+ ch_a[5] = values.next_double();
+ ch_a[6] = values.next_double();
+
+ utils::sfgets(FLERR, s, MAXLIB, fp, filename, error);
+ values = ValueTokenizer(s, " \t\n\r\f");
+
+ nsplpcn = values.next_int();
+ nsplrad = values.next_int();
+ nspltor = values.next_int();
+
+ utils::sfgets(FLERR, s, MAXLIB, fp, filename, error);
+ values = ValueTokenizer(s, " \t\n\r\f");
+
+ maxx = values.next_int();
+ maxy = values.next_int();
+ maxz = values.next_int();
- for (l=0; lone(FLERR, e.what());
+ }
}
k = 0;
@@ -584,211 +576,163 @@ void PairComb3::read_lib()
MPI_Bcast(&iin2[0][0],32,MPI_INT,0,world);
MPI_Bcast(&iin3[0][0],192,MPI_INT,0,world);
- delete [] words;
}
/* ---------------------------------------------------------------------- */
void PairComb3::read_file(char *file)
{
- int params_per_line = 74;
- char **words = new char*[params_per_line+1];
-
memory->sfree(params);
- params = NULL;
- nparams = maxparam = 0;
+ params = nullptr;
+ nparams = 0;
+ maxparam = 0;
// open file on proc 0
-
- FILE *fp;
if (comm->me == 0) {
- fp = force->open_potential(file);
- if (fp == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open COMB3 potential file %s",file);
- error->one(FLERR,str);
- }
- }
+ PotentialFileReader reader(lmp, file, "COMB3");
+ char * line;
- // read each line out of file, skipping blank lines or leading '#'
- // store line of params if all 3 element tags are in element list
+ while((line = reader.next_line(NPARAMS_PER_LINE))) {
+ try {
+ ValueTokenizer values(line, " \t\n\r\f");
- int n,nwords,ielement,jelement,kelement;
- char line[MAXLINE],*ptr;
- int eof = 0;
- nwords=0;
- while (1) {
- if (comm->me == 0) {
- ptr = fgets(line,MAXLINE,fp);
- if (ptr == NULL) {
- eof = 1;
- fclose(fp);
- } else n = strlen(line) + 1;
- }
- MPI_Bcast(&eof,1,MPI_INT,0,world);
- if (eof) break;
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
+ std::string iname = values.next_string();
+ std::string jname = values.next_string();
+ std::string kname = values.next_string();
- // strip comment, skip line if blank
+ // ielement,jelement,kelement = 1st args
+ // if all 3 args are in element list, then parse this line
+ // else skip to next line
+ int ielement, jelement, kelement;
- if ((ptr = strchr(line,'#'))) *ptr = '\0';
+ for (ielement = 0; ielement < nelements; ielement++)
+ if (iname == elements[ielement]) break;
+ if (ielement == nelements) continue;
+ for (jelement = 0; jelement < nelements; jelement++)
+ if (jname == elements[jelement]) break;
+ if (jelement == nelements) continue;
+ for (kelement = 0; kelement < nelements; kelement++)
+ if (kname == elements[kelement]) break;
+ if (kelement == nelements) continue;
- nwords = atom->count_words(line);
- if (nwords == 0) continue;
+ // load up parameter settings and error check their values
- // concatenate additional lines until have params_per_line words
+ if (nparams == maxparam) {
+ maxparam += DELTA;
+ params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
+ "pair:params");
+ }
- while (nwords < params_per_line) {
- n = strlen(line);
-
- if (comm->me == 0) {
- ptr = fgets(&line[n],MAXLINE-n,fp);
- if (ptr == NULL) {
- eof = 1;
- fclose(fp);
- } else n = strlen(line) + 1;
+ params[nparams].ielement = ielement;
+ params[nparams].jelement = jelement;
+ params[nparams].kelement = kelement;
+ params[nparams].ielementgp = values.next_int();
+ params[nparams].jelementgp = values.next_int();
+ params[nparams].kelementgp = values.next_int();
+ params[nparams].ang_flag = values.next_int();
+ params[nparams].pcn_flag = values.next_int();
+ params[nparams].rad_flag = values.next_int();
+ params[nparams].tor_flag = values.next_int();
+ params[nparams].vdwflag = values.next_double();
+ params[nparams].powerm = values.next_double();
+ params[nparams].veps = values.next_double();
+ params[nparams].vsig = values.next_double();
+ params[nparams].paaa = values.next_double();
+ params[nparams].pbbb = values.next_double();
+ params[nparams].lami = values.next_double();
+ params[nparams].alfi = values.next_double();
+ params[nparams].powern = values.next_double();
+ params[nparams].QL = values.next_double();
+ params[nparams].QU = values.next_double();
+ params[nparams].DL = values.next_double();
+ params[nparams].DU = values.next_double();
+ params[nparams].qmin = values.next_double();
+ params[nparams].qmax = values.next_double();
+ params[nparams].chi = values.next_double();
+ params[nparams].dj = values.next_double();
+ params[nparams].dk = values.next_double();
+ params[nparams].dl = values.next_double();
+ params[nparams].esm = values.next_double();
+ params[nparams].cmn1 = values.next_double();
+ params[nparams].cmn2 = values.next_double();
+ params[nparams].pcmn1 = values.next_double();
+ params[nparams].pcmn2 = values.next_double();
+ params[nparams].coulcut = values.next_double();
+ params[nparams].polz = values.next_double();
+ params[nparams].curl = values.next_double();
+ params[nparams].curlcut1 = values.next_double();
+ params[nparams].curlcut2 = values.next_double();
+ params[nparams].curl0 = values.next_double();
+ params[nparams].alpha1 = values.next_double();
+ params[nparams].bigB1 = values.next_double();
+ params[nparams].alpha2 = values.next_double();
+ params[nparams].bigB2 = values.next_double();
+ params[nparams].alpha3 = values.next_double();
+ params[nparams].bigB3 = values.next_double();
+ params[nparams].lambda = values.next_double();
+ params[nparams].bigA = values.next_double();
+ params[nparams].beta = values.next_double();
+ params[nparams].bigr = values.next_double();
+ params[nparams].bigd = values.next_double();
+ params[nparams].pcos6 = values.next_double();
+ params[nparams].pcos5 = values.next_double();
+ params[nparams].pcos4 = values.next_double();
+ params[nparams].pcos3 = values.next_double();
+ params[nparams].pcos2 = values.next_double();
+ params[nparams].pcos1 = values.next_double();
+ params[nparams].pcos0 = values.next_double();
+ params[nparams].pcna = values.next_double();
+ params[nparams].pcnb = values.next_double();
+ params[nparams].pcnc = values.next_double();
+ params[nparams].pcnd = values.next_double();
+ params[nparams].p6p0 = values.next_double();
+ params[nparams].p6p1 = values.next_double();
+ params[nparams].p6p2 = values.next_double();
+ params[nparams].p6p3 = values.next_double();
+ params[nparams].p6p4 = values.next_double();
+ params[nparams].p6p5 = values.next_double();
+ params[nparams].p6p6 = values.next_double();
+ params[nparams].ptork1 = values.next_double();
+ params[nparams].ptork2 = values.next_double();
+ params[nparams].addrepr = values.next_double();
+ params[nparams].addrep = values.next_double();
+ params[nparams].pcross = values.next_double();
+ params[nparams].powermint = int(params[nparams].powerm);
+ } catch (TokenizerException & e) {
+ error->one(FLERR, e.what());
}
- MPI_Bcast(&eof,1,MPI_INT,0,world);
- if (eof) break;
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
- if ((ptr = strchr(line,'#'))) *ptr = '\0';
- nwords = atom->count_words(line);
+
+ // parameter sanity checks
+
+ if (params[nparams].lambda < 0.0 || params[nparams].powern < 0.0 ||
+ params[nparams].beta < 0.0 || params[nparams].alpha1 < 0.0 ||
+ params[nparams].bigB1< 0.0 || params[nparams].bigA< 0.0 ||
+ params[nparams].bigB2< 0.0 || params[nparams].alpha2 <0.0 ||
+ params[nparams].bigB3< 0.0 || params[nparams].alpha3 <0.0 ||
+ params[nparams].bigr < 0.0 || params[nparams].bigd < 0.0 ||
+ params[nparams].bigd > params[nparams].bigr ||
+ params[nparams].powerm - params[nparams].powermint != 0.0 ||
+ params[nparams].addrepr < 0.0 || params[nparams].powermint < 1.0 ||
+ params[nparams].QL > 0.0 || params[nparams].QU < 0.0 ||
+ params[nparams].DL < 0.0 || params[nparams].DU > 0.0 ||
+ params[nparams].pcross < 0.0 ||
+ params[nparams].esm < 0.0 || params[nparams].veps < 0.0 ||
+ params[nparams].vsig < 0.0 || params[nparams].vdwflag < 0.0
+ )
+ error->one(FLERR,"Illegal COMB3 parameter");
+
+ nparams++;
}
- if (nwords != params_per_line) {
- error->all(FLERR,"Incorrect format in COMB3 potential file");
- }
- // words = ptrs to all words in line
-
- nwords = 0;
- words[nwords++] = strtok(line," \t\n\r\f");
- while ((nwords <= params_per_line)
- && (words[nwords++] = strtok(NULL," \t\n\r\f"))) {
- continue;
- }
-
- // ielement,jelement,kelement = 1st args
- // if all 3 args are in element list, then parse this line
- // else skip to next line
-
- for (ielement = 0; ielement < nelements; ielement++)
- if (strcmp(words[0],elements[ielement]) == 0) break;
- if (ielement == nelements) continue;
- for (jelement = 0; jelement < nelements; jelement++)
- if (strcmp(words[1],elements[jelement]) == 0) break;
- if (jelement == nelements) continue;
- for (kelement = 0; kelement < nelements; kelement++)
- if (strcmp(words[2],elements[kelement]) == 0) break;
- if (kelement == nelements) continue;
-
- // load up parameter settings and error check their values
-
- if (nparams == maxparam) {
- maxparam += DELTA;
- params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
- "pair:params");
- }
-
- params[nparams].ielement = ielement;
- params[nparams].jelement = jelement;
- params[nparams].kelement = kelement;
- params[nparams].ielementgp = atoi(words[3]);
- params[nparams].jelementgp = atoi(words[4]);
- params[nparams].kelementgp = atoi(words[5]);
- params[nparams].ang_flag = atoi(words[6]);
- params[nparams].pcn_flag = atoi(words[7]);
- params[nparams].rad_flag = atoi(words[8]);
- params[nparams].tor_flag = atoi(words[9]);
- params[nparams].vdwflag = atof(words[10]);
- params[nparams].powerm = atof(words[11]);
- params[nparams].veps = atof(words[12]);
- params[nparams].vsig = atof(words[13]);
- params[nparams].paaa = atof(words[14]);
- params[nparams].pbbb = atof(words[15]);
- params[nparams].lami = atof(words[16]);
- params[nparams].alfi = atof(words[17]);
- params[nparams].powern = atof(words[18]);
- params[nparams].QL = atof(words[19]);
- params[nparams].QU = atof(words[20]);
- params[nparams].DL = atof(words[21]);
- params[nparams].DU = atof(words[22]);
- params[nparams].qmin = atof(words[23]);
- params[nparams].qmax = atof(words[24]);
- params[nparams].chi = atof(words[25]);
- params[nparams].dj = atof(words[26]);
- params[nparams].dk = atof(words[27]);
- params[nparams].dl = atof(words[28]);
- params[nparams].esm = atof(words[29]);
- params[nparams].cmn1 = atof(words[30]);
- params[nparams].cmn2 = atof(words[31]);
- params[nparams].pcmn1 = atof(words[32]);
- params[nparams].pcmn2 = atof(words[33]);
- params[nparams].coulcut = atof(words[34]);
- params[nparams].polz = atof(words[35]);
- params[nparams].curl = atof(words[36]);
- params[nparams].curlcut1 = atof(words[37]);
- params[nparams].curlcut2 = atof(words[38]);
- params[nparams].curl0 = atof(words[39]);
- params[nparams].alpha1 = atof(words[40]);
- params[nparams].bigB1 = atof(words[41]);
- params[nparams].alpha2 = atof(words[42]);
- params[nparams].bigB2 = atof(words[43]);
- params[nparams].alpha3 = atof(words[44]);
- params[nparams].bigB3 = atof(words[45]);
- params[nparams].lambda = atof(words[46]);
- params[nparams].bigA = atof(words[47]);
- params[nparams].beta = atof(words[48]);
- params[nparams].bigr = atof(words[49]);
- params[nparams].bigd = atof(words[50]);
- params[nparams].pcos6 = atof(words[51]);
- params[nparams].pcos5 = atof(words[52]);
- params[nparams].pcos4 = atof(words[53]);
- params[nparams].pcos3 = atof(words[54]);
- params[nparams].pcos2 = atof(words[55]);
- params[nparams].pcos1 = atof(words[56]);
- params[nparams].pcos0 = atof(words[57]);
- params[nparams].pcna = atof(words[58]);
- params[nparams].pcnb = atof(words[59]);
- params[nparams].pcnc = atof(words[60]);
- params[nparams].pcnd = atof(words[61]);
- params[nparams].p6p0 = atof(words[62]);
- params[nparams].p6p1 = atof(words[63]);
- params[nparams].p6p2 = atof(words[64]);
- params[nparams].p6p3 = atof(words[65]);
- params[nparams].p6p4 = atof(words[66]);
- params[nparams].p6p5 = atof(words[67]);
- params[nparams].p6p6 = atof(words[68]);
- params[nparams].ptork1=atof(words[69]);
- params[nparams].ptork2=atof(words[70]);
- params[nparams].addrepr=atof(words[71]);
- params[nparams].addrep=atof(words[72]);
- params[nparams].pcross = atof(words[73]);
- params[nparams].powermint = int(params[nparams].powerm);
-
- // parameter sanity checks
-
- if (params[nparams].lambda < 0.0 || params[nparams].powern < 0.0 ||
- params[nparams].beta < 0.0 || params[nparams].alpha1 < 0.0 ||
- params[nparams].bigB1< 0.0 || params[nparams].bigA< 0.0 ||
- params[nparams].bigB2< 0.0 || params[nparams].alpha2 <0.0 ||
- params[nparams].bigB3< 0.0 || params[nparams].alpha3 <0.0 ||
- params[nparams].bigr < 0.0 || params[nparams].bigd < 0.0 ||
- params[nparams].bigd > params[nparams].bigr ||
- params[nparams].powerm - params[nparams].powermint != 0.0 ||
- params[nparams].addrepr < 0.0 || params[nparams].powermint < 1.0 ||
- params[nparams].QL > 0.0 || params[nparams].QU < 0.0 ||
- params[nparams].DL < 0.0 || params[nparams].DU > 0.0 ||
- params[nparams].pcross < 0.0 ||
- params[nparams].esm < 0.0 || params[nparams].veps < 0.0 ||
- params[nparams].vsig < 0.0 || params[nparams].vdwflag < 0.0
- )
- error->all(FLERR,"Illegal COMB3 parameter");
-
- nparams++;
}
- delete [] words;
+ MPI_Bcast(&nparams, 1, MPI_INT, 0, world);
+ MPI_Bcast(&maxparam, 1, MPI_INT, 0, world);
+
+ if(comm->me != 0) {
+ params = (Param *) memory->srealloc(params,maxparam*sizeof(Param), "pair:params");
+ }
+
+ MPI_Bcast(params, maxparam*sizeof(Param), MPI_BYTE, 0, world);
}
/* ---------------------------------------------------------------------- */
diff --git a/src/MANYBODY/pair_comb3.h b/src/MANYBODY/pair_comb3.h
index 210f843139..567859127b 100644
--- a/src/MANYBODY/pair_comb3.h
+++ b/src/MANYBODY/pair_comb3.h
@@ -37,8 +37,10 @@ class PairComb3 : public Pair {
virtual double combqeq(double *, int &);
double enegtot;
- // general potential parameters
+ static const int NPARAMS_PER_LINE = 74;
+
protected:
+ // general potential parameters
struct Param {
int ielement,jelement,kelement,powermint;
int ielementgp,jelementgp,kelementgp; //element group
diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp
index a1bc6e1eb4..0013a3271d 100644
--- a/src/MANYBODY/pair_eam.cpp
+++ b/src/MANYBODY/pair_eam.cpp
@@ -29,6 +29,8 @@
#include "error.h"
#include "update.h"
#include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
using namespace LAMMPS_NS;
@@ -462,54 +464,59 @@ void PairEAM::read_file(char *filename)
{
Funcfl *file = &funcfl[nfuncfl-1];
- int me = comm->me;
- FILE *fptr;
- char line[MAXLINE];
+ // read potential file
+ if(comm->me == 0) {
+ PotentialFileReader reader(lmp, filename, "EAM");
- if (me == 0) {
- fptr = force->open_potential(filename);
- if (fptr == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open EAM potential file %s",filename);
- error->one(FLERR,str);
+ try {
+ char * line = nullptr;
+
+ reader.skip_line();
+
+ line = reader.next_line(2);
+ ValueTokenizer values(line);
+ values.next_int(); // ignore
+ file->mass = values.next_double();
+
+ line = reader.next_line(5);
+ values = ValueTokenizer(line);
+ file->nrho = values.next_int();
+ file->drho = values.next_double();
+ file->nr = values.next_int();
+ file->dr = values.next_double();
+ file->cut = values.next_double();
+
+ if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
+ error->one(FLERR,"Invalid EAM potential file");
+
+ memory->create(file->frho, (file->nrho+1), "pair:frho");
+ memory->create(file->rhor, (file->nr+1), "pair:rhor");
+ memory->create(file->zr, (file->nr+1), "pair:zr");
+
+ reader.next_dvector(file->nrho, &file->frho[1]);
+ reader.next_dvector(file->nr, &file->zr[1]);
+ reader.next_dvector(file->nr, &file->rhor[1]);
+ } catch (TokenizerException & e) {
+ error->one(FLERR, e.what());
}
}
- int tmp,nwords;
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- sscanf(line,"%d %lg",&tmp,&file->mass);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- nwords = sscanf(line,"%d %lg %d %lg %lg",
- &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);
+ MPI_Bcast(&file->mass, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->nr, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world);
+
+ if(comm->me != 0) {
+ memory->create(file->frho, (file->nrho+1), "pair:frho");
+ memory->create(file->rhor, (file->nr+1), "pair:rhor");
+ memory->create(file->zr, (file->nr+1), "pair:zr");
}
- MPI_Bcast(&nwords,1,MPI_INT,0,world);
- MPI_Bcast(&file->mass,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->nrho,1,MPI_INT,0,world);
- MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->nr,1,MPI_INT,0,world);
- MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world);
-
- if ((nwords != 5) || (file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
- error->all(FLERR,"Invalid EAM potential file");
-
- memory->create(file->frho,(file->nrho+1),"pair:frho");
- memory->create(file->rhor,(file->nr+1),"pair:rhor");
- memory->create(file->zr,(file->nr+1),"pair:zr");
-
- if (me == 0) grab(fptr,file->nrho,&file->frho[1]);
- MPI_Bcast(&file->frho[1],file->nrho,MPI_DOUBLE,0,world);
-
- if (me == 0) grab(fptr,file->nr,&file->zr[1]);
- MPI_Bcast(&file->zr[1],file->nr,MPI_DOUBLE,0,world);
-
- if (me == 0) grab(fptr,file->nr,&file->rhor[1]);
- MPI_Bcast(&file->rhor[1],file->nr,MPI_DOUBLE,0,world);
-
- if (me == 0) fclose(fptr);
+ MPI_Bcast(&file->frho[1], file->nrho, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->zr[1], file->nr, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->rhor[1], file->nr, MPI_DOUBLE, 0, world);
}
/* ----------------------------------------------------------------------
@@ -782,26 +789,6 @@ void PairEAM::interpolate(int n, double delta, double *f, double **spline)
}
}
-/* ----------------------------------------------------------------------
- grab n values from file fp and put them in list
- values can be several to a line
- only called by proc 0
-------------------------------------------------------------------------- */
-
-void PairEAM::grab(FILE *fptr, int n, double *list)
-{
- char *ptr;
- char line[MAXLINE];
-
- int i = 0;
- while (i < n) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,NULL,error);
- ptr = strtok(line," \t\n\r\f");
- if (ptr) list[i++] = atof(ptr);
- while ((ptr = strtok(NULL," \t\n\r\f"))) list[i++] = atof(ptr);
- }
-}
-
/* ---------------------------------------------------------------------- */
double PairEAM::single(int i, int j, int itype, int jtype,
diff --git a/src/MANYBODY/pair_eam.h b/src/MANYBODY/pair_eam.h
index d98dba1e66..ed525a471c 100644
--- a/src/MANYBODY/pair_eam.h
+++ b/src/MANYBODY/pair_eam.h
@@ -107,7 +107,6 @@ class PairEAM : public Pair {
virtual void allocate();
virtual void array2spline();
void interpolate(int, double, double *, double **);
- void grab(FILE *, int, double *);
virtual void read_file(char *);
virtual void file2array();
diff --git a/src/MANYBODY/pair_eam_alloy.cpp b/src/MANYBODY/pair_eam_alloy.cpp
index a9622f9e07..65ed9b7d2f 100644
--- a/src/MANYBODY/pair_eam_alloy.cpp
+++ b/src/MANYBODY/pair_eam_alloy.cpp
@@ -24,11 +24,11 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
using namespace LAMMPS_NS;
-#define MAXLINE 1024
-
/* ---------------------------------------------------------------------- */
PairEAMAlloy::PairEAMAlloy(LAMMPS *lmp) : PairEAM(lmp)
@@ -61,7 +61,7 @@ void PairEAMAlloy::coeff(int narg, char **arg)
if (setfl) {
for (i = 0; i < setfl->nelements; i++) delete [] setfl->elements[i];
delete [] setfl->elements;
- delete [] setfl->mass;
+ memory->destroy(setfl->mass);
memory->destroy(setfl->frho);
memory->destroy(setfl->rhor);
memory->destroy(setfl->z2r);
@@ -117,98 +117,112 @@ void PairEAMAlloy::read_file(char *filename)
{
Setfl *file = setfl;
- // open potential file
+ // read potential file
+ if(comm->me == 0) {
+ PotentialFileReader reader(lmp, filename, "EAM");
- int me = comm->me;
- FILE *fptr;
- char line[MAXLINE];
+ try {
+ char * line = nullptr;
- if (me == 0) {
- fptr = force->open_potential(filename);
- if (fptr == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open EAM potential file %s",filename);
- error->one(FLERR,str);
+ reader.skip_line();
+ reader.skip_line();
+ reader.skip_line();
+
+ // extract element names from nelements line
+ line = reader.next_line(1);
+ ValueTokenizer values(line);
+ file->nelements = values.next_int();
+
+ if (values.count() != file->nelements + 1)
+ error->one(FLERR,"Incorrect element names in EAM potential file");
+
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) {
+ const std::string word = values.next_string();
+ const int n = word.length() + 1;
+ file->elements[i] = new char[n];
+ strcpy(file->elements[i], word.c_str());
+ }
+
+ //
+
+ line = reader.next_line(5);
+ values = ValueTokenizer(line);
+ file->nrho = values.next_int();
+ file->drho = values.next_double();
+ file->nr = values.next_int();
+ file->dr = values.next_double();
+ file->cut = values.next_double();
+
+ if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
+ error->one(FLERR,"Invalid EAM potential file");
+
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
+
+ for (int i = 0; i < file->nelements; i++) {
+ line = reader.next_line(2);
+ values = ValueTokenizer(line);
+ values.next_int(); // ignore
+ file->mass[i] = values.next_double();
+
+ reader.next_dvector(file->nrho, &file->frho[i][1]);
+ reader.next_dvector(file->nr, &file->rhor[i][1]);
+ }
+
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ reader.next_dvector(file->nr, &file->z2r[i][j][1]);
+ }
+ }
+ } catch (TokenizerException & e) {
+ error->one(FLERR, e.what());
}
}
- // read and broadcast header
- // extract element names from nelements line
+ // broadcast potential information
+ MPI_Bcast(&file->nelements, 1, MPI_INT, 0, world);
- int n;
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- n = strlen(line) + 1;
+ MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->nr, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world);
+
+ // allocate memory on other procs
+ if (comm->me != 0) {
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) file->elements[i] = nullptr;
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
}
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
- sscanf(line,"%d",&file->nelements);
- int nwords = atom->count_words(line);
- if (nwords != file->nelements + 1)
- error->all(FLERR,"Incorrect element names in EAM potential file");
-
- char **words = new char*[file->nelements+1];
- nwords = 0;
- strtok(line," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
-
- file->elements = new char*[file->nelements];
+ // broadcast file->elements string array
for (int i = 0; i < file->nelements; i++) {
- n = strlen(words[i]) + 1;
- file->elements[i] = new char[n];
- strcpy(file->elements[i],words[i]);
- }
- delete [] words;
-
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- nwords = sscanf(line,"%d %lg %d %lg %lg",
- &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);
+ int n;
+ if (comm->me == 0) n = strlen(file->elements[i]) + 1;
+ MPI_Bcast(&n, 1, MPI_INT, 0, world);
+ if (comm->me != 0) file->elements[i] = new char[n];
+ MPI_Bcast(file->elements[i], n, MPI_CHAR, 0, world);
}
- MPI_Bcast(&nwords,1,MPI_INT,0,world);
- MPI_Bcast(&file->nrho,1,MPI_INT,0,world);
- MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->nr,1,MPI_INT,0,world);
- MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world);
+ // broadcast file->mass, frho, rhor
+ for (int i = 0; i < file->nelements; i++) {
+ MPI_Bcast(&file->mass[i], 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->frho[i][1], file->nrho, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->rhor[i][1], file->nr, MPI_DOUBLE, 0, world);
+ }
- if ((nwords != 5) || (file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
- error->all(FLERR,"Invalid EAM potential file");
-
- file->mass = new double[file->nelements];
- memory->create(file->frho,file->nelements,file->nrho+1,"pair:frho");
- memory->create(file->rhor,file->nelements,file->nr+1,"pair:rhor");
- memory->create(file->z2r,file->nelements,file->nelements,file->nr+1,
- "pair:z2r");
-
- int i,j,tmp;
- for (i = 0; i < file->nelements; i++) {
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- sscanf(line,"%d %lg",&tmp,&file->mass[i]);
+ // broadcast file->z2r
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world);
}
- MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world);
-
- if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]);
- MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world);
- if (me == 0) grab(fptr,file->nr,&file->rhor[i][1]);
- MPI_Bcast(&file->rhor[i][1],file->nr,MPI_DOUBLE,0,world);
}
-
- for (i = 0; i < file->nelements; i++)
- for (j = 0; j <= i; j++) {
- if (me == 0) grab(fptr,file->nr,&file->z2r[i][j][1]);
- MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
- }
-
- // close the potential file
-
- if (me == 0) fclose(fptr);
}
/* ----------------------------------------------------------------------
diff --git a/src/MANYBODY/pair_eam_cd.cpp b/src/MANYBODY/pair_eam_cd.cpp
index 28004eae7f..099ffd77cd 100644
--- a/src/MANYBODY/pair_eam_cd.cpp
+++ b/src/MANYBODY/pair_eam_cd.cpp
@@ -28,6 +28,7 @@
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
+#include "tokenizer.h"
using namespace LAMMPS_NS;
@@ -512,16 +513,20 @@ void PairEAMCD::read_h_coeff(char *filename)
while(fgets(nextline, MAXLINE, fptr) != NULL) {
strcpy(line, nextline);
}
- char* ptr = strtok(line, " \t\n\r\f");
- int degree = atoi(ptr);
+
+ ValueTokenizer values(line, " \t\n\r\f");
+ int degree = values.next_int();
nhcoeff = degree+1;
+
+ if (values.count() != nhcoeff + 1 || nhcoeff < 1)
+ error->one(FLERR, "Failed to read h(x) function coefficients in EAM file.");
+
hcoeff = new double[nhcoeff];
+
int i = 0;
- while((ptr = strtok(NULL," \t\n\r\f")) != NULL && i < nhcoeff) {
- hcoeff[i++] = atof(ptr);
+ while(values.has_next()) {
+ hcoeff[i++] = values.next_double();
}
- if (i != nhcoeff || nhcoeff < 1)
- error->one(FLERR,"Failed to read h(x) function coefficients from EAM file.");
// Close the potential file.
diff --git a/src/MANYBODY/pair_eam_fs.cpp b/src/MANYBODY/pair_eam_fs.cpp
index c91e7b5298..fe4fa95268 100644
--- a/src/MANYBODY/pair_eam_fs.cpp
+++ b/src/MANYBODY/pair_eam_fs.cpp
@@ -24,11 +24,11 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
using namespace LAMMPS_NS;
-#define MAXLINE 1024
-
/* ---------------------------------------------------------------------- */
PairEAMFS::PairEAMFS(LAMMPS *lmp) : PairEAM(lmp)
@@ -61,7 +61,7 @@ void PairEAMFS::coeff(int narg, char **arg)
if (fs) {
for (i = 0; i < fs->nelements; i++) delete [] fs->elements[i];
delete [] fs->elements;
- delete [] fs->mass;
+ memory->destroy(fs->mass);
memory->destroy(fs->frho);
memory->destroy(fs->rhor);
memory->destroy(fs->z2r);
@@ -117,103 +117,118 @@ void PairEAMFS::read_file(char *filename)
{
Fs *file = fs;
- // open potential file
+ // read potential file
+ if(comm->me == 0) {
+ PotentialFileReader reader(lmp, filename, "EAM");
- int me = comm->me;
- FILE *fptr;
- char line[MAXLINE];
+ try {
+ char * line = nullptr;
- if (me == 0) {
- fptr = force->open_potential(filename);
- if (fptr == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open EAM potential file %s",filename);
- error->one(FLERR,str);
+ reader.skip_line();
+ reader.skip_line();
+ reader.skip_line();
+
+ // extract element names from nelements line
+ line = reader.next_line(1);
+ ValueTokenizer values(line);
+ file->nelements = values.next_int();
+
+ if (values.count() != file->nelements + 1)
+ error->one(FLERR,"Incorrect element names in EAM potential file");
+
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) {
+ const std::string word = values.next_string();
+ const int n = word.length() + 1;
+ file->elements[i] = new char[n];
+ strcpy(file->elements[i], word.c_str());
+ }
+
+ //
+
+ line = reader.next_line(5);
+ values = ValueTokenizer(line);
+ file->nrho = values.next_int();
+ file->drho = values.next_double();
+ file->nr = values.next_int();
+ file->dr = values.next_double();
+ file->cut = values.next_double();
+
+ if ((file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
+ error->one(FLERR,"Invalid EAM potential file");
+
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
+
+ for (int i = 0; i < file->nelements; i++) {
+ line = reader.next_line(2);
+ values = ValueTokenizer(line);
+ values.next_int(); // ignore
+ file->mass[i] = values.next_double();
+
+ reader.next_dvector(file->nrho, &file->frho[i][1]);
+
+ for (int j = 0; j < file->nelements; j++) {
+ reader.next_dvector(file->nr, &file->rhor[i][j][1]);
+ }
+ }
+
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ reader.next_dvector(file->nr, &file->z2r[i][j][1]);
+ }
+ }
+ } catch (TokenizerException & e) {
+ error->one(FLERR, e.what());
}
}
- // read and broadcast header
- // extract element names from nelements line
+ // broadcast potential information
+ MPI_Bcast(&file->nelements, 1, MPI_INT, 0, world);
- int n;
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- n = strlen(line) + 1;
+ MPI_Bcast(&file->nrho, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->drho, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->nr, 1, MPI_INT, 0, world);
+ MPI_Bcast(&file->dr, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->cut, 1, MPI_DOUBLE, 0, world);
+
+ // allocate memory on other procs
+ if (comm->me != 0) {
+ file->elements = new char*[file->nelements];
+ for (int i = 0; i < file->nelements; i++) file->elements[i] = nullptr;
+ memory->create(file->mass, file->nelements, "pair:mass");
+ memory->create(file->frho, file->nelements, file->nrho + 1, "pair:frho");
+ memory->create(file->rhor, file->nelements, file->nelements, file->nr + 1, "pair:rhor");
+ memory->create(file->z2r, file->nelements, file->nelements, file->nr + 1, "pair:z2r");
}
- MPI_Bcast(&n,1,MPI_INT,0,world);
- MPI_Bcast(line,n,MPI_CHAR,0,world);
- sscanf(line,"%d",&file->nelements);
- int nwords = atom->count_words(line);
- if (nwords != file->nelements + 1)
- error->all(FLERR,"Incorrect element names in EAM potential file");
-
- char **words = new char*[file->nelements+1];
- nwords = 0;
- strtok(line," \t\n\r\f");
- while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
-
- file->elements = new char*[file->nelements];
+ // broadcast file->elements string array
for (int i = 0; i < file->nelements; i++) {
- n = strlen(words[i]) + 1;
- file->elements[i] = new char[n];
- strcpy(file->elements[i],words[i]);
- }
- delete [] words;
-
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- nwords = sscanf(line,"%d %lg %d %lg %lg",
- &file->nrho,&file->drho,&file->nr,&file->dr,&file->cut);
+ int n;
+ if (comm->me == 0) n = strlen(file->elements[i]) + 1;
+ MPI_Bcast(&n, 1, MPI_INT, 0, world);
+ if (comm->me != 0) file->elements[i] = new char[n];
+ MPI_Bcast(file->elements[i], n, MPI_CHAR, 0, world);
}
- MPI_Bcast(&nwords,1,MPI_INT,0,world);
- MPI_Bcast(&file->nrho,1,MPI_INT,0,world);
- MPI_Bcast(&file->drho,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->nr,1,MPI_INT,0,world);
- MPI_Bcast(&file->dr,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&file->cut,1,MPI_DOUBLE,0,world);
+ // broadcast file->mass, frho, rhor
+ for (int i = 0; i < file->nelements; i++) {
+ MPI_Bcast(&file->mass[i], 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&file->frho[i][1], file->nrho, MPI_DOUBLE, 0, world);
- if ((nwords != 5) || (file->nrho <= 0) || (file->nr <= 0) || (file->dr <= 0.0))
- error->all(FLERR,"Invalid EAM potential file");
-
- file->mass = new double[file->nelements];
- memory->create(file->frho,file->nelements,file->nrho+1,
- "pair:frho");
- memory->create(file->rhor,file->nelements,file->nelements,
- file->nr+1,"pair:rhor");
- memory->create(file->z2r,file->nelements,file->nelements,
- file->nr+1,"pair:z2r");
-
- int i,j,tmp;
- for (i = 0; i < file->nelements; i++) {
- if (me == 0) {
- utils::sfgets(FLERR,line,MAXLINE,fptr,filename,error);
- sscanf(line,"%d %lg",&tmp,&file->mass[i]);
- }
- MPI_Bcast(&file->mass[i],1,MPI_DOUBLE,0,world);
-
- if (me == 0) grab(fptr,file->nrho,&file->frho[i][1]);
- MPI_Bcast(&file->frho[i][1],file->nrho,MPI_DOUBLE,0,world);
-
- for (j = 0; j < file->nelements; j++) {
- if (me == 0) grab(fptr,file->nr,&file->rhor[i][j][1]);
- MPI_Bcast(&file->rhor[i][j][1],file->nr,MPI_DOUBLE,0,world);
+ for (int j = 0; j < file->nelements; j++) {
+ MPI_Bcast(&file->rhor[i][j][1], file->nr, MPI_DOUBLE, 0, world);
}
}
- for (i = 0; i < file->nelements; i++)
- for (j = 0; j <= i; j++) {
- if (me == 0) grab(fptr,file->nr,&file->z2r[i][j][1]);
- MPI_Bcast(&file->z2r[i][j][1],file->nr,MPI_DOUBLE,0,world);
+ // broadcast file->z2r
+ for (int i = 0; i < file->nelements; i++) {
+ for (int j = 0; j <= i; j++) {
+ MPI_Bcast(&file->z2r[i][j][1], file->nr, MPI_DOUBLE, 0, world);
}
-
- // close the potential file
-
- if (me == 0) fclose(fptr);
+ }
}
/* ----------------------------------------------------------------------
diff --git a/src/MANYBODY/pair_eim.cpp b/src/MANYBODY/pair_eim.cpp
index dc1c7fa019..6b668eeade 100644
--- a/src/MANYBODY/pair_eim.cpp
+++ b/src/MANYBODY/pair_eim.cpp
@@ -27,11 +27,11 @@
#include "memory.h"
#include "error.h"
#include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
using namespace LAMMPS_NS;
-#define MAXLINE 1024
-
/* ---------------------------------------------------------------------- */
PairEIM::PairEIM(LAMMPS *lmp) : Pair(lmp)
@@ -451,20 +451,6 @@ double PairEIM::init_one(int i, int j)
void PairEIM::read_file(char *filename)
{
- // open potential file
-
- int me = comm->me;
- FILE *fptr;
-
- if (me == 0) {
- fptr = force->open_potential(filename);
- if (fptr == NULL) {
- char str[128];
- snprintf(str,128,"Cannot open EIM potential file %s",filename);
- error->one(FLERR,str);
- }
- }
-
int npair = nelements*(nelements+1)/2;
setfl->ielement = new int[nelements];
setfl->mass = new double[nelements];
@@ -488,52 +474,55 @@ void PairEIM::read_file(char *filename)
setfl->rs = new double[npair];
setfl->tp = new int[npair];
- if (me == 0)
- if (!grabglobal(fptr))
- error->one(FLERR,"Could not grab global entry from EIM potential file");
- MPI_Bcast(&setfl->division,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->rbig,1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->rsmall,1,MPI_DOUBLE,0,world);
+ // read potential file
+ if( comm->me == 0) {
+ EIMPotentialFileReader reader(lmp, filename);
- for (int i = 0; i < nelements; i++) {
- if (me == 0)
- if (!grabsingle(fptr,i))
- error->one(FLERR,"Could not grab element entry from EIM potential file");
- MPI_Bcast(&setfl->ielement[i],1,MPI_INT,0,world);
- MPI_Bcast(&setfl->mass[i],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->negativity[i],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->ra[i],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->ri[i],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->Ec[i],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->q0[i],1,MPI_DOUBLE,0,world);
- }
+ reader.get_global(setfl);
- for (int i = 0; i < nelements; i++) {
- for (int j = i; j < nelements; j++) {
- int ij;
- if (i == j) ij = i;
- else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
- else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
- if (me == 0)
- if (grabpair(fptr,i,j) == 0)
- error->one(FLERR,"Could not grab pair entry from EIM potential file");
- MPI_Bcast(&setfl->rcutphiA[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->rcutphiR[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->Eb[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->r0[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->alpha[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->beta[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->rcutq[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->Asigma[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->rq[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->rcutsigma[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->Ac[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->zeta[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->rs[ij],1,MPI_DOUBLE,0,world);
- MPI_Bcast(&setfl->tp[ij],1,MPI_INT,0,world);
+ for (int i = 0; i < nelements; i++) {
+ reader.get_element(setfl, i, elements[i]);
+ }
+
+ for (int i = 0; i < nelements; i++) {
+ for (int j = i; j < nelements; j++) {
+ int ij;
+ if (i == j) ij = i;
+ else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
+ else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
+ reader.get_pair(setfl, ij, elements[i], elements[j]);
+ }
}
}
+ // broadcast potential information to other procs
+ MPI_Bcast(&setfl->division, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&setfl->rbig, 1, MPI_DOUBLE, 0, world);
+ MPI_Bcast(&setfl->rsmall, 1, MPI_DOUBLE, 0, world);
+
+ MPI_Bcast(setfl->ielement, nelements, MPI_INT, 0, world);
+ MPI_Bcast(setfl->mass, nelements, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->negativity, nelements, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->ra, nelements, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->ri, nelements, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->Ec, nelements, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->q0, nelements, MPI_DOUBLE, 0, world);
+
+ MPI_Bcast(setfl->rcutphiA, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->rcutphiR, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->Eb, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->r0, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->alpha, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->beta, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->rcutq, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->Asigma, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->rq, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->rcutsigma, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->Ac, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->zeta, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->rs, npair, MPI_DOUBLE, 0, world);
+ MPI_Bcast(setfl->tp, npair, MPI_INT, 0, world);
+
setfl->nr = 5000;
setfl->cut = 0.0;
for (int i = 0; i < npair; i++) {
@@ -602,10 +591,6 @@ void PairEIM::read_file(char *filename)
}
}
}
-
- // close the potential file
-
- if (me == 0) fclose(fptr);
}
/* ----------------------------------------------------------------------
@@ -879,114 +864,6 @@ void PairEIM::interpolate(int n, double delta, double *f,
}
}
-/* ----------------------------------------------------------------------
- grab global line from file and store info in setfl
- return 0 if error
-------------------------------------------------------------------------- */
-
-int PairEIM::grabglobal(FILE *fptr)
-{
- char line[MAXLINE];
-
- char *pch = NULL, *data = NULL;
- while (pch == NULL) {
- if (fgets(line,MAXLINE,fptr) == NULL) break;
- pch = strstr(line,"global");
- if (pch != NULL) {
- data = strtok (line," \t\n\r\f");
- data = strtok (NULL,"?");
- sscanf(data,"%lg %lg %lg",&setfl->division,&setfl->rbig,&setfl->rsmall);
- }
- }
- if (pch == NULL) return 0;
- return 1;
-}
-
-/* ----------------------------------------------------------------------
- grab elemental line from file and store info in setfl
- return 0 if error
-------------------------------------------------------------------------- */
-
-int PairEIM::grabsingle(FILE *fptr, int i)
-{
- char line[MAXLINE];
-
- rewind(fptr);
-
- char *pch1 = NULL, *pch2 = NULL, *data = NULL;
- while (pch1 == NULL || pch2 == NULL) {
- if (fgets(line,MAXLINE,fptr) == NULL) break;
- pch1 = strtok (line," \t\n\r\f");
- pch1 = strstr(pch1,"element:");
- if (pch1 != NULL) {
- pch2 = strtok(NULL, " \t\n\r\f");
- if (pch2 != NULL) {
- data = strtok (NULL, "?");
- if (strcmp(pch2,elements[i]) == 0) {
- sscanf(data,"%d %lg %lg %lg %lg %lg %lg",&setfl->ielement[i],
- &setfl->mass[i],&setfl->negativity[i],&setfl->ra[i],
- &setfl->ri[i],&setfl->Ec[i],&setfl->q0[i]);
- } else pch2 = NULL;
- }
- }
- }
- if (pch1 == NULL || pch2 == NULL) return 0;
- return 1;
-}
-
-/* ----------------------------------------------------------------------
- grab pair line from file and store info in setfl
- return 0 if error
-------------------------------------------------------------------------- */
-
-int PairEIM::grabpair(FILE *fptr, int i, int j)
-{
- char line[MAXLINE];
-
- rewind(fptr);
-
- int ij;
- if (i == j) ij = i;
- else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
- else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
-
- char *pch1 = NULL, *pch2 = NULL, *pch3 = NULL, *data = NULL;
- while (pch1 == NULL || pch2 == NULL || pch3 == NULL) {
- if (fgets(line,MAXLINE,fptr) == NULL) break;
- pch1 = strtok (line," \t\n\r\f");
- pch1 = strstr(pch1,"pair:");
- if (pch1 != NULL) {
- pch2 = strtok (NULL, " \t\n\r\f");
- if (pch2 != NULL) pch3 = strtok (NULL, " \t\n\r\f");
- if (pch3 != NULL) data = strtok (NULL, "?");
- if ((pch2 != NULL) && (pch3 != NULL)) {
- if ((strcmp(pch2,elements[i]) == 0 &&
- strcmp(pch3,elements[j]) == 0) ||
- (strcmp(pch2,elements[j]) == 0 &&
- strcmp(pch3,elements[i]) == 0)) {
- sscanf(data,"%lg %lg %lg %lg %lg",
- &setfl->rcutphiA[ij],&setfl->rcutphiR[ij],
- &setfl->Eb[ij],&setfl->r0[ij],&setfl->alpha[ij]);
- utils::sfgets(FLERR,line,MAXLINE,fptr,NULL,error);
- sscanf(line,"%lg %lg %lg %lg %lg",
- &setfl->beta[ij],&setfl->rcutq[ij],&setfl->Asigma[ij],
- &setfl->rq[ij],&setfl->rcutsigma[ij]);
- utils::sfgets(FLERR,line,MAXLINE,fptr,NULL,error);
- sscanf(line,"%lg %lg %lg %d",
- &setfl->Ac[ij],&setfl->zeta[ij],&setfl->rs[ij],
- &setfl->tp[ij]);
- } else {
- pch1 = NULL;
- pch2 = NULL;
- pch3 = NULL;
- }
- }
- }
- }
- if (pch1 == NULL || pch2 == NULL || pch3 == NULL) return 0;
- return 1;
-}
-
/* ----------------------------------------------------------------------
cutoff function
------------------------------------------------------------------------- */
@@ -1171,3 +1048,221 @@ double PairEIM::memory_usage()
bytes += 2 * nmax * sizeof(double);
return bytes;
}
+
+EIMPotentialFileReader::EIMPotentialFileReader(LAMMPS * lmp, const std::string & filename) :
+ Pointers(lmp), filename(filename)
+{
+ if (comm->me != 0) {
+ error->one(FLERR, "EIMPotentialFileReader should only be called by proc 0!");
+ }
+
+ FILE * fp = force->open_potential(filename.c_str());
+
+ if (fp == NULL) {
+ char str[128];
+ snprintf(str, 128, "cannot open EIM potential file %s", filename.c_str());
+ error->one(FLERR, str);
+ }
+
+ parse(fp);
+
+ fclose(fp);
+}
+
+std::pair EIMPotentialFileReader::get_pair(const std::string & a, const std::string & b) {
+ if (a < b) {
+ return std::make_pair(a, b);
+ }
+ return std::make_pair(b, a);
+}
+
+char * EIMPotentialFileReader::next_line(FILE * fp) {
+ // concatenate lines if they end with '&'
+ // strip comments after '#'
+ int n = 0;
+ int nwords = 0;
+ bool concat = false;
+
+ char *ptr = fgets(line, MAXLINE, fp);
+
+ if (ptr == nullptr) {
+ // EOF
+ return nullptr;
+ }
+
+ // strip comment
+ if ((ptr = strchr(line, '#'))) *ptr = '\0';
+
+ // strip ampersand
+ if ((ptr = strrchr(line, '&'))) {
+ concat = true;
+ *ptr = '\0';
+ }
+
+ nwords = utils::count_words(line);
+
+ if (nwords > 0) {
+ n = strlen(line);
+ }
+
+ while(n == 0 || concat) {
+ char *ptr = fgets(&line[n], MAXLINE - n, fp);
+
+ if (ptr == nullptr) {
+ // EOF
+ return line;
+ }
+
+ // strip comment
+ if ((ptr = strchr(line, '#'))) *ptr = '\0';
+
+ // strip ampersand
+ if ((ptr = strrchr(line, '&'))) {
+ concat = true;
+ *ptr = '\0';
+ } else {
+ concat = false;
+ }
+
+ nwords = utils::count_words(line);
+
+ // skip line if blank
+ if (nwords > 0) {
+ n = strlen(line);
+ }
+ }
+
+ return line;
+}
+
+void EIMPotentialFileReader::parse(FILE * fp)
+{
+ char * line = nullptr;
+ bool found_global = false;
+
+ while((line = next_line(fp))) {
+ ValueTokenizer values(line, " \t\r\n\f");
+ std::string type = values.next_string();
+
+ if (type == "global:") {
+ if (values.count() != 4) {
+ error->one(FLERR, "Invalid global line in EIM potential file");
+ }
+
+ division = values.next_double();
+ rbig = values.next_double();
+ rsmall = values.next_double();
+
+ found_global = true;
+ } else if (type == "element:") {
+ if (values.count() != 9) {
+ error->one(FLERR, "Invalid element line in EIM potential file");
+ }
+
+ std::string name = values.next_string();
+
+ ElementData data;
+ data.ielement = values.next_int();
+ data.mass = values.next_double();
+ data.negativity = values.next_double();
+ data.ra = values.next_double();
+ data.ri = values.next_double();
+ data.Ec = values.next_double();
+ data.q0 = values.next_double();
+
+ if (elements.find(name) == elements.end()) {
+ elements[name] = data;
+ } else {
+ error->one(FLERR, "Duplicate pair line in EIM potential file");
+ }
+
+ } else if (type == "pair:") {
+ if (values.count() != 17) {
+ error->one(FLERR, "Invalid element line in EIM potential file");
+ }
+
+ std::string elementA = values.next_string();
+ std::string elementB = values.next_string();
+
+ PairData data;
+ data.rcutphiA = values.next_double();
+ data.rcutphiR = values.next_double();
+ data.Eb = values.next_double();
+ data.r0 = values.next_double();
+ data.alpha = values.next_double();
+ data.beta = values.next_double();
+ data.rcutq = values.next_double();
+ data.Asigma = values.next_double();
+ data.rq = values.next_double();
+ data.rcutsigma = values.next_double();
+ data.Ac = values.next_double();
+ data.zeta = values.next_double();
+ data.rs = values.next_double();
+
+ // should be next_int, but since existing potential files have 1.0000e+00 format
+ // we're doing this instead to keep compatibility
+ data.tp = (int)values.next_double();
+
+ auto p = get_pair(elementA, elementB);
+
+ if (pairs.find(p) == pairs.end()) {
+ pairs[p] = data;
+ } else {
+ error->one(FLERR, "Duplicate pair line in EIM potential file");
+ }
+ }
+ }
+
+ if (!found_global) {
+ error->one(FLERR, "Missing global line in EIM potential file");
+ }
+}
+
+void EIMPotentialFileReader::get_global(PairEIM::Setfl * setfl) {
+ setfl->division = division;
+ setfl->rbig = rbig;
+ setfl->rsmall = rsmall;
+}
+
+void EIMPotentialFileReader::get_element(PairEIM::Setfl * setfl, int i, const std::string & name) {
+ if (elements.find(name) == elements.end()) {
+ char str[128];
+ snprintf(str, 128, "Element %s not defined in EIM potential file", name.c_str());
+ error->one(FLERR, str);
+ }
+
+ ElementData & data = elements[name];
+ setfl->ielement[i] = data.ielement;
+ setfl->mass[i] = data.mass;
+ setfl->negativity[i] = data.negativity;
+ setfl->ra[i] = data.ra;
+ setfl->ri[i] = data.ri;
+ setfl->Ec[i] = data.Ec;
+ setfl->q0[i] = data.q0;
+}
+
+void EIMPotentialFileReader::get_pair(PairEIM::Setfl * setfl, int ij, const std::string & elemA, const std::string & elemB) {
+ auto p = get_pair(elemA, elemB);
+
+ if (pairs.find(p) == pairs.end()) {
+ char str[128];
+ snprintf(str, 128, "Pair (%s, %s) not defined in EIM potential file", elemA.c_str(), elemB.c_str());
+ error->one(FLERR, str);
+ }
+
+ PairData & data = pairs[p];
+ setfl->rcutphiA[ij] = data.rcutphiA;
+ setfl->rcutphiR[ij] = data.rcutphiR;
+ setfl->Eb[ij] = data.Eb;
+ setfl->r0[ij] = data.r0;
+ setfl->alpha[ij] = data.alpha;
+ setfl->beta[ij] = data.beta;
+ setfl->rcutq[ij] = data.rcutq;
+ setfl->Asigma[ij] = data.Asigma;
+ setfl->rq[ij] = data.rq;
+ setfl->rcutsigma[ij] = data.rcutsigma;
+ setfl->Ac[ij] = data.Ac;
+ setfl->zeta[ij] = data.zeta;
+ setfl->rs[ij] = data.rs;
+ setfl->tp[ij] = data.tp;
+}
diff --git a/src/MANYBODY/pair_eim.h b/src/MANYBODY/pair_eim.h
index f9fb2d5a77..986db52453 100644
--- a/src/MANYBODY/pair_eim.h
+++ b/src/MANYBODY/pair_eim.h
@@ -21,6 +21,8 @@ PairStyle(eim,PairEIM)
#define LMP_PAIR_EIM_H
#include "pair.h"
+#include
+#include