diff --git a/doc/.gitignore b/doc/.gitignore
index 35b5e99aee31762e1675c5a067f7628b7a191422..bcd8f6db717a6e3871380308172b0733b68951eb 100644
--- a/doc/.gitignore
+++ b/doc/.gitignore
@@ -1 +1,5 @@
/html
+/LAMMPS.epub
+/LAMMPS.mobi
+/Manual.pdf
+/Developer.pdf
diff --git a/doc/Makefile b/doc/Makefile
index dfae774c61190864ed918b512421f6d12b4c481b..3da7239c5456038fc4e1b53a6a1727646a3f5eef 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -22,7 +22,7 @@ endif
SOURCES=$(wildcard src/*.txt)
OBJECTS=$(SOURCES:src/%.txt=$(RSTDIR)/%.rst)
-.PHONY: help clean-all clean html pdf old venv
+.PHONY: help clean-all clean epub html pdf old venv
# ------------------------------------------
@@ -32,6 +32,7 @@ help:
@echo " pdf create Manual.pdf and Developer.pdf in this dir"
@echo " old create old-style HTML doc pages in old dir"
@echo " fetch fetch HTML and PDF files from LAMMPS web site"
+ @echo " epub create ePUB format manual for e-book readers"
@echo " clean remove all intermediate RST files"
@echo " clean-all reset the entire build environment"
@echo " txt2html build txt2html tool"
@@ -63,6 +64,20 @@ html: $(OBJECTS)
@rm -rf html/USER/*/*.[sg]*
@echo "Build finished. The HTML pages are in doc/html."
+epub: $(OBJECTS)
+ @mkdir -p epub
+ @rm -f LAMMPS.epub
+ @cp src/JPG/lammps-logo.png epub/
+ @(\
+ . $(VENV)/bin/activate ;\
+ cp -r src/* $(RSTDIR)/ ;\
+ sphinx-build -j 8 -b epub -c utils/sphinx-config -d $(BUILDDIR)/doctrees $(RSTDIR) epub ;\
+ deactivate ;\
+ )
+ @mv epub/LAMMPS.epub .
+ @rm -rf epub
+ @echo "Build finished. The ePUB manual file is created."
+
pdf: utils/txt2html/txt2html.exe
@(\
cd src; \
diff --git a/doc/README b/doc/README
index ea0edc0356e6ac4b145235dd03b897dddea77bc6..6db4ba3ca76d604a41a6d19097a91cc0d1ac3efd 100644
--- a/doc/README
+++ b/doc/README
@@ -1,13 +1,14 @@
LAMMPS Documentation
Depending on how you obtained LAMMPS, this directory has 2 or 3
-sub-directories and optionally 2 PDF files:
+sub-directories and optionally 2 PDF files and an ePUB file:
src content files for LAMMPS documentation
html HTML version of the LAMMPS manual (see html/Manual.html)
tools tools and settings for building the documentation
Manual.pdf large PDF version of entire manual
Developer.pdf small PDF with info about how LAMMPS is structured
+LAMMPS.epub Manual in ePUB format
If you downloaded LAMMPS as a tarball from the web site, all these
directories and files should be included.
@@ -49,6 +50,7 @@ make pdf # generate 2 PDF files (Manual.pdf,Developer.pdf)
make old # generate old-style HTML pages in old dir via txt2html
make fetch # fetch HTML doc pages and 2 PDF files from web site
# as a tarball and unpack into html dir and 2 PDFs
+make epub # generate LAMMPS.epub in ePUB format using Sphinx
make clean # remove intermediate RST files created by HTML build
make clean-all # remove entire build folder and any cached data
@@ -91,3 +93,23 @@ This will install virtualenv from the Python Package Index.
----------------
Installing prerequisites for PDF build
+
+[TBA]
+
+----------------
+
+Installing prerequisites for epub build
+
+## ePUB
+
+Same as for HTML. This uses the same tools and configuration
+files as the HTML tree.
+
+For converting the generated ePUB file to a mobi format file
+(for e-book readers like Kindle, that cannot read ePUB), you
+also need to have the 'ebook-convert' tool from the "calibre"
+software installed. http://calibre-ebook.com/
+You first create the ePUB file with 'make epub' and then do:
+
+ebook-convert LAMMPS.epub LAMMPS.mobi
+
diff --git a/doc/src/JPG/lammps-logo.png b/doc/src/JPG/lammps-logo.png
new file mode 100644
index 0000000000000000000000000000000000000000..ae5ba2d8ad56ec55dc7a54ee8be0a21d44003b30
Binary files /dev/null and b/doc/src/JPG/lammps-logo.png differ
diff --git a/doc/src/Manual.txt b/doc/src/Manual.txt
index a957470735540488c6f6de1abff3ea5c3d835645..8590e33acec4e03028a9078392b81f992122067c 100644
--- a/doc/src/Manual.txt
+++ b/doc/src/Manual.txt
@@ -1,7 +1,7 @@
<!-- HTML_ONLY -->
<HEAD>
<TITLE>LAMMPS Users Manual</TITLE>
-<META NAME="docnumber" CONTENT="19 Oct 2016 version">
+<META NAME="docnumber" CONTENT="27 Oct 2016 version">
<META NAME="author" CONTENT="http://lammps.sandia.gov - Sandia National Laboratories">
<META NAME="copyright" CONTENT="Copyright (2003) Sandia Corporation. This software and manual is distributed under the GNU General Public License.">
</HEAD>
@@ -21,7 +21,7 @@
<H1></H1>
LAMMPS Documentation :c,h3
-19 Oct 2016 version :c,h4
+27 Oct 2016 version :c,h4
Version info: :h4
diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index f1b4cae8bbd3dfc7fe49bbc6c62d13c44dda0de2..5acaf8260c66e8d3269b6742fac7760b563eec63 100644
Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ
diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt
index f44f61362da733f1130995f11b93704ff59575e6..4f949648f0bd95abf71d953b859e8aaa51c5f305 100644
--- a/doc/src/Section_commands.txt
+++ b/doc/src/Section_commands.txt
@@ -106,7 +106,7 @@ the $. Thus $\{myTemp\} and $x refer to variable names "myTemp" and
"x".
How the variable is converted to a text string depends on what style
-of variable it is; see the "variable"_variable doc page for details.
+of variable it is; see the "variable"_variable.html doc page for details.
It can be a variable that stores multiple text strings, and return one
of them. The returned text string can be multiple "words" (space
separated) which will then be interpreted as multiple arguments in the
diff --git a/doc/src/Section_errors.txt b/doc/src/Section_errors.txt
index c7ea5844b56e394e4774febe1a4ff8a522f16c38..d3ae9d94b75c7446cfdac8f7d5e1ce9e86aa8595 100644
--- a/doc/src/Section_errors.txt
+++ b/doc/src/Section_errors.txt
@@ -8116,11 +8116,11 @@ boundary of a processor's sub-domain has moved more than 1/2 the
rebuilt and atoms being migrated to new processors. This also means
you may be missing pairwise interactions that need to be computed.
The solution is to change the re-neighboring criteria via the
-"neigh_modify"_neigh_modify command. The safest settings are "delay 0
-every 1 check yes". Second, it may mean that an atom has moved far
-outside a processor's sub-domain or even the entire simulation box.
-This indicates bad physics, e.g. due to highly overlapping atoms, too
-large a timestep, etc. :dd
+"neigh_modify"_neigh_modify.html command. The safest settings are
+"delay 0 every 1 check yes". Second, it may mean that an atom has
+moved far outside a processor's sub-domain or even the entire
+simulation box. This indicates bad physics, e.g. due to highly
+overlapping atoms, too large a timestep, etc. :dd
{Out of range atoms - cannot compute PPPM} :dt
@@ -8132,11 +8132,11 @@ boundary of a processor's sub-domain has moved more than 1/2 the
rebuilt and atoms being migrated to new processors. This also means
you may be missing pairwise interactions that need to be computed.
The solution is to change the re-neighboring criteria via the
-"neigh_modify"_neigh_modify command. The safest settings are "delay 0
-every 1 check yes". Second, it may mean that an atom has moved far
-outside a processor's sub-domain or even the entire simulation box.
-This indicates bad physics, e.g. due to highly overlapping atoms, too
-large a timestep, etc. :dd
+"neigh_modify"_neigh_modify.html command. The safest settings are
+"delay 0 every 1 check yes". Second, it may mean that an atom has
+moved far outside a processor's sub-domain or even the entire
+simulation box. This indicates bad physics, e.g. due to highly
+overlapping atoms, too large a timestep, etc. :dd
{Out of range atoms - cannot compute PPPMDisp} :dt
@@ -8148,11 +8148,11 @@ boundary of a processor's sub-domain has moved more than 1/2 the
rebuilt and atoms being migrated to new processors. This also means
you may be missing pairwise interactions that need to be computed.
The solution is to change the re-neighboring criteria via the
-"neigh_modify"_neigh_modify command. The safest settings are "delay 0
-every 1 check yes". Second, it may mean that an atom has moved far
-outside a processor's sub-domain or even the entire simulation box.
-This indicates bad physics, e.g. due to highly overlapping atoms, too
-large a timestep, etc. :dd
+"neigh_modify"_neigh_modify.html command. The safest settings are
+"delay 0 every 1 check yes". Second, it may mean that an atom has
+moved far outside a processor's sub-domain or even the entire
+simulation box. This indicates bad physics, e.g. due to highly
+overlapping atoms, too large a timestep, etc. :dd
{Overflow of allocated fix vector storage} :dt
diff --git a/doc/src/Section_howto.txt b/doc/src/Section_howto.txt
index 7b8ed2d9a640b433fe8b241ddf304f1212cbd5b7..33cbfa95865b8a2222abf94f7a7f4529459cd08e 100644
--- a/doc/src/Section_howto.txt
+++ b/doc/src/Section_howto.txt
@@ -1864,8 +1864,8 @@ void lammps_close(void *)
int lammps_version(void *)
void lammps_file(void *, char *)
char *lammps_command(void *, char *)
-char *lammps_commands_list(void *, int, char **)
-char *lammps_commands_concatenated(void *, char *)
+void lammps_commands_list(void *, int, char **)
+void lammps_commands_string(void *, char *)
void lammps_free(void *) :pre
The lammps_open() function is used to initialize LAMMPS, passing in a
@@ -1901,51 +1901,66 @@ changes to the LAMMPS command syntax between versions. The returned
LAMMPS version code is an integer (e.g. 2 Sep 2015 results in
20150902) that grows with every new LAMMPS version.
-The lammps_file(), lammps_command(), lammps_commands_list(),
-lammps_commands_concatenated() functions are used to pass one or more
-commands to LAMMPS to execute as if they were coming from an input
-script.
+The lammps_file(), lammps_command(), lammps_commands_list(), and
+lammps_commands_string() functions are used to pass one or more
+commands to LAMMPS to execute, the same as if they were coming from an
+input script.
+
+Via these functions, the calling code can read or generate a series of
+LAMMPS commands one or multiple at a time and pass it thru the library
+interface to setup a problem and then run it in stages. The caller
+can interleave the command function calls with operations it performs,
+calls to extract information from or set information within LAMMPS, or
+calls to another code's library.
The lammps_file() function passes the filename of an input script.
The lammps_command() function passes a single command as a string.
-The lammps_commands_multiple() function passes multiple commands in a
-char** list. The lammps_commands_concatentaed() function passes
-multiple commands concatenated into one long string, separated by
-newline characters. A single command can be spread across multiple
-lines, if the last printable character of all but the last line is
-"&", the same as if the lines appeared in an input script.
-
-Via these library functions, the calling code can read or generate a
-series of LAMMPS commands one or multiple at a time and pass it thru
-the library interface to setup a problem and then run it, interleaving
-the command function calls with operations performed within the
-faller, calls to extract information from LAMMPS, or calls to another
-code's library.
-
-The lammps_free() is a clean-up function to free memory that the library
-allocated previously.
+The lammps_commands_list() function passes multiple commands in a
+char** list. In both lammps_command() and lammps_commands_list(),
+individual commands may or may not have a trailing newline. The
+lammps_commands_string() function passes multiple commands
+concatenated into one long string, separated by newline characters.
+In both lammps_commands_list() and lammps_commands_string(), a single
+command can be spread across multiple lines, if the last printable
+character of all but the last line is "&", the same as if the lines
+appeared in an input script.
+
+The lammps_free() function is a clean-up function to free memory that
+the library allocated previously via other function calls. See
+comments in src/library.cpp file for which other functions need this
+clean-up.
Library.cpp also contains these functions for extracting information
from LAMMPS and setting value within LAMMPS. Again, see the
-documentation in the src/library.cpp file for details:
+documentation in the src/library.cpp file for details, including
+which quantities can be queried by name:
void *lammps_extract_global(void *, char *)
void *lammps_extract_atom(void *, char *)
void *lammps_extract_compute(void *, char *, int, int)
void *lammps_extract_fix(void *, char *, int, int, int, int)
-void *lammps_extract_variable(void *, char *, char *)
+void *lammps_extract_variable(void *, char *, char *) :pre
+
int lammps_set_variable(void *, char *, char *)
-double lammps_get_thermo(void *, char *)
+double lammps_get_thermo(void *, char *) :pre
+
int lammps_get_natoms(void *)
void lammps_gather_atoms(void *, double *)
void lammps_scatter_atoms(void *, double *) :pre
+void lammps_create_atoms(void *, int, tagint *, int *, double *, double *) :pre
+
+The extract functions return a pointer to various global or per-atom
+quantities stored in LAMMPS or to values calculated by a compute, fix,
+or variable. The pointer returned by the extract_global() function
+can be used as a permanent reference to a value which may change. For
+the other extract functions, the underlying storage may be reallocated
+as LAMMPS runs, so you need to re-call the function to assure a
+current pointer or returned value(s).
-These functions can extract various global or per-atom quantities from
-LAMMPS as well as values calculated by a compute, fix, or variable.
-The lammps_set_variable() function can set an existing string-style variable
-to a new value, so that subsequent LAMMPS commands can access the
-variable. The lammps_get_thermo() function returns the current
-value of a thermo keyword.
+The lammps_set_variable() function can set an existing string-style
+variable to a new string value, so that subsequent LAMMPS commands can
+access the variable. The lammps_get_thermo() function returns the
+current value of a thermo keyword as a double.
The lammps_get_natoms() function returns the total number of atoms in
the system and can be used by the caller to allocate space for the
@@ -1956,15 +1971,22 @@ returns a full list to each calling processor. The scatter function
does the inverse. It distributes the same kinds of values,
passed by the caller, to each atom owned by individual processors.
+The lammps_create_atoms() function takes a list of N atoms as input
+with atom types and coords (required), an optionally atom IDs and
+velocities. It uses the coords of each atom to assign it as a new
+atom to the processor that owns it. Additional properties for the new
+atoms can be assigned via the lammps_scatter_atoms() or
+lammps_extract_atom() functions.
+
The examples/COUPLE and python directories have example C++ and C and
Python codes which show how a driver code can link to LAMMPS as a
library, run LAMMPS on a subset of processors, grab data from LAMMPS,
change it, and put it back into LAMMPS.
-NOTE: You can write your own additional functions as needed to define
+NOTE: You can write code for additional functions as needed to define
how your code talks to LAMMPS and add them to src/library.cpp and
src/library.h, as well as to the "Python
-interface"_Section_python.html. The routines you add can access or
+interface"_Section_python.html. The added functions can access or
change any LAMMPS data you wish.
:line
diff --git a/doc/src/Section_packages.txt b/doc/src/Section_packages.txt
index 27ab64d50bf0544ef74d31b6e69b4642285abc93..341483d7a7231ef39ca2595ce5d2c4aa6d3d4ac7 100644
--- a/doc/src/Section_packages.txt
+++ b/doc/src/Section_packages.txt
@@ -1153,6 +1153,7 @@ Package, Description, Author(s), Doc page, Example, Pic/movie, Library
"USER-MISC"_#USER-MISC, single-file contributions, USER-MISC/README, USER-MISC/README, -, -, -
"USER-MANIFOLD"_#USER-MANIFOLD, motion on 2d surface, Stefan Paquay (Eindhoven U of Technology), "fix manifoldforce"_fix_manifoldforce.html, USER/manifold, "manifold"_manifold, -
"USER-MOLFILE"_#USER-MOLFILE, "VMD"_VMD molfile plug-ins, Axel Kohlmeyer (Temple U), "dump molfile"_dump_molfile.html, -, -, VMD-MOLFILE
+"USER-NC-DUMP"_#USER-NC-DUMP, dump output via NetCDF, Lars Pastewka (Karlsruhe Institute of Technology, KIT), "dump nc, dump nc/mpiio"_dump_nc.html, -, -, lib/netcdf
"USER-OMP"_#USER-OMP, OpenMP threaded styles, Axel Kohlmeyer (Temple U), "Section 5.3.4"_accelerate_omp.html, -, -, -
"USER-PHONON"_#USER-PHONON, phonon dynamical matrix, Ling-Ti Kong (Shanghai Jiao Tong U), "fix phonon"_fix_phonon.html, USER/phonon, -, -
"USER-QMMM"_#USER-QMMM, QM/MM coupling, Axel Kohlmeyer (Temple U), "fix qmmm"_fix_qmmm.html, USER/qmmm, -, lib/qmmm
@@ -1598,6 +1599,29 @@ The person who created this package is Axel Kohlmeyer at Temple U
:line
+USER-NC-DUMP package :link(USER-NC-DUMP),h5
+
+Contents: Dump styles for writing NetCDF format files. NetCDF is a binary,
+portable, self-describing file format on top of HDF5. The file format
+contents follow the AMBER NetCDF trajectory conventions
+(http://ambermd.org/netcdf/nctraj.xhtml), but include extensions to this
+convention. This package implements a "dump nc"_dump_nc.html command
+and a "dump nc/mpiio"_dump_nc.html command to output LAMMPS snapshots
+in this format. See src/USER-NC-DUMP/README for more details.
+
+NetCDF files can be directly visualized with the following tools:
+Ovito (http://www.ovito.org/). Ovito supports the AMBER convention
+ and all of the above extensions. :ulb,l
+VMD (http://www.ks.uiuc.edu/Research/vmd/) :l
+AtomEye (http://www.libatoms.org/). The libAtoms version of AtomEye contains
+ a NetCDF reader that is not present in the standard distribution of AtomEye :l,ule
+
+The person who created these files is Lars Pastewka at
+Karlsruhe Institute of Technology (lars.pastewka at kit.edu).
+Contact him directly if you have questions.
+
+:line
+
USER-OMP package :link(USER-OMP),h5
Supporting info:
diff --git a/doc/src/Section_python.txt b/doc/src/Section_python.txt
index 352c1fa5927b50c3a595295ce612385c174f00af..b3b5171e74ac291429c23768b63f0da2706d14d8 100644
--- a/doc/src/Section_python.txt
+++ b/doc/src/Section_python.txt
@@ -550,7 +550,8 @@ version = lmp.version() # return the numerical version id, e.g. LAMMPS 2 Sep 20
lmp.file(file) # run an entire input script, file = "in.lj"
lmp.command(cmd) # invoke a single LAMMPS command, cmd = "run 100" :pre
-lmp.commands_list(cmdlist) # invoke commands in cmdlist
+lmp.commands_list(cmdlist) # invoke commands in cmdlist = ["run 10", "run 20"]
+lmp.commands_string(multicmd) # invoke commands in multicmd = "run 10\nrun 20"
xlo = lmp.extract_global(name,type) # extract a global quantity
# name = "boxxlo", "nlocal", etc
@@ -631,8 +632,9 @@ lmp2 = lammps()
lmp1.file("in.file1")
lmp2.file("in.file2") :pre
-The file() and command() methods allow an input script or single
-commands to be invoked.
+The file(), command(), commands_list(), commands_string() methods
+allow an input script, a single command, or multiple commands to be
+invoked.
The extract_global(), extract_atom(), extract_compute(),
extract_fix(), and extract_variable() methods return values or
diff --git a/doc/src/Section_start.txt b/doc/src/Section_start.txt
index 4d46ce82ae6132dbe3cc31fcec120a0a5187ee5d..a0cc792455564ce06057789b5444fd4bd44eda79 100644
--- a/doc/src/Section_start.txt
+++ b/doc/src/Section_start.txt
@@ -1601,9 +1601,9 @@ implementations, either by environment variables that specify how to
order physical processors, or by config files that specify what
physical processors to assign to each MPI rank. The -reorder switch
simply gives you a portable way to do this without relying on MPI
-itself. See the "processors out"_processors command for how to output
-info on the final assignment of physical processors to the LAMMPS
-simulation domain.
+itself. See the "processors out"_processors.html command for how
+to output info on the final assignment of physical processors to
+the LAMMPS simulation domain.
-screen file :pre
diff --git a/doc/src/accelerate_intel.txt b/doc/src/accelerate_intel.txt
index cda5a80c4e463b6e7ae1fc44315be0d275ad7477..3e94e7762acb8a516122d38dec690db43ecf0262 100644
--- a/doc/src/accelerate_intel.txt
+++ b/doc/src/accelerate_intel.txt
@@ -151,7 +151,7 @@ can start running so that the CPU pipeline is still being used
efficiently. Although benefits can be seen by launching a MPI task
for every hardware thread, for multinode simulations, we recommend
that OpenMP threads are used for SMT instead, either with the
-USER-INTEL package, "USER-OMP package"_accelerate_omp.html", or
+USER-INTEL package, "USER-OMP package"_accelerate_omp.html, or
"KOKKOS package"_accelerate_kokkos.html. In the example above, up
to 36X speedups can be observed by using all 36 physical cores with
LAMMPS. By using all 72 hardware threads, an additional 10-30%
@@ -343,7 +343,7 @@ when using offload.
Not all styles are supported in the USER-INTEL package. You can mix
the USER-INTEL package with styles from the "OPT"_accelerate_opt.html
-package or the "USER-OMP package"_accelerate_omp.html". Of course,
+package or the "USER-OMP package"_accelerate_omp.html. Of course,
this requires that these packages were installed at build time. This
can performed automatically by using "-sf hybrid intel opt" or
"-sf hybrid intel omp" command-line options. Alternatively, the "opt"
diff --git a/doc/src/atom_style.txt b/doc/src/atom_style.txt
index e632d0018d8843e52f1db461339e5784d421b2b9..9ba817a663b477987c86be89a61be95a23e97984 100644
--- a/doc/src/atom_style.txt
+++ b/doc/src/atom_style.txt
@@ -166,7 +166,7 @@ stores a per-particle mass and size and orientation (i.e. the corner
points of the triangle).
The {template} style allows molecular topolgy (bonds,angles,etc) to be
-defined via a molecule template using the "molecule"_molecule.txt
+defined via a molecule template using the "molecule"_molecule.html
command. The template stores one or more molecules with a single copy
of the topology info (bonds,angles,etc) of each. Individual atoms
only store a template index and template atom to identify which
diff --git a/doc/src/bond_fene.txt b/doc/src/bond_fene.txt
index d7ab09c48f7ee9e7422a9283edef4c831a70c645..a4dd393d8d82afa1ca8194d25c4ff7f8f6dbc455 100644
--- a/doc/src/bond_fene.txt
+++ b/doc/src/bond_fene.txt
@@ -73,7 +73,7 @@ This bond style can only be used if LAMMPS was built with the
MOLECULE package (which it is by default). See the "Making
LAMMPS"_Section_start.html#start_3 section for more info on packages.
-You typically should specify "special_bonds fene"_special_bonds.html"
+You typically should specify "special_bonds fene"_special_bonds.html
or "special_bonds lj/coul 0 1 1"_special_bonds.html to use this bond
style. LAMMPS will issue a warning it that's not the case.
diff --git a/doc/src/bond_fene_expand.txt b/doc/src/bond_fene_expand.txt
index 4342231d20184657b73f8b5e9f11a62b95681840..6ddd6d4876a343cb53525cdaf351a81386e60d54 100644
--- a/doc/src/bond_fene_expand.txt
+++ b/doc/src/bond_fene_expand.txt
@@ -76,7 +76,7 @@ This bond style can only be used if LAMMPS was built with the
MOLECULE package (which it is by default). See the "Making
LAMMPS"_Section_start.html#start_3 section for more info on packages.
-You typically should specify "special_bonds fene"_special_bonds.html"
+You typically should specify "special_bonds fene"_special_bonds.html
or "special_bonds lj/coul 0 1 1"_special_bonds.html to use this bond
style. LAMMPS will issue a warning it that's not the case.
diff --git a/doc/src/compute_fep.txt b/doc/src/compute_fep.txt
index 1db74fef064795f52b9dd57d9ff150853d0ae015..9bbae7b20f1c835c67d7808867d4a398a79b4541 100644
--- a/doc/src/compute_fep.txt
+++ b/doc/src/compute_fep.txt
@@ -236,7 +236,7 @@ LAMMPS"_Section_start.html#start_3 section for more info.
[Related commands:]
"fix adapt/fep"_fix_adapt_fep.html, "fix ave/time"_fix_ave_time.html,
-"pair_lj_soft_coul_soft"_pair_lj_soft_coul_soft.txt
+"pair_style lj/soft/coul/soft"_pair_lj_soft.html
[Default:]
diff --git a/doc/src/compute_orientorder_atom.txt b/doc/src/compute_orientorder_atom.txt
index 58378b009bdaa5693fc41f8059b2428d594c7c0b..c5ecef3cb29ee3656a79f01965ddcf62a0687baa 100644
--- a/doc/src/compute_orientorder_atom.txt
+++ b/doc/src/compute_orientorder_atom.txt
@@ -15,7 +15,7 @@ compute ID group-ID orientorder/atom keyword values ... :pre
ID, group-ID are documented in "compute"_compute.html command :ulb,l
orientorder/atom = style name of this compute command :l
one or more keyword/value pairs may be appended :l
-keyword = {cutoff} or {nnn} or {ql}
+keyword = {cutoff} or {nnn} or {degrees}
{cutoff} value = distance cutoff
{nnn} value = number of nearest neighbors
{degrees} values = nlvalues, l1, l2,... :pre
@@ -111,7 +111,7 @@ options.
[Default:]
-The option defaults are {cutoff} = pair style cutoff, {nnn} = 12, {degrees} = 5 4 6 8 9 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12.
+The option defaults are {cutoff} = pair style cutoff, {nnn} = 12, {degrees} = 5 4 6 8 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12.
:line
diff --git a/doc/src/compute_property_local.txt b/doc/src/compute_property_local.txt
index bc2ecca46671b8e656bbd9a1902974b009358c4f..f7851e864bfbb14cef1a8e055592018a00ba609c 100644
--- a/doc/src/compute_property_local.txt
+++ b/doc/src/compute_property_local.txt
@@ -78,7 +78,7 @@ defined by the "pair_style"_pair_style.html command for the types of
the two atoms is used. For the {radius} setting, the sum of the radii
of the two particles is used as a cutoff. For example, this is
appropriate for granular particles which only interact when they are
-overlapping, as computed by "granular pair styles"_pair_gran.txt.
+overlapping, as computed by "granular pair styles"_pair_gran.html.
If the inputs are bond, angle, etc attributes, the local data is
generated by looping over all the atoms owned on a processor and
diff --git a/doc/src/fix_phonon.txt b/doc/src/fix_phonon.txt
index 254cf4740806844540ac3b5a4e9dff123356b7c6..10a0c95c80191ed8fd5b75097d68ac68aff20e6a 100644
--- a/doc/src/fix_phonon.txt
+++ b/doc/src/fix_phonon.txt
@@ -184,7 +184,7 @@ This fix requires LAMMPS be built with an FFT library. See the
[Default:]
The option defaults are sysdim = the same dimemsion as specified by
-the "dimension"_dimension command, and nasr = 20.
+the "dimension"_dimension.html command, and nasr = 20.
:line
diff --git a/doc/src/fix_property_atom.txt b/doc/src/fix_property_atom.txt
index 10dd08738f24a93e4b6f6f9f9e478e956f8675bb..8c9164c25d91d533739d412ab4f0dad1bdc5fc72 100644
--- a/doc/src/fix_property_atom.txt
+++ b/doc/src/fix_property_atom.txt
@@ -177,7 +177,7 @@ their values. This means that the values can be output via the "dump
custom"_dump.html command, accessed by fixes like "fix
ave/atom"_fix_ave_atom.html, accessed by other computes like "compute
reduce"_compute_reduce.html, or used in "atom-style
-variables"_variables.
+variables"_variable.html.
For example, these commands will output two new properties to a custom
dump file:
diff --git a/doc/src/fix_rigid.txt b/doc/src/fix_rigid.txt
index 685e2694d3a2283d6b550d5ec644e484e5072c45..a38948e560048aca0915b4b0f532f37ba876b63f 100644
--- a/doc/src/fix_rigid.txt
+++ b/doc/src/fix_rigid.txt
@@ -620,7 +620,7 @@ rigid styles for the rigid bodies. :l
Use "fix press/berendsen"_fix_press_berendsen.html to compute the
pressure and change the box dimensions. Use one of the 4 NVE or 2 NVT
-rigid styles for the rigid bodies. Use "fix nvt"_fix_nh.thml (or any
+rigid styles for the rigid bodies. Use "fix nvt"_fix_nh.html (or any
other thermostat) for the non-rigid particles. :l
:ule
diff --git a/doc/src/fix_ti_spring.txt b/doc/src/fix_ti_spring.txt
index 64ffe2fc8ac0cac931612aaa64878fa776b5a921..40e595e21e41d4a901d9d2afd78e97d771413371 100644
--- a/doc/src/fix_ti_spring.txt
+++ b/doc/src/fix_ti_spring.txt
@@ -104,7 +104,7 @@ the Nose-Hoover thermostat ("fix nvt"_fix_nh.html) is {NOT}
recommended due to its well documented issues with the canonical
sampling of harmonic degrees of freedom (notice that the {chain}
option will {NOT} solve this problem). The Langevin thermostat ("fix
-langevin"_fix_langevin.html") correctly thermostats the system and we
+langevin"_fix_langevin.html) correctly thermostats the system and we
advise its usage with ti/spring command.
[Restart, fix_modify, output, run start/stop, minimize info:]
diff --git a/doc/src/neb.txt b/doc/src/neb.txt
index a2f8161ee7c8fad6dabe7799ffdafbd19c7f3ed5..17bd8544b780ad44ad12dfaa1e6833875f7f979d 100644
--- a/doc/src/neb.txt
+++ b/doc/src/neb.txt
@@ -83,9 +83,9 @@ replica. Conceptually, the non-NEB atoms provide a background force
field for the NEB atoms. They can be allowed to move during the NEB
minimiation procedure (which will typically induce different
coordinates for non-NEB atoms in different replicas), or held fixed
-using other LAMMPS commands such as "fix setforce"_fix_setforce. Note
-that the "partition"_partition.html command can be used to invoke a
-command on a subset of the replicas, e.g. if you wish to hold NEB or
+using other LAMMPS commands such as "fix setforce"_fix_setforce.html.
+Note that the "partition"_partition.html command can be used to invoke
+a command on a subset of the replicas, e.g. if you wish to hold NEB or
non-NEB atoms fixed in only the end-point replicas.
The initial atomic configuration for each of the replicas can be
diff --git a/doc/src/pair_hbond_dreiding.txt b/doc/src/pair_hbond_dreiding.txt
index caaa715a7bd542a5c7f7e1eafb28e2bf38fe2b98..b3c9b5b5d50d6fc4cbdb01e45b6a9c332b4e9289 100644
--- a/doc/src/pair_hbond_dreiding.txt
+++ b/doc/src/pair_hbond_dreiding.txt
@@ -138,8 +138,8 @@ angle cutoff (degrees) :ul
A single hydrogen atom type K can be specified, or a wild-card
asterisk can be used in place of or in conjunction with the K
arguments to select multiple types as hydrogens. This takes the form
-"*" or "*n" or "n*" or "m*n". See the "pair_coeff"_pair_coeff command
-doc page for details.
+"*" or "*n" or "n*" or "m*n". See the "pair_coeff"_pair_coeff.html
+command doc page for details.
If the donor flag is {i}, then the atom of type I in the pair_coeff
command is treated as the donor, and J is the acceptor. If the donor
diff --git a/doc/src/pair_line_lj.txt b/doc/src/pair_line_lj.txt
index 6a72e296a6e7671cc39e5b59df115a0f311dfca1..bb1146b64e16c5cf220481a032265b6ac2a6895f 100644
--- a/doc/src/pair_line_lj.txt
+++ b/doc/src/pair_line_lj.txt
@@ -60,8 +60,8 @@ pair_style command or overridden with an optional argument in the
pair_coeff command for a type pair as discussed below. The distance
between the centers of 2 line segments, or the center of a line
segment and a point particle, must be less than this distance (plus
-the neighbor skin; see the "neighbor"_neighbor command), for the pair
-of particles to be included in the neighbor list.
+the neighbor skin; see the "neighbor"_neighbor.html command), for
+the pair of particles to be included in the neighbor list.
NOTE: This means that a too-short value for the {cutoff} setting can
exclude a pair of particles from the neighbor list even if pairs of
diff --git a/doc/src/pair_lubricateU.txt b/doc/src/pair_lubricateU.txt
index fb1ce8d89bd69fe8d334efc31eb81a9375832023..52e63ac35016ffbda357193527092a3ddfc2b8bc 100644
--- a/doc/src/pair_lubricateU.txt
+++ b/doc/src/pair_lubricateU.txt
@@ -119,7 +119,7 @@ of walls (whether moving or stationary) will affect the volume
fraction available to colloidal particles. This is currently accounted
for with the following types of walls: "wall/lj93"_fix_wall.html,
"wall/lj126"_fix_wall.html, "wall/colloid"_fix_wall.html, and
-"wall/harmonic_fix_wall.html". For these wall styles, the correct
+"wall/harmonic"_fix_wall.html. For these wall styles, the correct
volume fraction will be used when walls do not coincide with the box
boundary, as well as when walls move and thereby cause a change in the
volume fraction. To use these wall styles with pair_style {lubricateU}
diff --git a/doc/src/pair_thole.txt b/doc/src/pair_thole.txt
index cf43e6967393cd358e05e5231630f07839f7c090..923485eaa7a3d144cd327ca6267082bb096eb089 100644
--- a/doc/src/pair_thole.txt
+++ b/doc/src/pair_thole.txt
@@ -180,7 +180,7 @@ package.
langevin/drude"_fix_langevin_drude.html, "fix
drude/transform"_fix_drude_transform.html, "compute
temp/drude"_compute_temp_drude.html
-"pair_style lj/cut/coul/long"_pair_lj_cut_coul_long
+"pair_style lj/cut/coul/long"_pair_lj.html
[Default:] none
diff --git a/doc/src/read_data.txt b/doc/src/read_data.txt
index ce617d5d0db09deccb0ed18b8cb7f0463bd6d4f3..9ddc9d89218bb79d9496fce812f8eb306c8a3800 100644
--- a/doc/src/read_data.txt
+++ b/doc/src/read_data.txt
@@ -431,8 +431,8 @@ Atoms # sphere
Pair Coeffs # lj/cut :pre
will check if the currently-defined "atom_style"_atom_style.html is
-{sphere}, and the current "pair_style"_pair_style is {lj/cut}. If
-not, LAMMPS will issue a warning to indicate that the data file
+{sphere}, and the current "pair_style"_pair_style.html is {lj/cut}.
+If not, LAMMPS will issue a warning to indicate that the data file
section likely does not contain the correct number or type of
parameters expected for the currently-defined style.
diff --git a/doc/src/read_dump.txt b/doc/src/read_dump.txt
index 30345fd6cdae1de023d777a8c2577fdee422bdd4..f6424f65291793443f54707bfac985c20fb431bd 100644
--- a/doc/src/read_dump.txt
+++ b/doc/src/read_dump.txt
@@ -322,3 +322,6 @@ They are only enabled if LAMMPS was built with that packages. See the
The option defaults are box = yes, replace = yes, purge = no, trim =
no, add = no, scaled = no, wrapped = yes, and format = native.
+
+:link(vmd,http://www.ks.uiuc.edu/Research/vmd)
+
diff --git a/doc/src/read_restart.txt b/doc/src/read_restart.txt
index c79ca4460df1f25793abe14a2dcf3f3ede7aeea6..d59b5313fdf785d5816400f66f424255d5897dbc 100644
--- a/doc/src/read_restart.txt
+++ b/doc/src/read_restart.txt
@@ -141,11 +141,11 @@ these settings after the restart file is read.
"units"_units.html
"newton bond"_newton.html (see discussion of newton command below)
"atom style"_atom_style.html and "atom_modify"_atom_modify.html settings id, map, sort
-"comm style"_comm_style.html and "comm_modify"_comm_modify settings mode, cutoff, vel
+"comm style"_comm_style.html and "comm_modify"_comm_modify.html settings mode, cutoff, vel
"timestep"_timestep.html
simulation box size and shape and "boundary"_boundary.html settings
atom "group"_group.html definitions
-per-type atom settings such as "mass"_mass.thml
+per-type atom settings such as "mass"_mass.html
per-atom attributes including their group assignments and molecular topology attributes (bonds, angles, etc)
force field styles ("pair"_pair_style.html, "bond"_bond_style.html, "angle"_angle_style.html, etc)
force field coefficients ("pair"_pair_coeff.html, "bond"_bond_coeff.html, "angle"_angle_coeff.html, etc) in some cases (see below)
diff --git a/doc/src/thermo_style.txt b/doc/src/thermo_style.txt
index 5e662e85b90a6d6d9be675e361fbb8e6ef5270a6..6cdb8a3cc87871b88c0583bf8766ba910e12ee21 100644
--- a/doc/src/thermo_style.txt
+++ b/doc/src/thermo_style.txt
@@ -108,7 +108,7 @@ Style {custom} is the most general setting and allows you to specify
which of the keywords listed above you want printed on each
thermodynamic timestep. Note that the keywords c_ID, f_ID, v_name are
references to "computes"_compute.html, "fixes"_fix.html, and
-equal-style "variables"_variable.html" that have been defined
+equal-style "variables"_variable.html that have been defined
elsewhere in the input script or can even be new styles which users
have added to LAMMPS (see the "Section 10"_Section_modify.html
section of the documentation). Thus the {custom} style provides a
diff --git a/doc/src/write_dump.txt b/doc/src/write_dump.txt
index cf8302121c71efe321b3010f635f7db02a517f5b..ae32a94e840d32d0c25bd384bcca70a4fc51e596 100644
--- a/doc/src/write_dump.txt
+++ b/doc/src/write_dump.txt
@@ -26,7 +26,7 @@ write_dump all atom dump.atom
write_dump subgroup atom dump.run.bin
write_dump all custom dump.myforce.* id type x y vx fx
write_dump flow custom dump.%.myforce id type c_myF\[3\] v_ke modify sort id
-write_dump all xyz system.xyz modify sort id elements O H
+write_dump all xyz system.xyz modify sort id element O H
write_dump all image snap*.jpg type type size 960 960 modify backcolor white
write_dump all image snap*.jpg element element &
bond atom 0.3 shiny 0.1 ssao yes 6345 0.2 size 1600 1600 &
diff --git a/doc/utils/sphinx-config/conf.py b/doc/utils/sphinx-config/conf.py
index d89217834623a5791559c07bc95407f9c2f68af8..2b4950519a44d41d10501e5a0683242442cec09b 100644
--- a/doc/utils/sphinx-config/conf.py
+++ b/doc/utils/sphinx-config/conf.py
@@ -276,4 +276,27 @@ texinfo_documents = [
# If true, do not generate a @detailmenu in the "Top" node's menu.
#texinfo_no_detailmenu = False
+# -- Options for ePUB output ----------------------------------------------
+
+epub_title = 'LAMMPS Documentation - ' + get_lammps_version()
+
+epub_cover = ('lammps-logo.png', '')
+
+epub_description = """
+This is the Manual for the LAMMPS software package.
+
+LAMMPS stands for Large-scale Atomic/Molecular Massively Parallel
+Simulator and is a classical molecular dynamics simulation code
+designed to run efficiently on parallel computers. It was developed
+at Sandia National Laboratories, a US Department of Energy facility,
+with funding from the DOE. It is an open-source code, distributed
+freely under the terms of the GNU Public License (GPL).
+
+The primary author of the code is Steve Plimpton, who can be emailed
+at sjplimp@sandia.gov. The LAMMPS WWW Site at lammps.sandia.gov has
+more information about the code and its uses.
+"""
+
+epub_author = 'The LAMMPS Developers'
+
diff --git a/examples/COUPLE/simple/in.lj b/examples/COUPLE/simple/in.lj
index 5f6ebdc4f2d414babe4a03a38f3180ce84e2737e..3dc80bdfea750f6b7e12cf81929d4b0304c11929 100644
--- a/examples/COUPLE/simple/in.lj
+++ b/examples/COUPLE/simple/in.lj
@@ -20,4 +20,6 @@ neigh_modify delay 0 every 20 check no
fix 1 all nve
+variable fx atom fx
+
run 10
diff --git a/examples/COUPLE/simple/simple.c b/examples/COUPLE/simple/simple.c
index c50f289b7af63c2d595cec213b7764f957be9270..cc813d78fe785cd695344e8f939c356013a7cd3a 100644
--- a/examples/COUPLE/simple/simple.c
+++ b/examples/COUPLE/simple/simple.c
@@ -71,8 +71,8 @@ int main(int narg, char **arg)
(could just send it to proc 0 of comm_lammps and let it Bcast)
all LAMMPS procs call lammps_command() on the line */
- void *ptr;
- if (lammps == 1) lammps_open(0,NULL,comm_lammps,&ptr);
+ void *lmp = NULL;
+ if (lammps == 1) lammps_open(0,NULL,comm_lammps,&lmp);
int n;
char line[1024];
@@ -85,7 +85,7 @@ int main(int narg, char **arg)
MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
if (n == 0) break;
MPI_Bcast(line,n,MPI_CHAR,0,MPI_COMM_WORLD);
- if (lammps == 1) lammps_command(ptr,line);
+ if (lammps == 1) lammps_command(lmp,line);
}
/* run 10 more steps
@@ -94,23 +94,72 @@ int main(int narg, char **arg)
put coords back into LAMMPS
run a single step with changed coords */
+ double *x = NULL;
+ double *v = NULL;
+
if (lammps == 1) {
- lammps_command(ptr,"run 10");
+ lammps_command(lmp,"run 10");
- int natoms = lammps_get_natoms(ptr);
- double *x = (double *) malloc(3*natoms*sizeof(double));
- lammps_gather_atoms(ptr,"x",1,3,x);
+ int natoms = lammps_get_natoms(lmp);
+ x = (double *) malloc(3*natoms*sizeof(double));
+ lammps_gather_atoms(lmp,"x",1,3,x);
+ v = (double *) malloc(3*natoms*sizeof(double));
+ lammps_gather_atoms(lmp,"v",1,3,v);
double epsilon = 0.1;
x[0] += epsilon;
- lammps_scatter_atoms(ptr,"x",1,3,x);
- free(x);
+ lammps_scatter_atoms(lmp,"x",1,3,x);
+
+ lammps_command(lmp,"run 1");
+ }
+
+ // extract force on single atom two different ways
+
+ if (lammps == 1) {
+ double **f = (double **) lammps_extract_atom(lmp,"f");
+ printf("Force on 1 atom via extract_atom: %g\n",f[0][0]);
+
+ double *fx = (double *) lammps_extract_variable(lmp,"fx","all");
+ printf("Force on 1 atom via extract_variable: %g\n",fx[0]);
+ }
- lammps_command(ptr,"run 1");
+ /* use commands_string() and commands_list() to invoke more commands */
+
+ char *strtwo = "run 10\nrun 20";
+ if (lammps == 1) lammps_commands_string(lmp,strtwo);
+
+ char *cmds[2];
+ cmds[0] = "run 10";
+ cmds[1] = "run 20";
+ if (lammps == 1) lammps_commands_list(lmp,2,cmds);
+
+ /* delete all atoms
+ create_atoms() to create new ones with old coords, vels
+ initial thermo should be same as step 20 */
+
+ int *type = NULL;
+
+ if (lammps == 1) {
+ int natoms = lammps_get_natoms(lmp);
+ type = (int *) malloc(natoms*sizeof(double));
+ int i;
+ for (i = 0; i < natoms; i++) type[i] = 1;
+
+ lammps_command(lmp,"delete_atoms group all");
+ lammps_create_atoms(lmp,natoms,NULL,type,x,v);
+ lammps_command(lmp,"run 10");
}
- if (lammps == 1) lammps_close(ptr);
+ if (x) free(x);
+ if (v) free(v);
+ if (type) free(type);
+
+ // close down LAMMPS
+
+ lammps_close(lmp);
/* close down MPI */
+ if (lammps == 1) MPI_Comm_free(&comm_lammps);
+ MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
}
diff --git a/examples/COUPLE/simple/simple.cpp b/examples/COUPLE/simple/simple.cpp
index a440bcb1ae9c711dad0cf9dbe9d046071fb62b95..b279ae704cb9febd3401fe420ec0a7db292c5189 100644
--- a/examples/COUPLE/simple/simple.cpp
+++ b/examples/COUPLE/simple/simple.cpp
@@ -77,7 +77,7 @@ int main(int narg, char **arg)
// (could just send it to proc 0 of comm_lammps and let it Bcast)
// all LAMMPS procs call input->one() on the line
- LAMMPS *lmp;
+ LAMMPS *lmp = NULL;
if (lammps == 1) lmp = new LAMMPS(0,NULL,comm_lammps);
int n;
@@ -91,7 +91,7 @@ int main(int narg, char **arg)
MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);
if (n == 0) break;
MPI_Bcast(line,n,MPI_CHAR,0,MPI_COMM_WORLD);
- if (lammps == 1) lmp->input->one(line);
+ if (lammps == 1) lammps_command(lmp,line);
}
// run 10 more steps
@@ -100,23 +100,74 @@ int main(int narg, char **arg)
// put coords back into LAMMPS
// run a single step with changed coords
+ double *x = NULL;
+ double *v = NULL;
+
if (lammps == 1) {
lmp->input->one("run 10");
int natoms = static_cast<int> (lmp->atom->natoms);
- double *x = new double[3*natoms];
+ x = new double[3*natoms];
+ v = new double[3*natoms];
lammps_gather_atoms(lmp,"x",1,3,x);
+ lammps_gather_atoms(lmp,"v",1,3,v);
double epsilon = 0.1;
x[0] += epsilon;
lammps_scatter_atoms(lmp,"x",1,3,x);
- delete [] x;
+ // these 2 lines are the same
+
+ // lammps_command(lmp,"run 1");
lmp->input->one("run 1");
}
- if (lammps == 1) delete lmp;
+ // extract force on single atom two different ways
+
+ if (lammps == 1) {
+ double **f = (double **) lammps_extract_atom(lmp,"f");
+ printf("Force on 1 atom via extract_atom: %g\n",f[0][0]);
+
+ double *fx = (double *) lammps_extract_variable(lmp,"fx","all");
+ printf("Force on 1 atom via extract_variable: %g\n",fx[0]);
+ }
+
+ // use commands_string() and commands_list() to invoke more commands
+
+ char *strtwo = "run 10\nrun 20";
+ if (lammps == 1) lammps_commands_string(lmp,strtwo);
+
+ char *cmds[2];
+ cmds[0] = "run 10";
+ cmds[1] = "run 20";
+ if (lammps == 1) lammps_commands_list(lmp,2,cmds);
+
+ // delete all atoms
+ // create_atoms() to create new ones with old coords, vels
+ // initial thermo should be same as step 20
+
+ int *type = NULL;
+
+ if (lammps == 1) {
+ int natoms = static_cast<int> (lmp->atom->natoms);
+ type = new int[natoms];
+ for (int i = 0; i < natoms; i++) type[i] = 1;
+
+ lmp->input->one("delete_atoms group all");
+ lammps_create_atoms(lmp,natoms,NULL,type,x,v);
+ lmp->input->one("run 10");
+ }
+
+ delete [] x;
+ delete [] v;
+ delete [] type;
+
+ // close down LAMMPS
+
+ delete lmp;
// close down MPI
+ if (lammps == 1) MPI_Comm_free(&comm_lammps);
+ MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
}
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index dbf5fa5219ac04b1469425070e251ee545832f06..78aa2608cfd170cae7c1b159e7a617d9b34d771f 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -1144,9 +1144,9 @@ cvm::real colvar::update_forces_energy()
// For a periodic colvar, both walls may be applicable at the same time
// in which case we pick the closer one
if ( (!is_enabled(f_cv_upper_wall)) ||
- (this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) {
+ (this->dist2(x_reported, lower_wall) < this->dist2(x_reported, upper_wall)) ) {
- cvm::real const grad = this->dist2_lgrad(x, lower_wall);
+ cvm::real const grad = this->dist2_lgrad(x_reported, lower_wall);
if (grad < 0.0) {
fw = -0.5 * lower_wall_k * grad;
f += fw;
@@ -1157,7 +1157,7 @@ cvm::real colvar::update_forces_energy()
} else {
- cvm::real const grad = this->dist2_lgrad(x, upper_wall);
+ cvm::real const grad = this->dist2_lgrad(x_reported, upper_wall);
if (grad > 0.0) {
fw = -0.5 * upper_wall_k * grad;
f += fw;
@@ -1177,15 +1177,21 @@ cvm::real colvar::update_forces_energy()
// atoms only feel the harmonic force
// fr: bias force on extended variable (without harmonic spring), for output in trajectory
// f_ext: total force on extended variable (including harmonic spring)
- // f: - initially, external biasing force
+ // f: - initially, external biasing force (including wall forces)
// - after this code block, colvar force to be applied to atomic coordinates, ie. spring force
fr = f;
f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x);
- // The total force acting on the extended variable is f_ext
- // This will be used in the next timestep
- ft_reported = f_ext;
+ if (is_enabled(f_cv_subtract_applied_force)) {
+ // Report a "system" force without the biases on this colvar
+ // that is, just the spring force
+ ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x);
+ } else {
+ // The total force acting on the extended variable is f_ext
+ // This will be used in the next timestep
+ ft_reported = f_ext;
+ }
// leapfrog: starting from x_i, f_i, v_(i-1/2)
vr += (0.5 * dt) * f_ext / ext_mass;
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index 0d0a339a69382cdae29ceb4a839201a58dc2869a..e8046a97d5d2db92b3048d6106b01da8d2ff0452 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -367,6 +367,7 @@ int cvm::atom_group::parse(std::string const &conf)
cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", INPUT_ERROR);
}
+ // NOTE: calls to add_atom() and/or add_atom_id() are in the proxy-implemented function
cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
}
}
@@ -403,11 +404,21 @@ int cvm::atom_group::parse(std::string const &conf)
}
}
+ // We need to know the fitting options to decide whether the group is scalable
+ parse_error |= parse_fitting_options(group_conf);
+
+ if (is_available(f_ag_scalable_com) && !b_rotate) {
+ enable(f_ag_scalable_com);
+ enable(f_ag_scalable);
+ }
+
if (is_enabled(f_ag_scalable) && !b_dummy) {
+ cvm::log("Enabling scalable calculation for group \""+this->key+"\".\n");
index = (cvm::proxy)->init_atom_group(atoms_ids);
}
- parse_error |= parse_fitting_options(group_conf);
+ bool b_print_atom_ids = false;
+ get_keyval(group_conf, "printAtomIDs", b_print_atom_ids, false, colvarparse::parse_silent);
// TODO move this to colvarparse object
check_keywords(group_conf, key.c_str());
@@ -427,6 +438,10 @@ int cvm::atom_group::parse(std::string const &conf)
cvm::to_str(total_mass)+", total charge = "+
cvm::to_str(total_charge)+".\n");
+ if (b_print_atom_ids) {
+ cvm::log("Internal definition of the atom group:\n");
+ }
+
cvm::decrease_depth();
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@@ -583,6 +598,21 @@ int cvm::atom_group::add_atom_name_residue_range(std::string const &psf_segid,
}
+std::string const cvm::atom_group::print_atom_ids() const
+{
+ size_t line_count = 0;
+ std::ostringstream os("");
+ for (size_t i = 0; i < atoms_ids.size(); i++) {
+ os << " " << std::setw(9) << atoms_ids[i];
+ if (++line_count == 7) {
+ os << "\n";
+ line_count = 0;
+ }
+ }
+ return os.str();
+}
+
+
int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
{
bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false);
@@ -1118,8 +1148,7 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)
log("Communicating a colvar force from atom group to the MD engine.\n");
}
- if (b_dummy)
- return;
+ if (b_dummy) return;
if (noforce) {
cvm::error("Error: sending a force to a group that has "
@@ -1161,17 +1190,21 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)
void cvm::atom_group::apply_force(cvm::rvector const &force)
{
- if (b_dummy)
- return;
+ if (cvm::debug()) {
+ log("Communicating a colvar force from atom group to the MD engine.\n");
+ }
+
+ if (b_dummy) return;
if (noforce) {
cvm::error("Error: sending a force to a group that has "
- "\"disableForces\" defined.\n");
+ "\"enableForces\" set to off.\n");
return;
}
if (is_enabled(f_ag_scalable)) {
(cvm::proxy)->apply_atom_group_force(index, force);
+ return;
}
if (b_rotate) {
diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h
index 7a4d031daa0c5136fa70f1d36bc85de328bc391b..a2a771a8a3f3f8a181c2f05035a826cc68b6a2a9 100644
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@@ -253,6 +253,8 @@ public:
return atoms.size();
}
+ std::string const print_atom_ids() const;
+
/// \brief If this option is on, this group merely acts as a wrapper
/// for a fixed position; any calls to atoms within or to
/// functions that return disaggregated data will fail
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index e3690a4eadfbbc02c37af8ef526864d46df9bcba..1665ea61a9c577afea8cbe3a3d8c680729581aa2 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -80,6 +80,7 @@ int colvarbias_abf::init(std::string const &conf)
if (update_bias) {
// Request calculation of total force (which also checks for availability)
+ // TODO - change this to a dependency - needs ABF-specific features
if(enable(f_cvb_get_total_force)) return cvm::get_error();
}
@@ -133,6 +134,10 @@ int colvarbias_abf::init(std::string const &conf)
// Data for eABF z-based estimator
if (b_extended) {
+ // CZAR output files for stratified eABF
+ get_keyval(conf, "writeCZARwindowFile", b_czar_window_file, false,
+ colvarparse::parse_silent);
+
z_bin.assign(colvars.size(), 0);
z_samples = new colvar_grid_count(colvars);
z_samples->request_actual_value();
@@ -241,7 +246,7 @@ int colvarbias_abf::update()
for (size_t i = 0; i < colvars.size(); i++) {
// get total forces (lagging by 1 timestep) from colvars
- // and subtract previous ABF force
+ // and subtract previous ABF force if necessary
update_system_force(i);
}
gradients->acc_force(force_bin, system_force);
@@ -457,28 +462,30 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
if (z_gradients) {
// Write eABF-related quantities
+
std::string z_samples_out_name = prefix + ".zcount";
- std::string z_gradients_out_name = prefix + ".zgrad";
- std::string czar_gradients_out_name = prefix + ".czar";
cvm::ofstream z_samples_os;
- cvm::ofstream z_gradients_os;
- cvm::ofstream czar_gradients_os;
if (!append) cvm::backup_file(z_samples_out_name.c_str());
z_samples_os.open(z_samples_out_name.c_str(), mode);
if (!z_samples_os.is_open()) {
- cvm::error("Error opening ABF z sample file " + z_samples_out_name + " for writing");
+ cvm::error("Error opening eABF z-histogram file " + z_samples_out_name + " for writing");
}
z_samples->write_multicol(z_samples_os);
z_samples_os.close();
- if (!append) cvm::backup_file(z_gradients_out_name.c_str());
- z_gradients_os.open(z_gradients_out_name.c_str(), mode);
- if (!z_gradients_os.is_open()) {
- cvm::error("Error opening ABF z gradient file " + z_gradients_out_name + " for writing");
+ if (b_czar_window_file) {
+ std::string z_gradients_out_name = prefix + ".zgrad";
+ cvm::ofstream z_gradients_os;
+
+ if (!append) cvm::backup_file(z_gradients_out_name.c_str());
+ z_gradients_os.open(z_gradients_out_name.c_str(), mode);
+ if (!z_gradients_os.is_open()) {
+ cvm::error("Error opening eABF z-gradient file " + z_gradients_out_name + " for writing");
+ }
+ z_gradients->write_multicol(z_gradients_os);
+ z_gradients_os.close();
}
- z_gradients->write_multicol(z_gradients_os);
- z_gradients_os.close();
// Calculate CZAR estimator of gradients
for (std::vector<int> ix = czar_gradients->new_index();
@@ -490,6 +497,9 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
}
}
+ std::string czar_gradients_out_name = prefix + ".czar.grad";
+ cvm::ofstream czar_gradients_os;
+
if (!append) cvm::backup_file(czar_gradients_out_name.c_str());
czar_gradients_os.open(czar_gradients_out_name.c_str(), mode);
if (!czar_gradients_os.is_open()) {
@@ -499,17 +509,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
czar_gradients_os.close();
if (colvars.size() == 1) {
- std::string z_pmf_out_name = prefix + ".zpmf";
- if (!append) cvm::backup_file(z_pmf_out_name.c_str());
- cvm::ofstream z_pmf_os;
- // Do numerical integration and output a PMF
- z_pmf_os.open(z_pmf_out_name.c_str(), mode);
- if (!z_pmf_os.is_open()) cvm::error("Error opening z pmf file " + z_pmf_out_name + " for writing");
- z_gradients->write_1D_integral(z_pmf_os);
- z_pmf_os << std::endl;
- z_pmf_os.close();
-
- std::string czar_pmf_out_name = prefix + ".czarpmf";
+ std::string czar_pmf_out_name = prefix + ".czar.pmf";
if (!append) cvm::backup_file(czar_pmf_out_name.c_str());
cvm::ofstream czar_pmf_os;
// Do numerical integration and output a PMF
@@ -520,8 +520,6 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
czar_pmf_os.close();
}
}
-
-
return;
}
@@ -559,7 +557,7 @@ void colvarbias_abf::read_gradients_samples()
std::ifstream is;
- cvm::log("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name);
+ cvm::log("Reading sample count from " + samples_in_name + " and gradient from " + gradients_in_name);
is.open(samples_in_name.c_str());
if (!is.is_open()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading");
samples->read_multicol(is, true);
@@ -572,17 +570,18 @@ void colvarbias_abf::read_gradients_samples()
is.close();
if (z_gradients) {
- cvm::log("Reading z sample count from " + z_samples_in_name + " and z gradients from " + z_gradients_in_name);
+ // Read eABF z-averaged data for CZAR
+ cvm::log("Reading z-histogram from " + z_samples_in_name + " and z-gradient from " + z_gradients_in_name);
is.clear();
is.open(z_samples_in_name.c_str());
- if (!is.is_open()) cvm::error("Error opening ABF z samples file " + z_samples_in_name + " for reading");
+ if (!is.is_open()) cvm::error("Error opening eABF z-histogram file " + z_samples_in_name + " for reading");
z_samples->read_multicol(is, true);
is.close();
is.clear();
is.open(z_gradients_in_name.c_str());
- if (!is.is_open()) cvm::error("Error opening ABF z gradient file " + z_gradients_in_name + " for reading");
+ if (!is.is_open()) cvm::error("Error opening eABF z-gradient file " + z_gradients_in_name + " for reading");
z_gradients->read_multicol(is, true);
is.close();
}
diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h
index a706bbb5ceb0dd5cc9ad504c2f55de32411b10ec..ee7624096835e52bf27dc02cbb1b96d10c6e9c81 100644
--- a/lib/colvars/colvarbias_abf.h
+++ b/lib/colvars/colvarbias_abf.h
@@ -40,6 +40,8 @@ private:
int output_freq;
/// Write combined files with a history of all output data?
bool b_history_files;
+ /// Write CZAR output file for stratified eABF (.zgrad)
+ bool b_czar_window_file;
size_t history_freq;
/// Cap applied biasing force?
diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp
index 596972253b6dd1cc4eab86d3ee3b5ed534b2b0a2..2aa88bd12a7f408804f23416dc2ba4828dcb13a6 100644
--- a/lib/colvars/colvarbias_meta.cpp
+++ b/lib/colvars/colvarbias_meta.cpp
@@ -159,6 +159,23 @@ int colvarbias_meta::init(std::string const &conf)
cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\n");
}
+
+ // for ebmeta
+ target_dist = NULL;
+ get_keyval(conf, "ebMeta", ebmeta, false);
+ if(ebmeta){
+ target_dist = new colvar_grid_scalar();
+ target_dist->init_from_colvars(colvars);
+ get_keyval(conf, "targetdistfile", target_dist_file);
+ std::ifstream targetdiststream(target_dist_file.c_str());
+ target_dist->read_multicol(targetdiststream);
+ // normalize target distribution and multiply by effective volume = exp(differential entropy)
+ target_dist->multiply_constant(1.0/target_dist->integral());
+ cvm::real volume = std::exp(target_dist->entropy());
+ target_dist->multiply_constant(volume);
+ get_keyval(conf, "ebMetaEquilSteps", ebmeta_equil_steps, 0);
+ }
+
if (cvm::debug())
cvm::log("Done initializing the metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
@@ -186,6 +203,11 @@ colvarbias_meta::~colvarbias_meta()
if (hills_traj_os.is_open())
hills_traj_os.close();
+ if(target_dist) {
+ delete target_dist;
+ target_dist = NULL;
+ }
+
if (cvm::n_meta_biases > 0)
cvm::n_meta_biases -= 1;
}
@@ -221,9 +243,7 @@ colvarbias_meta::create_hill(colvarbias_meta::hill const &h)
// output to trajectory (if specified)
if (hills_traj_os.is_open()) {
hills_traj_os << (hills.back()).output_traj();
- if (cvm::debug()) {
- hills_traj_os.flush();
- }
+ hills_traj_os.flush();
}
has_data = true;
@@ -258,8 +278,7 @@ colvarbias_meta::delete_hill(hill_iter &h)
hills_traj_os << "# DELETED this hill: "
<< (hills.back()).output_traj()
<< "\n";
- if (cvm::debug())
- hills_traj_os.flush();
+ hills_traj_os.flush();
}
return hills.erase(h);
@@ -381,38 +400,37 @@ int colvarbias_meta::update()
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": adding a new hill at step "+cvm::to_str(cvm::step_absolute())+".\n");
- switch (comm) {
+ cvm::real hills_scale=1.0;
- case single_replica:
- if (well_tempered) {
- cvm::real hills_energy_sum_here = 0.0;
- if (use_grids) {
- std::vector<int> curr_bin = hills_energy->get_colvars_index();
- hills_energy_sum_here = hills_energy->value(curr_bin);
- } else {
- calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here);
- }
- cvm::real const exp_weight = std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
- create_hill(hill((hill_weight*exp_weight), colvars, hill_width));
+ if (ebmeta) {
+ hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
+ if(cvm::step_absolute() <= ebmeta_equil_steps) {
+ cvm::real const hills_lambda=(cvm::real(ebmeta_equil_steps - cvm::step_absolute()))/(cvm::real(ebmeta_equil_steps));
+ hills_scale = hills_lambda + (1-hills_lambda)*hills_scale;
+ }
+ }
+
+ if (well_tempered) {
+ cvm::real hills_energy_sum_here = 0.0;
+ if (use_grids) {
+ std::vector<int> curr_bin = hills_energy->get_colvars_index();
+ hills_energy_sum_here = hills_energy->value(curr_bin);
} else {
- create_hill(hill(hill_weight, colvars, hill_width));
+ calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here);
}
+ hills_scale *= std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
+ }
+
+ switch (comm) {
+
+ case single_replica:
+
+ create_hill(hill(hill_weight*hills_scale, colvars, hill_width));
+
break;
case multiple_replicas:
- if (well_tempered) {
- cvm::real hills_energy_sum_here = 0.0;
- if (use_grids) {
- std::vector<int> curr_bin = hills_energy->get_colvars_index();
- hills_energy_sum_here = hills_energy->value(curr_bin);
- } else {
- calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here);
- }
- cvm::real const exp_weight = std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
- create_hill(hill((hill_weight*exp_weight), colvars, hill_width, replica_id));
- } else {
- create_hill(hill(hill_weight, colvars, hill_width, replica_id));
- }
+ create_hill(hill(hill_weight*hills_scale, colvars, hill_width, replica_id));
if (replica_hills_os.is_open()) {
replica_hills_os << hills.back();
} else {
@@ -1507,7 +1525,9 @@ int colvarbias_meta::setup_output()
("."+replica_id) :
("") )+
".hills.traj");
- hills_traj_os.open(traj_file_name.c_str());
+ if (!hills_traj_os.is_open()) {
+ hills_traj_os.open(traj_file_name.c_str());
+ }
if (!hills_traj_os.is_open())
cvm::error("Error: in opening hills output file \"" +
traj_file_name+"\".\n", FILE_ERROR);
diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h
index c021b2ffeba4ce587e1de6b00b48183d3c5e6434..6936e9795454d702a8d825f8f2f0a278f28c1b42 100644
--- a/lib/colvars/colvarbias_meta.h
+++ b/lib/colvars/colvarbias_meta.h
@@ -147,6 +147,14 @@ protected:
/// \brief Biasing temperature in well-tempered metadynamics
cvm::real bias_temperature;
+ // EBmeta parameters
+ bool ebmeta;
+ colvar_grid_scalar* target_dist;
+ std::string target_dist_file;
+ cvm::real target_dist_volume;
+ size_t ebmeta_equil_steps;
+
+
/// \brief Try to read the restart information by allocating new
/// grids before replacing the current ones (used e.g. in
/// multiple_replicas)
diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index 73b1dbd81c132b5bcfc4f7f88b30933839f2096a..050e03f78ed520d7d4068a16734c28eb24ebffc4 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -84,15 +84,12 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
if (is_available(f_cvc_scalable_com) && is_available(f_cvc_com_based)) {
enable(f_cvc_scalable_com);
enable(f_cvc_scalable);
- group->enable(f_ag_scalable_com);
- group->enable(f_ag_scalable);
+ // The CVC makes the feature available;
+ // the atom group will enable it unless it needs to compute a rotational fit
+ group->provide(f_ag_scalable_com);
}
// TODO check for other types of parallelism here
-
- if (is_enabled(f_cvc_scalable)) {
- cvm::log("Will enable scalable calculation for group \""+group->key+"\".\n");
- }
}
if (group->parse(conf) == COLVARS_OK) {
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index 7736ea14bd1f5f66b786050f3d9102a365a3be6d..a44f86c2925661d410f4cc6a321ed1ccbdcf4f7b 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -449,8 +449,10 @@ void colvardeps::init_ag_requires() {
// Features that are implemented (or not) by all atom groups
feature_states[f_ag_active]->available = true;
- feature_states[f_ag_scalable_com]->available = (cvm::proxy->scalable_group_coms() == COLVARS_OK);
- feature_states[f_ag_scalable]->available = feature_states[f_ag_scalable_com]->available;
+ // f_ag_scalable_com is provided by the CVC iff it is COM-based
+ feature_states[f_ag_scalable_com]->available = false;
+ // TODO make f_ag_scalable depend on f_ag_scalable_com (or something else)
+ feature_states[f_ag_scalable]->available = true;
}
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index ee5db06e300b6e006cd5a0a9fd19bff6d49a3c28..8e110275b21501208c725c270175ff46e1f11813 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -156,6 +156,12 @@ int colvarmodule::parse_global_params(std::string const &conf)
read_index_file(index_file_name.c_str());
}
+ if (parse->get_keyval(conf, "smp", proxy->b_smp_active, proxy->b_smp_active)) {
+ if (proxy->b_smp_active == false) {
+ cvm::log("SMP parallelism has been disabled.\n");
+ }
+ }
+
parse->get_keyval(conf, "analysis", b_analysis, b_analysis);
parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size,
@@ -810,6 +816,7 @@ int colvarmodule::analyze()
int colvarmodule::setup()
{
+ if (this->size() == 0) return cvm::get_error();
// loop over all components of all colvars to reset masses of all groups
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end(); cvi++) {
@@ -867,25 +874,35 @@ int colvarmodule::reset()
int colvarmodule::setup_input()
{
- // name of input state file
- restart_in_name = proxy->input_prefix().size() ?
- std::string(proxy->input_prefix()+".colvars.state") :
- std::string("") ;
+ if (this->size() == 0) return cvm::get_error();
+
+ std::string restart_in_name("");
// read the restart configuration, if available
- if (restart_in_name.size()) {
+ if (proxy->input_prefix().size()) {
// read the restart file
+ restart_in_name = proxy->input_prefix();
std::ifstream input_is(restart_in_name.c_str());
if (!input_is.good()) {
- cvm::error("Error: in opening restart file \""+
+ // try by adding the suffix
+ input_is.clear();
+ restart_in_name = restart_in_name+std::string(".colvars.state");
+ input_is.open(restart_in_name.c_str());
+ }
+
+ if (!input_is.good()) {
+ cvm::error("Error: in opening input file \""+
std::string(restart_in_name)+"\".\n",
FILE_ERROR);
return COLVARS_ERROR;
} else {
+ cvm::log(cvm::line_marker);
cvm::log("Restarting from file \""+restart_in_name+"\".\n");
read_restart(input_is);
if (cvm::get_error() != COLVARS_OK) {
return COLVARS_ERROR;
+ } else {
+ proxy->input_prefix().clear();
}
cvm::log(cvm::line_marker);
}
@@ -897,7 +914,9 @@ int colvarmodule::setup_input()
int colvarmodule::setup_output()
{
- int error_code = 0;
+ if (this->size() == 0) return cvm::get_error();
+
+ int error_code = COLVARS_OK;
// output state file (restart)
restart_out_name = proxy->restart_output_prefix().size() ?
@@ -1545,11 +1564,7 @@ std::list<std::string> colvarmodule::index_group_names;
std::list<std::vector<int> > colvarmodule::index_groups;
bool colvarmodule::use_scripted_forces = false;
bool colvarmodule::scripting_after_biases = true;
-
-// file name prefixes
std::string colvarmodule::output_prefix = "";
-std::string colvarmodule::restart_in_name = "";
-
// i/o constants
size_t const colvarmodule::it_width = 12;
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 0b19d1f779dc95a404a01fe8266bcc9fb4e62ccc..c8df6c7b5eb9f54dc1e752a7fb3c72559c95ac19 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -4,7 +4,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-10-05"
+#define COLVARS_VERSION "2016-10-21"
#endif
#ifndef COLVARS_DEBUG
@@ -12,6 +12,9 @@
#endif
/*! \mainpage Main page
+This is the Developer's documentation for the Collective Variables Module.
+
+You can browse the class hierarchy or the list of source files.
*/
/// \file colvarmodule.h
@@ -154,9 +157,6 @@ public:
/// Prefix for all output files for this run
static std::string output_prefix;
- /// input restart file name (determined from input_prefix)
- static std::string restart_in_name;
-
/// Array of collective variables
static std::vector<colvar *> colvars;
@@ -197,6 +197,11 @@ public:
return COLVARS_DEBUG;
}
+ /// \brief How many objects are configured yet?
+ inline size_t const size() const
+ {
+ return colvars.size() + biases.size();
+ }
/// \brief Constructor \param config_name Configuration file name
/// \param restart_name (optional) Restart file name
diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp
index e222fa9ecf705b6cca4a4942527885380e79dd0e..86cd00ad29e504138f05e3959ca03fb38e62f6da 100644
--- a/lib/colvars/colvarparse.cpp
+++ b/lib/colvars/colvarparse.cpp
@@ -644,9 +644,9 @@ bool colvarparse::key_lookup(std::string const &conf,
// find the matching closing brace
- if (cvm::debug()) {
- cvm::log("Multi-line value, config is now \""+line+"\".\n");
- }
+// if (cvm::debug()) {
+// cvm::log("Multi-line value, config is now \""+line+"\".\n");
+// }
int brace_count = 1;
@@ -689,9 +689,9 @@ bool colvarparse::key_lookup(std::string const &conf,
line_end = nl;
line.append(conf, line_begin, (line_end-line_begin));
- if (cvm::debug()) {
- cvm::log("Added a new line, config is now \""+line+"\".\n");
- }
+// if (cvm::debug()) {
+// cvm::log("Added a new line, config is now \""+line+"\".\n");
+// }
}
if (brace_count < 0) {
diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h
index 3d41cb58a7f3107b639de5ed8f10c12b9467e38f..59f6dc9bb21a57f02ced60741752c147e779aa84 100644
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@@ -25,7 +25,7 @@ public:
colvarmodule *colvars;
/// Default constructor
- inline colvarproxy() : script(NULL) {}
+ inline colvarproxy() : script(NULL), b_smp_active(true) {}
/// Default destructor
virtual ~colvarproxy() {}
@@ -80,7 +80,7 @@ public:
/// configuration)
std::string input_prefix_str, output_prefix_str, restart_output_prefix_str;
- inline std::string input_prefix()
+ inline std::string & input_prefix()
{
return input_prefix_str;
}
@@ -116,12 +116,15 @@ public:
// ***************** SHARED-MEMORY PARALLELIZATION *****************
- /// Whether or not threaded parallelization is available
+ /// Whether threaded parallelization is available (TODO: make this a cvm::deps feature)
virtual int smp_enabled()
{
return COLVARS_NOT_IMPLEMENTED;
}
+ /// Whether threaded parallelization should be used (TODO: make this a cvm::deps feature)
+ bool b_smp_active;
+
/// Distribute calculation of colvars (and their components) across threads
virtual int smp_colvars_loop()
{
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index 86b52ac819824a1496f1f15af2790a5a2bcb5bbd..a0269aa7a911d8a1f3321e8df9bd40a6283308b0 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -14,6 +14,36 @@ colvarscript::colvarscript(colvarproxy *p)
{
}
+
+extern "C" {
+
+ // Generic hooks; NAMD and VMD have Tcl-specific versions in the respective proxies
+
+ int run_colvarscript_command(int argc, const char **argv)
+ {
+ colvarproxy *cvp = cvm::proxy;
+ if (!cvp) {
+ return -1;
+ }
+ if (!cvp->script) {
+ cvm::error("Called run_colvarscript_command without a script object initialized.\n");
+ return -1;
+ }
+ return cvp->script->run(argc, argv);
+ }
+
+ const char * get_colvarscript_result()
+ {
+ colvarproxy *cvp = cvm::proxy;
+ if (!cvp->script) {
+ cvm::error("Called run_colvarscript_command without a script object initialized.\n");
+ return "";
+ }
+ return cvp->script->result.c_str();
+ }
+}
+
+
/// Run method based on given arguments
int colvarscript::run(int argc, char const *argv[]) {
@@ -125,7 +155,7 @@ int colvarscript::run(int argc, char const *argv[]) {
result = "Missing arguments\n" + help_string();
return COLVARSCRIPT_ERROR;
}
- proxy->input_prefix_str = argv[2];
+ proxy->input_prefix() = argv[2];
if (colvars->setup_input() == COLVARS_OK) {
return COLVARS_OK;
} else {
diff --git a/python/examples/simple.py b/python/examples/simple.py
index 9315a555fca0e64c368b2a8d9ecf8576a779a710..6336c89de3e6171d16ed82968a99ed2bf8d39a3d 100755
--- a/python/examples/simple.py
+++ b/python/examples/simple.py
@@ -13,6 +13,8 @@
from __future__ import print_function
import sys
+import numpy as np
+import ctypes
# parse command line
@@ -51,17 +53,39 @@ for line in lines: lmp.command(line)
lmp.command("run 10")
x = lmp.gather_atoms("x",1,3)
+v = lmp.gather_atoms("v",1,3)
epsilon = 0.1
x[0] += epsilon
lmp.scatter_atoms("x",1,3,x)
lmp.command("run 1");
+# extract force on single atom two different ways
+
f = lmp.extract_atom("f",3)
print("Force on 1 atom via extract_atom: ",f[0][0])
fx = lmp.extract_variable("fx","all",1)
print("Force on 1 atom via extract_variable:",fx[0])
+# use commands_string() and commands_list() to invoke more commands
+
+strtwo = "run 10\nrun 20"
+lmp.commands_string(strtwo)
+
+cmds = ["run 10","run 20"]
+lmp.commands_list(cmds)
+
+# delete all atoms
+# create_atoms() to create new ones with old coords, vels
+# initial thermo should be same as step 20
+
+natoms = lmp.get_natoms()
+type = natoms*[1]
+
+lmp.command("delete_atoms group all");
+lmp.create_atoms(natoms,None,type,x,v);
+lmp.command("run 10");
+
# uncomment if running in parallel via Pypar
#print("Proc %d out of %d procs has" % (me,nprocs), lmp)
#pypar.finalize()
diff --git a/python/lammps.py b/python/lammps.py
index 8c5792806450e6d73119531d41ca0cd5d1a33f6e..7344a16be30fb8ab47e3ffc1f6cd939dae95b4fe 100644
--- a/python/lammps.py
+++ b/python/lammps.py
@@ -172,6 +172,8 @@ class lammps(object):
if file: file = file.encode()
self.lib.lammps_file(self.lmp,file)
+ # send a single command
+
def command(self,cmd):
if cmd: cmd = cmd.encode()
self.lib.lammps_command(self.lmp,cmd)
@@ -185,6 +187,19 @@ class lammps(object):
raise MPIAbortException(error_msg)
raise Exception(error_msg)
+ # send a list of commands
+
+ def commands_list(self,cmdlist):
+ args = (c_char_p * len(cmdlist))(*cmdlist)
+ self.lib.lammps_commands_list(self.lmp,len(cmdlist),args)
+
+ # send a string of commands
+
+ def commands_string(self,multicmd):
+ self.lib.lammps_commands_string(self.lmp,c_char_p(multicmd))
+
+ # extract global info
+
def extract_global(self,name,type):
if name: name = name.encode()
if type == 0:
@@ -195,6 +210,8 @@ class lammps(object):
ptr = self.lib.lammps_extract_global(self.lmp,name)
return ptr[0]
+ # extract per-atom info
+
def extract_atom(self,name,type):
if name: name = name.encode()
if type == 0:
@@ -209,6 +226,8 @@ class lammps(object):
ptr = self.lib.lammps_extract_atom(self.lmp,name)
return ptr
+ # extract compute info
+
def extract_compute(self,id,style,type):
if id: id = id.encode()
if type == 0:
@@ -226,6 +245,7 @@ class lammps(object):
return ptr
return None
+ # extract fix info
# in case of global datum, free memory for 1 double via lammps_free()
# double was allocated by library interface function
@@ -249,6 +269,7 @@ class lammps(object):
else:
return None
+ # extract variable info
# free memory for 1 double or 1 vector of doubles via lammps_free()
# for vector, must copy nlocal returned values to local c_double vector
# memory was allocated by library interface function
@@ -296,7 +317,14 @@ class lammps(object):
return self.lib.lammps_get_natoms(self.lmp)
# return vector of atom properties gathered across procs, ordered by atom ID
-
+ # name = atom property recognized by LAMMPS in atom->extract()
+ # type = 0 for integer values, 1 for double values
+ # count = number of per-atom valus, 1 for type or charge, 3 for x or f
+ # returned data is a 1d vector - doc how it is ordered?
+ # NOTE: how could we insure are converting to correct Python type
+ # e.g. for Python list or NumPy, etc
+ # ditto for extact_atom() above
+
def gather_atoms(self,name,type,count):
if name: name = name.encode()
natoms = self.lib.lammps_get_natoms(self.lmp)
@@ -310,12 +338,38 @@ class lammps(object):
return data
# scatter vector of atom properties across procs, ordered by atom ID
- # assume vector is of correct type and length, as created by gather_atoms()
-
+ # name = atom property recognized by LAMMPS in atom->extract()
+ # type = 0 for integer values, 1 for double values
+ # count = number of per-atom valus, 1 for type or charge, 3 for x or f
+ # assume data is of correct type and length, as created by gather_atoms()
+ # NOTE: how could we insure are passing correct type to LAMMPS
+ # e.g. for Python list or NumPy, etc
+
def scatter_atoms(self,name,type,count,data):
if name: name = name.encode()
self.lib.lammps_scatter_atoms(self.lmp,name,type,count,data)
+ # create N atoms on all procs
+ # N = global number of atoms
+ # id = ID of each atom (optional, can be None)
+ # type = type of each atom (1 to Ntypes) (required)
+ # x = coords of each atom as (N,3) array (required)
+ # v = velocity of each atom as (N,3) array (optional, can be None)
+ # NOTE: how could we insure are passing correct type to LAMMPS
+ # e.g. for Python list or NumPy, etc
+ # ditto for gather_atoms() above
+
+ def create_atoms(self,n,id,type,x,v):
+ if id:
+ id_lmp = (c_int * n)()
+ id_lmp[:] = id
+ else: id_lmp = id
+ type_lmp = (c_int * n)()
+ type_lmp[:] = type
+ self.lib.lammps_create_atoms(self.lmp,n,id_lmp,type_lmp,x,v)
+
+ # document this?
+
@property
def uses_exceptions(self):
try:
diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh
index c45be8c973502b9fa536dc6ba7665bb56748184f..41f9304f8bf6cdd6d1f76e83589334336efe9094 100644
--- a/src/KOKKOS/Install.sh
+++ b/src/KOKKOS/Install.sh
@@ -83,6 +83,10 @@ action fix_nvt_kokkos.cpp
action fix_nvt_kokkos.h
action fix_qeq_reax_kokkos.cpp fix_qeq_reax.cpp
action fix_qeq_reax_kokkos.h fix_qeq_reax.h
+action fix_reaxc_bonds_kokkos.cpp fix_reaxc_bonds.cpp
+action fix_reaxc_bonds_kokkos.h fix_reaxc_bonds.h
+action fix_reaxc_species_kokkos.cpp fix_reaxc_species.cpp
+action fix_reaxc_species_kokkos.h fix_reaxc_species.h
action fix_setforce_kokkos.cpp
action fix_setforce_kokkos.h
action fix_wall_reflect_kokkos.cpp
diff --git a/src/KOKKOS/fix_qeq_reax_kokkos.cpp b/src/KOKKOS/fix_qeq_reax_kokkos.cpp
index e54be7412402736e71243395442839753c01bb95..d39d5cc84a2b8ff322fc411a3246c364bfe6c6f1 100644
--- a/src/KOKKOS/fix_qeq_reax_kokkos.cpp
+++ b/src/KOKKOS/fix_qeq_reax_kokkos.cpp
@@ -431,8 +431,8 @@ void FixQEqReaxKokkos<DeviceType>::compute_h_item(int ii, int &m_fill, const boo
const F_FLOAT rsq = delx*delx + dely*dely + delz*delz;
if (rsq > cutsq) continue;
- const F_FLOAT r = sqrt(rsq);
if (final) {
+ const F_FLOAT r = sqrt(rsq);
d_jlist(m_fill) = j;
const F_FLOAT shldij = d_shield(itype,jtype);
d_val(m_fill) = calculate_H_k(r,shldij);
diff --git a/src/KOKKOS/fix_reaxc_bonds_kokkos.cpp b/src/KOKKOS/fix_reaxc_bonds_kokkos.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6e96c491968d6d90056de3c07b8c4e56ca652b0f
--- /dev/null
+++ b/src/KOKKOS/fix_reaxc_bonds_kokkos.cpp
@@ -0,0 +1,124 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing author: Stan Moore (Sandia)
+------------------------------------------------------------------------- */
+
+#include <stdlib.h>
+#include <string.h>
+#include "fix_ave_atom.h"
+#include "fix_reaxc_bonds_kokkos.h"
+#include "atom.h"
+#include "update.h"
+#include "pair_reax_c_kokkos.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "comm.h"
+#include "force.h"
+#include "compute.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+#include "reaxc_list.h"
+#include "reaxc_types.h"
+#include "reaxc_defs.h"
+#include "atom_masks.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixReaxCBondsKokkos::FixReaxCBondsKokkos(LAMMPS *lmp, int narg, char **arg) :
+ FixReaxCBonds(lmp, narg, arg)
+{
+ kokkosable = 1;
+ atomKK = (AtomKokkos *) atom;
+
+ datamask_read = EMPTY_MASK;
+ datamask_modify = EMPTY_MASK;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixReaxCBondsKokkos::~FixReaxCBondsKokkos()
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCBondsKokkos::init()
+{
+ Pair *pair_kk = force->pair_match("reax/c/kk",1);
+ if (pair_kk == NULL) error->all(FLERR,"Cannot use fix reax/c/bonds without "
+ "pair_style reax/c/kk");
+
+ FixReaxCBonds::init();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCBondsKokkos::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp)
+
+{
+ int i, j;
+ int nbuf_local;
+ int nlocal_max, numbonds, numbonds_max;
+ double *buf;
+ DAT::tdual_ffloat_1d k_buf;
+
+ int nlocal = atom->nlocal;
+ int nlocal_tot = static_cast<int> (atom->natoms);
+
+ numbonds = 0;
+ if (reaxc->execution_space == Device)
+ ((PairReaxCKokkos<LMPDeviceType>*) reaxc)->FindBond(numbonds);
+ else
+ ((PairReaxCKokkos<LMPHostType>*) reaxc)->FindBond(numbonds);
+
+ // allocate a temporary buffer for the snapshot info
+ MPI_Allreduce(&numbonds,&numbonds_max,1,MPI_INT,MPI_MAX,world);
+ MPI_Allreduce(&nlocal,&nlocal_max,1,MPI_INT,MPI_MAX,world);
+
+ nbuf = 1+(numbonds_max*2+10)*nlocal_max;
+ memory->create_kokkos(k_buf,buf,nbuf,"reax/c/bonds:buf");
+
+ // Pass information to buffer
+ if (reaxc->execution_space == Device)
+ ((PairReaxCKokkos<LMPDeviceType>*) reaxc)->PackBondBuffer(k_buf,nbuf_local);
+ else
+ ((PairReaxCKokkos<LMPHostType>*) reaxc)->PackBondBuffer(k_buf,nbuf_local);
+
+ // Receive information from buffer for output
+ RecvBuffer(buf, nbuf, nbuf_local, nlocal_tot, numbonds_max);
+
+ memory->destroy_kokkos(k_buf,buf);
+}
+
+double FixReaxCBondsKokkos::memory_usage()
+{
+ double bytes;
+
+ bytes = nbuf*sizeof(double);
+ // These are accounted for in PairReaxCKokkos:
+ //bytes += nmax*sizeof(int);
+ //bytes += 1.0*nmax*MAXREAXBOND*sizeof(double);
+ //bytes += 1.0*nmax*MAXREAXBOND*sizeof(int);
+
+ return bytes;
+}
diff --git a/src/KOKKOS/fix_reaxc_bonds_kokkos.h b/src/KOKKOS/fix_reaxc_bonds_kokkos.h
new file mode 100644
index 0000000000000000000000000000000000000000..f27e558fddcd835efe3171a9127f462f6a9d435e
--- /dev/null
+++ b/src/KOKKOS/fix_reaxc_bonds_kokkos.h
@@ -0,0 +1,42 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(reax/c/bonds/kk,FixReaxCBondsKokkos)
+
+#else
+
+#ifndef LMP_FIX_REAXC_BONDS_KOKKOS_H
+#define LMP_FIX_REAXC_BONDS_KOKKOS_H
+
+#include "fix_reaxc_bonds.h"
+#include "kokkos_type.h"
+
+namespace LAMMPS_NS {
+
+class FixReaxCBondsKokkos : public FixReaxCBonds {
+ public:
+ FixReaxCBondsKokkos(class LAMMPS *, int, char **);
+ virtual ~FixReaxCBondsKokkos();
+ void init();
+
+ private:
+ double nbuf;
+ void Output_ReaxC_Bonds(bigint, FILE *);
+ double memory_usage();
+};
+}
+
+#endif
+#endif
diff --git a/src/KOKKOS/fix_reaxc_species_kokkos.cpp b/src/KOKKOS/fix_reaxc_species_kokkos.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d3de0c998c31b937cdebb2c4876bcd972879feb7
--- /dev/null
+++ b/src/KOKKOS/fix_reaxc_species_kokkos.cpp
@@ -0,0 +1,157 @@
+/* ----------------------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+ Contributing authors: Stan Moore (Sandia)
+------------------------------------------------------------------------- */
+
+#include <stdlib.h>
+#include <math.h>
+#include "atom.h"
+#include <string.h>
+#include "fix_ave_atom.h"
+#include "fix_reaxc_species_kokkos.h"
+#include "domain.h"
+#include "update.h"
+#include "pair_reax_c_kokkos.h"
+#include "modify.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.h"
+#include "comm.h"
+#include "force.h"
+#include "compute.h"
+#include "input.h"
+#include "variable.h"
+#include "memory.h"
+#include "error.h"
+#include "reaxc_list.h"
+#include "atom_masks.h"
+
+using namespace LAMMPS_NS;
+using namespace FixConst;
+
+/* ---------------------------------------------------------------------- */
+
+FixReaxCSpeciesKokkos::FixReaxCSpeciesKokkos(LAMMPS *lmp, int narg, char **arg) :
+ FixReaxCSpecies(lmp, narg, arg)
+{
+ kokkosable = 1;
+ atomKK = (AtomKokkos *) atom;
+
+ // NOTE: Could improve performance if a Kokkos version of ComputeSpecAtom is added
+
+ datamask_read = X_MASK | V_MASK | Q_MASK | MASK_MASK;
+ datamask_modify = EMPTY_MASK;
+}
+
+/* ---------------------------------------------------------------------- */
+
+FixReaxCSpeciesKokkos::~FixReaxCSpeciesKokkos()
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCSpeciesKokkos::init()
+{
+ Pair* pair_kk = force->pair_match("reax/c/kk",1);
+ if (pair_kk == NULL) error->all(FLERR,"Cannot use fix reax/c/species/kk without "
+ "pair_style reax/c/kk");
+
+ FixReaxCSpecies::init();
+
+ int irequest = neighbor->request(this,instance_me);
+ neighbor->requests[irequest]->pair = 0;
+ neighbor->requests[irequest]->fix = 1;
+ neighbor->requests[irequest]->newton = 2;
+ neighbor->requests[irequest]->ghost = 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void FixReaxCSpeciesKokkos::FindMolecule()
+{
+ int i,j,ii,jj,inum,itype,jtype,loop,looptot;
+ int change,done,anychange;
+ int *mask = atom->mask;
+ int *ilist;
+ double bo_tmp,bo_cut;
+ double **spec_atom = f_SPECBOND->array_atom;
+
+ inum = list->inum;
+ ilist = list->ilist;
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ if (mask[i] & groupbit) {
+ clusterID[i] = atom->tag[i];
+ x0[i].x = spec_atom[i][1];
+ x0[i].y = spec_atom[i][2];
+ x0[i].z = spec_atom[i][3];
+ }
+ else clusterID[i] = 0.0;
+ }
+
+ loop = 0;
+ while (1) {
+ comm->forward_comm_fix(this);
+ loop ++;
+
+ change = 0;
+ while (1) {
+ done = 1;
+
+ for (ii = 0; ii < inum; ii++) {
+ i = ilist[ii];
+ if (!(mask[i] & groupbit)) continue;
+
+ itype = atom->type[i];
+
+ for (jj = 0; jj < MAXSPECBOND; jj++) {
+ j = reaxc->tmpid[i][jj];
+
+ if (j < i) continue;
+ if (!(mask[j] & groupbit)) continue;
+
+ if (clusterID[i] == clusterID[j] && PBCconnected[i] == PBCconnected[j]
+ && x0[i].x == x0[j].x && x0[i].y == x0[j].y && x0[i].z == x0[j].z) continue;
+
+ jtype = atom->type[j];
+ bo_cut = BOCut[itype][jtype];
+ bo_tmp = spec_atom[i][jj+7];
+
+ if (bo_tmp > bo_cut) {
+ clusterID[i] = clusterID[j] = MIN(clusterID[i], clusterID[j]);
+ PBCconnected[i] = PBCconnected[j] = MAX(PBCconnected[i], PBCconnected[j]);
+ x0[i] = x0[j] = chAnchor(x0[i], x0[j]);
+ if ((fabs(spec_atom[i][1] - spec_atom[j][1]) > reaxc->control->bond_cut)
+ || (fabs(spec_atom[i][2] - spec_atom[j][2]) > reaxc->control->bond_cut)
+ || (fabs(spec_atom[i][3] - spec_atom[j][3]) > reaxc->control->bond_cut))
+ PBCconnected[i] = PBCconnected[j] = 1;
+ done = 0;
+ }
+ }
+ }
+ if (!done) change = 1;
+ if (done) break;
+ }
+ MPI_Allreduce(&change,&anychange,1,MPI_INT,MPI_MAX,world);
+ if (!anychange) break;
+
+ MPI_Allreduce(&loop,&looptot,1,MPI_INT,MPI_SUM,world);
+ if (looptot >= 400*nprocs) break;
+
+ }
+}
\ No newline at end of file
diff --git a/src/KOKKOS/fix_reaxc_species_kokkos.h b/src/KOKKOS/fix_reaxc_species_kokkos.h
new file mode 100644
index 0000000000000000000000000000000000000000..64de060006caf0a9965f40c33bad609cfdd468e2
--- /dev/null
+++ b/src/KOKKOS/fix_reaxc_species_kokkos.h
@@ -0,0 +1,41 @@
+/* -*- c++ -*- ----------------------------------------------------------
+ LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+ http://lammps.sandia.gov, Sandia National Laboratories
+ Steve Plimpton, sjplimp@sandia.gov
+
+ Copyright (2003) Sandia Corporation. Under the terms of Contract
+ DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+ certain rights in this software. This software is distributed under
+ the GNU General Public License.
+
+ See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef FIX_CLASS
+
+FixStyle(reax/c/species/kk,FixReaxCSpeciesKokkos)
+
+#else
+
+#ifndef LMP_FIX_REAXC_SPECIES_KOKKOS_H
+#define LMP_FIX_REAXC_SPECIES_KOKKOS_H
+
+#include "fix_reaxc_species.h"
+
+#define BUFLEN 1000
+
+namespace LAMMPS_NS {
+
+class FixReaxCSpeciesKokkos : public FixReaxCSpecies {
+ public:
+ FixReaxCSpeciesKokkos(class LAMMPS *, int, char **);
+ virtual ~FixReaxCSpeciesKokkos();
+ void init();
+
+ private:
+ void FindMolecule();
+};
+}
+
+#endif
+#endif
diff --git a/src/KOKKOS/neighbor_kokkos.cpp b/src/KOKKOS/neighbor_kokkos.cpp
index 3bb8b0dce4a20bfea8f88b4c7d5cd8b4facb32dd..9622129e703c13f1fdea34c56bb799a2c71aeca2 100644
--- a/src/KOKKOS/neighbor_kokkos.cpp
+++ b/src/KOKKOS/neighbor_kokkos.cpp
@@ -181,6 +181,9 @@ int NeighborKokkos::init_lists_kokkos()
void NeighborKokkos::init_list_flags1_kokkos(int i)
{
+ if (style != BIN)
+ error->all(FLERR,"KOKKOS package only supports 'bin' neighbor lists");
+
if (lists_host[i]) {
lists_host[i]->buildflag = 1;
if (pair_build_host[i] == NULL) lists_host[i]->buildflag = 0;
diff --git a/src/KOKKOS/neighbor_kokkos.h b/src/KOKKOS/neighbor_kokkos.h
index 3693ebd9a475df28ca450463e09fdc7eacd8eb66..8c097139a76afa1d9164ff7aaba40aa1601610ef 100644
--- a/src/KOKKOS/neighbor_kokkos.h
+++ b/src/KOKKOS/neighbor_kokkos.h
@@ -434,6 +434,10 @@ class NeighborKokkos : public Neighbor {
/* ERROR/WARNING messages:
+E: KOKKOS package only supports 'bin' neighbor lists
+
+Self-explanatory.
+
E: Too many local+ghost atoms for neighbor list
The number of nlocal + nghost atoms on a processor
diff --git a/src/KOKKOS/pair_kokkos.h b/src/KOKKOS/pair_kokkos.h
index b880654fbf64933cdf2f7d76f0f37f785c503d83..3710c460c0246e12e10571aa5a9680026597eeb7 100644
--- a/src/KOKKOS/pair_kokkos.h
+++ b/src/KOKKOS/pair_kokkos.h
@@ -611,7 +611,7 @@ struct PairComputeFunctor<PairStyle,N2,STACKPARAMS,Specialisation> {
// pair_compute_neighlist and pair_compute_fullcluster will match - either the dummy version
// or the real one further below.
template<class PairStyle, unsigned NEIGHFLAG, class Specialisation>
-EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if<!(NEIGHFLAG&PairStyle::EnabledNeighFlags), NeighListKokkos<typename PairStyle::device_type>*>::type list) {
+EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if<!((NEIGHFLAG&PairStyle::EnabledNeighFlags) != 0), NeighListKokkos<typename PairStyle::device_type>*>::type list) {
EV_FLOAT ev;
(void) fpair;
(void) list;
@@ -620,7 +620,7 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable
}
template<class PairStyle, class Specialisation>
-EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if<!(FULLCLUSTER&PairStyle::EnabledNeighFlags), NeighListKokkos<typename PairStyle::device_type>*>::type list) {
+EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if<!((FULLCLUSTER&PairStyle::EnabledNeighFlags) != 0), NeighListKokkos<typename PairStyle::device_type>*>::type list) {
EV_FLOAT ev;
(void) fpair;
(void) list;
@@ -630,7 +630,7 @@ EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enab
// Submit ParallelFor for NEIGHFLAG=HALF,HALFTHREAD,FULL,N2
template<class PairStyle, unsigned NEIGHFLAG, class Specialisation>
-EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if<NEIGHFLAG&PairStyle::EnabledNeighFlags, NeighListKokkos<typename PairStyle::device_type>*>::type list) {
+EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if<(NEIGHFLAG&PairStyle::EnabledNeighFlags) != 0, NeighListKokkos<typename PairStyle::device_type>*>::type list) {
EV_FLOAT ev;
if(fpair->atom->ntypes > MAX_TYPES_STACKPARAMS) {
PairComputeFunctor<PairStyle,NEIGHFLAG,false,Specialisation > ff(fpair,list);
@@ -646,7 +646,7 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable
// Submit ParallelFor for NEIGHFLAG=FULLCLUSTER
template<class PairStyle, class Specialisation>
-EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if<FULLCLUSTER&PairStyle::EnabledNeighFlags, NeighListKokkos<typename PairStyle::device_type>*>::type list) {
+EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if<(FULLCLUSTER&PairStyle::EnabledNeighFlags) != 0, NeighListKokkos<typename PairStyle::device_type>*>::type list) {
EV_FLOAT ev;
if(fpair->atom->ntypes > MAX_TYPES_STACKPARAMS) {
typedef PairComputeFunctor<PairStyle,FULLCLUSTER,false,Specialisation >
diff --git a/src/KOKKOS/pair_reax_c_kokkos.cpp b/src/KOKKOS/pair_reax_c_kokkos.cpp
index 4fcb3652e30f2c9881d2bb6e352e3bac585bd761..894c3ab53c9235a7c9d8cbf5e287f7883f536bc0 100644
--- a/src/KOKKOS/pair_reax_c_kokkos.cpp
+++ b/src/KOKKOS/pair_reax_c_kokkos.cpp
@@ -44,7 +44,7 @@
/* ---------------------------------------------------------------------- */
-namespace LAMMPS_NS{
+namespace LAMMPS_NS {
using namespace MathConst;
using namespace MathSpecial;
@@ -69,6 +69,9 @@ PairReaxCKokkos<DeviceType>::PairReaxCKokkos(LAMMPS *lmp) : PairReaxC(lmp)
nmax = 0;
maxbo = 1;
maxhb = 1;
+
+ k_error_flag = DAT::tdual_int_scalar("pair:error_flag");
+ k_nbuf_local = DAT::tdual_int_scalar("pair:nbuf_local");
}
/* ---------------------------------------------------------------------- */
@@ -76,10 +79,15 @@ PairReaxCKokkos<DeviceType>::PairReaxCKokkos(LAMMPS *lmp) : PairReaxC(lmp)
template<class DeviceType>
PairReaxCKokkos<DeviceType>::~PairReaxCKokkos()
{
- if (!copymode) {
- memory->destroy_kokkos(k_eatom,eatom);
- memory->destroy_kokkos(k_vatom,vatom);
- }
+ if (copymode) return;
+
+ memory->destroy_kokkos(k_eatom,eatom);
+ memory->destroy_kokkos(k_vatom,vatom);
+
+ memory->destroy_kokkos(k_tmpid,tmpid);
+ tmpid = NULL;
+ memory->destroy_kokkos(k_tmpbo,tmpbo);
+ tmpbo = NULL;
}
/* ---------------------------------------------------------------------- */
@@ -306,6 +314,11 @@ void PairReaxCKokkos<DeviceType>::setup()
bo_cut = 0.01 * gp[29];
thb_cut = control->thb_cut;
thb_cutsq = 0.000010; //thb_cut*thb_cut;
+
+ if (atom->nmax > nmax) {
+ nmax = atom->nmax;
+ allocate_array();
+ }
}
/* ---------------------------------------------------------------------- */
@@ -980,6 +993,9 @@ void PairReaxCKokkos<DeviceType>::compute(int eflag_in, int vflag_in)
k_vatom.template sync<LMPHostType>();
}
+ if (fixspecies_flag)
+ FindBondSpecies();
+
copymode = 0;
}
@@ -1346,6 +1362,19 @@ void PairReaxCKokkos<DeviceType>::allocate_array()
d_Delta_lp_temp = typename AT::t_ffloat_1d("reax/c/kk:Delta_lp_temp",nmax);
d_CdDelta = typename AT::t_ffloat_1d("reax/c/kk:CdDelta",nmax);
d_sum_ovun = typename AT::t_ffloat_2d_dl("reax/c/kk:sum_ovun",nmax,3);
+
+ // FixReaxCSpecies
+ if (fixspecies_flag) {
+ memory->destroy_kokkos(k_tmpid,tmpid);
+ memory->destroy_kokkos(k_tmpbo,tmpbo);
+ memory->create_kokkos(k_tmpid,tmpid,nmax,MAXSPECBOND,"pair:tmpid");
+ memory->create_kokkos(k_tmpbo,tmpbo,nmax,MAXSPECBOND,"pair:tmpbo");
+ }
+
+ // FixReaxCBonds
+ d_abo = typename AT::t_ffloat_2d("reax/c/kk:abo",nmax,maxbo);
+ d_neighid = typename AT::t_tagint_2d("reax/c/kk:neighid",nmax,maxbo);
+ d_numneigh_bonds = typename AT::t_int_1d("reax/c/kk:numneigh_bonds",nmax);
}
/* ---------------------------------------------------------------------- */
@@ -3954,11 +3983,214 @@ double PairReaxCKokkos<DeviceType>::memory_usage()
bytes += nmax*17*sizeof(F_FLOAT);
bytes += maxbo*nmax*34*sizeof(F_FLOAT);
+ // FixReaxCSpecies
+ if (fixspecies_flag) {
+ bytes += MAXSPECBOND*nmax*sizeof(tagint);
+ bytes += MAXSPECBOND*nmax*sizeof(F_FLOAT);
+ }
+
+ // FixReaxCBonds
+ bytes += maxbo*nmax*sizeof(tagint);
+ bytes += maxbo*nmax*sizeof(F_FLOAT);
+ bytes += nmax*sizeof(int);
+
return bytes;
}
/* ---------------------------------------------------------------------- */
+template<class DeviceType>
+void PairReaxCKokkos<DeviceType>::FindBond(int &numbonds)
+{
+ copymode = 1;
+ Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, PairReaxFindBondZero>(0,nmax),*this);
+ DeviceType::fence();
+
+ bo_cut_bond = control->bg_cut;
+
+ atomKK->sync(execution_space,TAG_MASK);
+ tag = atomKK->k_tag.view<DeviceType>();
+
+ const int inum = list->inum;
+ NeighListKokkos<DeviceType>* k_list = static_cast<NeighListKokkos<DeviceType>*>(list);
+ d_ilist = k_list->d_ilist;
+ k_list->clean_copy();
+
+ numbonds = 0;
+ PairReaxCKokkosFindBondFunctor<DeviceType> find_bond_functor(this);
+ Kokkos::parallel_reduce(inum,find_bond_functor,numbonds);
+ DeviceType::fence();
+ copymode = 0;
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::operator()(PairReaxFindBondZero, const int &i) const {
+ d_numneigh_bonds[i] = 0;
+ for (int j = 0; j < maxbo; j++) {
+ d_neighid(i,j) = 0;
+ d_abo(i,j) = 0.0;
+ }
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::calculate_find_bond_item(int ii, int &numbonds) const
+{
+ const int i = d_ilist[ii];
+ int nj = 0;
+
+ const int j_start = d_bo_first[i];
+ const int j_end = j_start + d_bo_num[i];
+ for (int jj = j_start; jj < j_end; jj++) {
+ int j = d_bo_list[jj];
+ j &= NEIGHMASK;
+ const tagint jtag = tag[j];
+ const int j_index = jj - j_start;
+ double bo_tmp = d_BO(i,j_index);
+
+ if (bo_tmp > bo_cut_bond) {
+ d_neighid(i,nj) = jtag;
+ d_abo(i,nj) = bo_tmp;
+ nj++;
+ }
+ }
+ d_numneigh_bonds[i] = nj;
+ if (nj > numbonds) numbonds = nj;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void PairReaxCKokkos<DeviceType>::PackBondBuffer(DAT::tdual_ffloat_1d k_buf, int &nbuf_local)
+{
+ d_buf = k_buf.view<DeviceType>();
+ k_params_sing.template sync<DeviceType>();
+ atomKK->sync(execution_space,TAG_MASK|TYPE_MASK|Q_MASK|MOLECULE_MASK);
+
+ tag = atomKK->k_tag.view<DeviceType>();
+ type = atomKK->k_type.view<DeviceType>();
+ q = atomKK->k_q.view<DeviceType>();
+ if (atom->molecule)
+ molecule = atomKK->k_molecule.view<DeviceType>();
+
+ copymode = 1;
+ nlocal = atomKK->nlocal;
+ PairReaxCKokkosPackBondBufferFunctor<DeviceType> pack_bond_buffer_functor(this);
+ Kokkos::parallel_scan(nlocal,pack_bond_buffer_functor);
+ DeviceType::fence();
+ copymode = 0;
+
+ k_buf.modify<DeviceType>();
+ k_nbuf_local.modify<DeviceType>();
+
+ k_buf.sync<LMPHostType>();
+ k_nbuf_local.sync<LMPHostType>();
+ nbuf_local = k_nbuf_local.h_view();
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::pack_bond_buffer_item(int i, int &j, const bool &final) const
+{
+ if (i == 0)
+ j += 2;
+
+ if (final) {
+ d_buf[j-1] = tag[i];
+ d_buf[j+0] = type[i];
+ d_buf[j+1] = d_total_bo[i];
+ d_buf[j+2] = paramssing(type[i]).nlp_opt - d_Delta_lp[i];
+ d_buf[j+3] = q[i];
+ d_buf[j+4] = d_numneigh_bonds[i];
+ }
+ const int numbonds = d_numneigh_bonds[i];
+
+ if (final) {
+ for (int k = 5; k < 5+numbonds; k++) {
+ d_buf[j+k] = d_neighid(i,k-5);
+ }
+ }
+ j += (5+numbonds);
+
+ if (final) {
+ if (!molecule.data()) d_buf[j] = 0.0;
+ else d_buf[j] = molecule[i];
+ }
+ j++;
+
+ if (final) {
+ for (int k = 0; k < numbonds; k++) {
+ d_buf[j+k] = d_abo(i,k);
+ }
+ }
+ j += (1+numbonds);
+
+ if (final && i == nlocal-1)
+ k_nbuf_local.view<DeviceType>()() = j - 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+template<class DeviceType>
+void PairReaxCKokkos<DeviceType>::FindBondSpecies()
+{
+ copymode = 1;
+ Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, PairReaxFindBondSpeciesZero>(0,nmax),*this);
+ DeviceType::fence();
+
+ nlocal = atomKK->nlocal;
+ Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, PairReaxFindBondSpecies>(0,nlocal),*this);
+ DeviceType::fence();
+ copymode = 0;
+
+ // NOTE: Could improve performance if a Kokkos version of ComputeSpecAtom is added
+
+ k_tmpbo.modify<DeviceType>();
+ k_tmpid.modify<DeviceType>();
+ k_error_flag.modify<DeviceType>();
+
+ k_tmpbo.sync<LMPHostType>();
+ k_tmpid.sync<LMPHostType>();
+ k_error_flag.sync<LMPHostType>();
+
+ if (k_error_flag.h_view())
+ error->all(FLERR,"Increase MAXSPECBOND in reaxc_defs.h");
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::operator()(PairReaxFindBondSpeciesZero, const int &i) const {
+ for (int j = 0; j < MAXSPECBOND; j++) {
+ k_tmpbo.view<DeviceType>()(i,j) = 0.0;
+ k_tmpid.view<DeviceType>()(i,j) = 0;
+ }
+}
+
+template<class DeviceType>
+KOKKOS_INLINE_FUNCTION
+void PairReaxCKokkos<DeviceType>::operator()(PairReaxFindBondSpecies, const int &i) const {
+ int nj = 0;
+
+ const int j_start = d_bo_first[i];
+ const int j_end = j_start + d_bo_num[i];
+ for (int jj = j_start; jj < j_end; jj++) {
+ int j = d_bo_list[jj];
+ j &= NEIGHMASK;
+ if (j < i) continue;
+ const int j_index = jj - j_start;
+
+ double bo_tmp = d_BO(i,j_index);
+
+ if (bo_tmp >= 0.10 ) { // Why is this a hardcoded value?
+ k_tmpid.view<DeviceType>()(i,nj) = j;
+ k_tmpbo.view<DeviceType>()(i,nj) = bo_tmp;
+ nj++;
+ if (nj > MAXSPECBOND) k_error_flag.view<DeviceType>()() = 1;
+ }
+ }
+}
+
template class PairReaxCKokkos<LMPDeviceType>;
#ifdef KOKKOS_HAVE_CUDA
template class PairReaxCKokkos<LMPHostType>;
diff --git a/src/KOKKOS/pair_reax_c_kokkos.h b/src/KOKKOS/pair_reax_c_kokkos.h
index 8c07ee2a0355a7870c88e48aeb8522f6e4a3538d..204ad7732f9cd56d54b30b6d442f9b0f32e2c2fa 100644
--- a/src/KOKKOS/pair_reax_c_kokkos.h
+++ b/src/KOKKOS/pair_reax_c_kokkos.h
@@ -115,6 +115,13 @@ struct PairReaxComputeTorsion{};
template<int NEIGHFLAG, int EVFLAG>
struct PairReaxComputeHydrogen{};
+struct PairReaxFindBondZero{};
+
+struct PairReaxFindBondSpeciesZero{};
+
+struct PairReaxFindBondSpecies{};
+
+
template<class DeviceType>
class PairReaxCKokkos : public PairReaxC {
public:
@@ -132,6 +139,9 @@ class PairReaxCKokkos : public PairReaxC {
void *extract(const char *, int &);
void init_style();
double memory_usage();
+ void FindBond(int &);
+ void PackBondBuffer(DAT::tdual_ffloat_1d, int &);
+ void FindBondSpecies();
template<int NEIGHFLAG, int EVFLAG>
KOKKOS_INLINE_FUNCTION
@@ -245,6 +255,21 @@ class PairReaxCKokkos : public PairReaxC {
KOKKOS_INLINE_FUNCTION
void operator()(PairReaxComputeHydrogen<NEIGHFLAG,EVFLAG>, const int&) const;
+ KOKKOS_INLINE_FUNCTION
+ void operator()(PairReaxFindBondZero, const int&) const;
+
+ KOKKOS_INLINE_FUNCTION
+ void calculate_find_bond_item(int, int&) const;
+
+ KOKKOS_INLINE_FUNCTION
+ void pack_bond_buffer_item(int, int&, const bool&) const;
+
+ KOKKOS_INLINE_FUNCTION
+ void operator()(PairReaxFindBondSpeciesZero, const int&) const;
+
+ KOKKOS_INLINE_FUNCTION
+ void operator()(PairReaxFindBondSpecies, const int&) const;
+
struct params_sing{
KOKKOS_INLINE_FUNCTION
params_sing(){mass=0;chi=0;eta=0;r_s=0;r_pi=0;r_pi2=0;valency=0;valency_val=0;valency_e=0;valency_boc=0;nlp_opt=0;
@@ -359,8 +384,9 @@ class PairReaxCKokkos : public PairReaxC {
typename AT::t_x_array_randomread x;
typename AT::t_f_array f;
typename AT::t_int_1d_randomread type;
- typename AT::t_tagint_1d tag;
+ typename AT::t_tagint_1d_randomread tag;
typename AT::t_float_1d_randomread q;
+ typename AT::t_tagint_1d_randomread molecule;
DAT::tdual_efloat_1d k_eatom;
typename AT::t_efloat_1d v_eatom;
@@ -405,6 +431,7 @@ class PairReaxCKokkos : public PairReaxC {
int neighflag,newton_pair, maxnumneigh, maxhb, maxbo;
int nlocal,nall,eflag,vflag;
F_FLOAT cut_nbsq, cut_hbsq, cut_bosq, bo_cut, thb_cut, thb_cutsq;
+ F_FLOAT bo_cut_bond;
int vdwflag, lgflag;
F_FLOAT gp[39], p_boc1, p_boc2;
@@ -418,6 +445,49 @@ class PairReaxCKokkos : public PairReaxC {
tdual_LR_lookup_table_kk_2d k_LR;
t_LR_lookup_table_kk_2d d_LR;
+
+ DAT::tdual_int_2d k_tmpid;
+ DAT::tdual_ffloat_2d k_tmpbo;
+ DAT::tdual_int_scalar k_error_flag;
+
+ typename AT::t_int_1d d_numneigh_bonds;
+ typename AT::t_tagint_2d d_neighid;
+ typename AT::t_ffloat_2d d_abo;
+
+ typename AT::t_ffloat_1d d_buf;
+ DAT::tdual_int_scalar k_nbuf_local;
+};
+
+template <class DeviceType>
+struct PairReaxCKokkosFindBondFunctor {
+ typedef DeviceType device_type;
+ typedef int value_type;
+ PairReaxCKokkos<DeviceType> c;
+ PairReaxCKokkosFindBondFunctor(PairReaxCKokkos<DeviceType>* c_ptr):c(*c_ptr) {};
+
+ KOKKOS_INLINE_FUNCTION
+ void join(volatile int &dst,
+ const volatile int &src) const {
+ dst = MAX(dst,src);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ void operator()(const int ii, int &numbonds) const {
+ c.calculate_find_bond_item(ii,numbonds);
+ }
+};
+
+template <class DeviceType>
+struct PairReaxCKokkosPackBondBufferFunctor {
+ typedef DeviceType device_type;
+ typedef int value_type;
+ PairReaxCKokkos<DeviceType> c;
+ PairReaxCKokkosPackBondBufferFunctor(PairReaxCKokkos<DeviceType>* c_ptr):c(*c_ptr) {};
+
+ KOKKOS_INLINE_FUNCTION
+ void operator()(const int ii, int &j, const bool &final) const {
+ c.pack_bond_buffer_item(ii,j,final);
+ }
};
}
diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp
index 3a9f7cc4d02a85b867aa69ff818fef5f186110a9..858bc83a3665fb0b324221b80b99034e91430efa 100644
--- a/src/MOLECULE/bond_table.cpp
+++ b/src/MOLECULE/bond_table.cpp
@@ -585,17 +585,26 @@ double BondTable::splint(double *xa, double *ya, double *y2a, int n, double x)
/* ----------------------------------------------------------------------
calculate potential u and force f at distance x
- insure x is between bond min/max
+ insure x is between bond min/max, exit with error if not
------------------------------------------------------------------------- */
void BondTable::uf_lookup(int type, double x, double &u, double &f)
{
int itable;
double fraction,a,b;
+ char estr[128];
Table *tb = &tables[tabindex[type]];
- x = MAX(x,tb->lo);
- x = MIN(x,tb->hi);
+ if (x < tb->lo) {
+ sprintf(estr,"Bond length < table inner cutoff: "
+ "type %d length %g",type,x);
+ error->one(FLERR,estr);
+ }
+ if (x > tb->hi) {
+ sprintf(estr,"Bond length > table outer cutoff: "
+ "type %d length %g",type,x);
+ error->one(FLERR,estr);
+ }
if (tabstyle == LINEAR) {
itable = static_cast<int> ((x - tb->lo) * tb->invdelta);
diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp
index dc97db194253e4a2ed8d6367da0b839ecdc51e31..350009f6bf9488b3f225caebe9064821343e1129 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.cpp
+++ b/src/USER-COLVARS/colvarproxy_lammps.cpp
@@ -333,7 +333,10 @@ int colvarproxy_lammps::backup_file(char const *filename)
int colvarproxy_lammps::smp_enabled()
{
- return COLVARS_OK;
+ if (b_smp_active) {
+ return COLVARS_OK;
+ }
+ return COLVARS_ERROR;
}
diff --git a/src/USER-REAXC/compute_spec_atom.cpp b/src/USER-REAXC/compute_spec_atom.cpp
index 38cf9ad6e6524f495b177f60749b66c0110b2906..4af8efcae71b324da7d5b8a872f13e5e491dae59 100644
--- a/src/USER-REAXC/compute_spec_atom.cpp
+++ b/src/USER-REAXC/compute_spec_atom.cpp
@@ -44,6 +44,8 @@ ComputeSpecAtom::ComputeSpecAtom(LAMMPS *lmp, int narg, char **arg) :
// Initiate reaxc
reaxc = (PairReaxC *) force->pair_match("reax/c",1);
+ if (reaxc == NULL)
+ reaxc = (PairReaxC *) force->pair_match("reax/c/kk",1);
pack_choice = new FnPtrPack[nvalues];
diff --git a/src/USER-REAXC/fix_reaxc_bonds.cpp b/src/USER-REAXC/fix_reaxc_bonds.cpp
index 5f653557a5dabdd3a03f3330c0bf3533f2a8ab54..543669de766cce7c3a1466988623af1f564b2dd0 100644
--- a/src/USER-REAXC/fix_reaxc_bonds.cpp
+++ b/src/USER-REAXC/fix_reaxc_bonds.cpp
@@ -107,11 +107,10 @@ void FixReaxCBonds::setup(int vflag)
void FixReaxCBonds::init()
{
- Pair *pair_kk = force->pair_match("reax/c/kk",1);
- if (pair_kk != NULL) error->all(FLERR,"Cannot (yet) use fix reax/c/bonds with "
- "pair_style reax/c/kk");
-
reaxc = (PairReaxC *) force->pair_match("reax/c",1);
+ if (reaxc == NULL)
+ reaxc = (PairReaxC *) force->pair_match("reax/c/kk",1);
+
if (reaxc == NULL) error->all(FLERR,"Cannot use fix reax/c/bonds without "
"pair_style reax/c");
diff --git a/src/USER-REAXC/fix_reaxc_bonds.h b/src/USER-REAXC/fix_reaxc_bonds.h
index b91085163b51c1f0eb73811c447668e3dd52175c..e1dcd82c0e43ed96323668a918dfd9259e05939e 100644
--- a/src/USER-REAXC/fix_reaxc_bonds.h
+++ b/src/USER-REAXC/fix_reaxc_bonds.h
@@ -29,13 +29,13 @@ namespace LAMMPS_NS {
class FixReaxCBonds : public Fix {
public:
FixReaxCBonds(class LAMMPS *, int, char **);
- ~FixReaxCBonds();
+ virtual ~FixReaxCBonds();
int setmask();
- void init();
+ virtual void init();
void setup(int);
void end_of_step();
- private:
+ protected:
int me, nprocs, nmax, ntypes, maxsize;
int *numneigh;
tagint **neighid;
@@ -44,12 +44,12 @@ class FixReaxCBonds : public Fix {
void allocate();
void destroy();
- void Output_ReaxC_Bonds(bigint, FILE *);
+ virtual void Output_ReaxC_Bonds(bigint, FILE *);
void FindBond(struct _reax_list*, int &);
void PassBuffer(double *, int &);
void RecvBuffer(double *, int, int, int, int);
int nint(const double &);
- double memory_usage();
+ virtual double memory_usage();
bigint nvalid, nextvalid();
struct _reax_list *lists;
diff --git a/src/USER-REAXC/fix_reaxc_species.cpp b/src/USER-REAXC/fix_reaxc_species.cpp
index c8acc5f0e0a803fbab1fc347686ba2f5ea36cd37..ead73f02a1328cef35a4af076ee8e1f0818aa490 100644
--- a/src/USER-REAXC/fix_reaxc_species.cpp
+++ b/src/USER-REAXC/fix_reaxc_species.cpp
@@ -284,11 +284,10 @@ void FixReaxCSpecies::init()
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use fix reax/c/species unless atoms have IDs");
- Pair *pair_kk = force->pair_match("reax/c/kk",1);
- if (pair_kk != NULL) error->all(FLERR,"Cannot (yet) use fix reax/c/species with "
- "pair_style reax/c/kk");
-
reaxc = (PairReaxC *) force->pair_match("reax/c",1);
+ if (reaxc == NULL)
+ reaxc = (PairReaxC *) force->pair_match("reax/c/kk",1);
+
if (reaxc == NULL) error->all(FLERR,"Cannot use fix reax/c/species without "
"pair_style reax/c");
@@ -485,7 +484,7 @@ void FixReaxCSpecies::Output_ReaxC_Bonds(bigint ntimestep, FILE *fp)
/* ---------------------------------------------------------------------- */
-AtomCoord chAnchor(AtomCoord in1, AtomCoord in2)
+AtomCoord FixReaxCSpecies::chAnchor(AtomCoord in1, AtomCoord in2)
{
if (in1.x < in2.x)
return in1;
diff --git a/src/USER-REAXC/fix_reaxc_species.h b/src/USER-REAXC/fix_reaxc_species.h
index a3a47b29ddfa6b7b42d7377372eb8f1cd2db2354..872ea2528f06df00cd27ee0e76fd19fa66fdad51 100644
--- a/src/USER-REAXC/fix_reaxc_species.h
+++ b/src/USER-REAXC/fix_reaxc_species.h
@@ -38,15 +38,15 @@ typedef struct {
class FixReaxCSpecies : public Fix {
public:
FixReaxCSpecies(class LAMMPS *, int, char **);
- ~FixReaxCSpecies();
+ virtual ~FixReaxCSpecies();
int setmask();
- void init();
+ virtual void init();
void init_list(int, class NeighList *);
void setup(int);
void post_integrate();
double compute_vector(int);
- private:
+ protected:
int me, nprocs, nmax, nlocal, ntypes, ntotal;
int nrepeat, nfreq, posfreq;
int Nmoltype, vector_nmole, vector_nspec;
@@ -67,7 +67,8 @@ class FixReaxCSpecies : public Fix {
void Output_ReaxC_Bonds(bigint, FILE *);
void create_compute();
void create_fix();
- void FindMolecule();
+ AtomCoord chAnchor(AtomCoord, AtomCoord);
+ virtual void FindMolecule();
void SortMolecule(int &);
void FindSpecies(int, int &);
void WriteFormulas(int, int);
diff --git a/src/atom.cpp b/src/atom.cpp
index 8e48611284b92700942c4ffcb38955af2a0252fe..e1ec60fd99fc1dcc229e48f1e68a155437cb4402 100644
--- a/src/atom.cpp
+++ b/src/atom.cpp
@@ -26,11 +26,14 @@
#include "force.h"
#include "modify.h"
#include "fix.h"
+#include "compute.h"
#include "output.h"
#include "thermo.h"
#include "update.h"
#include "domain.h"
#include "group.h"
+#include "input.h"
+#include "variable.h"
#include "molecule.h"
#include "atom_masks.h"
#include "math_const.h"
@@ -1388,6 +1391,33 @@ void Atom::data_bodies(int n, char *buf, AtomVecBody *avec_body,
delete [] dvalues;
}
+/* ----------------------------------------------------------------------
+ init per-atom fix/compute/variable values for newly created atoms
+ called from create_atoms, read_data, read_dump,
+ lib::lammps_create_atoms()
+ fixes, computes, variables may or may not exist when called
+------------------------------------------------------------------------- */
+
+void Atom::data_fix_compute_variable(int nprev, int nnew)
+{
+ for (int m = 0; m < modify->nfix; m++) {
+ Fix *fix = modify->fix[m];
+ if (fix->create_attribute)
+ for (int i = nprev; i < nnew; i++)
+ fix->set_arrays(i);
+ }
+
+ for (int m = 0; m < modify->ncompute; m++) {
+ Compute *compute = modify->compute[m];
+ if (compute->create_attribute)
+ for (int i = nprev; i < nnew; i++)
+ compute->set_arrays(i);
+ }
+
+ for (int i = nprev; i < nnew; i++)
+ input->variable->set_arrays(i);
+}
+
/* ----------------------------------------------------------------------
allocate arrays of length ntypes
only done after ntypes is set
diff --git a/src/atom.h b/src/atom.h
index 61cd9673bfc838e730d30698250780411fe7b15a..748d43cddd1c4e685986c7c9e0c048e7d9ebff1d 100644
--- a/src/atom.h
+++ b/src/atom.h
@@ -225,15 +225,14 @@ class Atom : protected Pointers {
void data_atoms(int, char *, tagint, int, int, double *);
void data_vels(int, char *, tagint);
-
void data_bonds(int, char *, int *, tagint, int);
void data_angles(int, char *, int *, tagint, int);
void data_dihedrals(int, char *, int *, tagint, int);
void data_impropers(int, char *, int *, tagint, int);
-
void data_bonus(int, char *, class AtomVec *, tagint);
void data_bodies(int, char *, class AtomVecBody *, tagint);
-
+ void data_fix_compute_variable(int, int);
+
virtual void allocate_type_arrays();
void set_mass(const char *, int);
void set_mass(int, double);
diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp
index 82315e7dc5f813374fdeace3d79c53531d60c4dc..776a94c8b933d58960cecf5ac6ba1178f780ba75 100644
--- a/src/create_atoms.cpp
+++ b/src/create_atoms.cpp
@@ -359,30 +359,9 @@ void CreateAtoms::command(int narg, char **arg)
else if (style == RANDOM) add_random();
else add_lattice();
- // invoke set_arrays() for fixes/computes/variables
- // that need initialization of attributes of new atoms
- // don't use modify->create_attributes() since would be inefficient
- // for large number of atoms
- // note that for typical early use of create_atoms,
- // no fixes/computes/variables exist yet
-
- int nlocal = atom->nlocal;
- for (int m = 0; m < modify->nfix; m++) {
- Fix *fix = modify->fix[m];
- if (fix->create_attribute)
- for (int i = nlocal_previous; i < nlocal; i++)
- fix->set_arrays(i);
- }
-
- for (int m = 0; m < modify->ncompute; m++) {
- Compute *compute = modify->compute[m];
- if (compute->create_attribute)
- for (int i = nlocal_previous; i < nlocal; i++)
- compute->set_arrays(i);
- }
-
- for (int i = nlocal_previous; i < nlocal; i++)
- input->variable->set_arrays(i);
+ // init per-atom fix/compute/variable values for created atoms
+
+ atom->data_fix_compute_variable(nlocal_previous,atom->nlocal);
// set new total # of atoms and error check
diff --git a/src/create_box.cpp b/src/create_box.cpp
index c7fd98b2b2898cb3b208c3dc7353f7408a82831c..3a131099813e086d1dd961a28748824959cfe50e 100644
--- a/src/create_box.cpp
+++ b/src/create_box.cpp
@@ -159,6 +159,7 @@ void CreateBox::command(int narg, char **arg)
} else if (strcmp(arg[iarg],"extra/special/per/atom") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal create_box command");
force->special_extra = force->inumeric(FLERR,arg[iarg+1]);
+ atom->maxspecial += force->special_extra;
iarg += 2;
} else error->all(FLERR,"Illegal create_box command");
}
diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp
index 4ac7c0262203c89db1f7f9208a40606622d92722..116c2a2f1e057431644c82b2747d7d95eac6095c 100644
--- a/src/delete_atoms.cpp
+++ b/src/delete_atoms.cpp
@@ -59,7 +59,9 @@ void DeleteAtoms::command(int narg, char **arg)
bigint ndihedrals_previous = atom->ndihedrals;
bigint nimpropers_previous = atom->nimpropers;
- // delete the atoms
+ // flag atoms for deletion
+
+ allflag = 0;
if (strcmp(arg[0],"group") == 0) delete_group(narg,arg);
else if (strcmp(arg[0],"region") == 0) delete_region(narg,arg);
@@ -67,36 +69,44 @@ void DeleteAtoms::command(int narg, char **arg)
else if (strcmp(arg[0],"porosity") == 0) delete_porosity(narg,arg);
else error->all(FLERR,"Illegal delete_atoms command");
- // optionally delete additional bonds or atoms in molecules
+ // if allflag = 1, just reset atom->nlocal
+ // else delete atoms one by one
- if (bond_flag) delete_bond();
- if (mol_flag) delete_molecule();
+ if (allflag) atom->nlocal = 0;
+ else {
- // delete local atoms flagged in dlist
- // reset nlocal
+ // optionally delete additional bonds or atoms in molecules
- AtomVec *avec = atom->avec;
- int nlocal = atom->nlocal;
+ if (bond_flag) delete_bond();
+ if (mol_flag) delete_molecule();
- int i = 0;
- while (i < nlocal) {
- if (dlist[i]) {
- avec->copy(nlocal-1,i,1);
- dlist[i] = dlist[nlocal-1];
- nlocal--;
- } else i++;
- }
+ // delete local atoms flagged in dlist
+ // reset nlocal
- atom->nlocal = nlocal;
- memory->destroy(dlist);
+ AtomVec *avec = atom->avec;
+ int nlocal = atom->nlocal;
+ int i = 0;
+ while (i < nlocal) {
+ if (dlist[i]) {
+ avec->copy(nlocal-1,i,1);
+ dlist[i] = dlist[nlocal-1];
+ nlocal--;
+ } else i++;
+ }
+
+ atom->nlocal = nlocal;
+ memory->destroy(dlist);
+ }
+
// if non-molecular system and compress flag set,
// reset atom tags to be contiguous
// set all atom IDs to 0, call tag_extend()
if (atom->molecular == 0 && compress_flag) {
tagint *tag = atom->tag;
- for (i = 0; i < nlocal; i++) tag[i] = 0;
+ int nlocal = atom->nlocal;
+ for (int i = 0; i < nlocal; i++) tag[i] = 0;
atom->tag_extend();
}
@@ -185,6 +195,13 @@ void DeleteAtoms::delete_group(int narg, char **arg)
if (igroup == -1) error->all(FLERR,"Could not find delete_atoms group ID");
options(narg-2,&arg[2]);
+ // check for special case of group = all
+
+ if (strcmp(arg[1],"all") == 0) {
+ allflag = 1;
+ return;
+ }
+
// allocate and initialize deletion list
int nlocal = atom->nlocal;
diff --git a/src/delete_atoms.h b/src/delete_atoms.h
index 6fc2ea1f90b8134c1322f37e463f2aa7fa6cdf3e..62ba47d715d47efab94161362d21c5dcfa542847 100644
--- a/src/delete_atoms.h
+++ b/src/delete_atoms.h
@@ -32,7 +32,7 @@ class DeleteAtoms : protected Pointers {
private:
int *dlist;
- int compress_flag,bond_flag,mol_flag;
+ int allflag,compress_flag,bond_flag,mol_flag;
std::map<tagint,int> *hash;
void delete_group(int, char **);
diff --git a/src/domain.cpp b/src/domain.cpp
index 2c3f61401d9112a3504afa0e504e2db9edb3fecc..c348167dfb623e4925d858db7b649ff6de491d8d 100644
--- a/src/domain.cpp
+++ b/src/domain.cpp
@@ -1496,6 +1496,28 @@ void Domain::image_flip(int m, int n, int p)
}
}
+/* ----------------------------------------------------------------------
+ return 1 if this proc owns atom with coords x, else return 0
+ x is returned remapped into periodic box
+------------------------------------------------------------------------- */
+
+int Domain::ownatom(double *x)
+{
+ double lamda[3];
+ double *coord;
+
+ remap(x);
+ if (triclinic) {
+ x2lamda(x,lamda);
+ coord = lamda;
+ } else coord = x;
+
+ if (coord[0] >= sublo[0] && coord[0] < subhi[0] &&
+ coord[1] >= sublo[1] && coord[1] < subhi[1] &&
+ coord[2] >= sublo[2] && coord[2] < subhi[2]) return 1;
+ return 0;
+}
+
/* ----------------------------------------------------------------------
create a lattice
------------------------------------------------------------------------- */
@@ -1929,5 +1951,3 @@ void Domain::lamda_box_corners(double *lo, double *hi)
corners[7][0] = hi[0]; corners[7][1] = hi[1]; corners[7][2] = subhi_lamda[2];
lamda2x(corners[7],corners[7]);
}
-
-
diff --git a/src/domain.h b/src/domain.h
index ef3de4207177438854ba13281334442d3bf8a2d4..e2425f458576da20197b8e9675476d33aaabcdaf 100644
--- a/src/domain.h
+++ b/src/domain.h
@@ -121,6 +121,8 @@ class Domain : protected Pointers {
void unmap(double *, imageint);
void unmap(double *, imageint, double *);
void image_flip(int, int, int);
+ int ownatom(double *);
+
void set_lattice(int, char **);
void add_region(int, char **);
void delete_region(int, char **);
diff --git a/src/fix_halt.h b/src/fix_halt.h
index bc7f0055a5b5a64d011188c5294cd6181250a7dd..156397a4b5a76cb979e3b524bbf544edbcf446db 100644
--- a/src/fix_halt.h
+++ b/src/fix_halt.h
@@ -36,7 +36,7 @@ class FixHalt : public Fix {
private:
int attribute,operation,eflag,ivar;
- double attvalue,value;
+ double value;
char *idvar;
double bondmax();
diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp
index 01ba81e5547c09c2dfc84418f4824e90206c3b90..1799a00558c039f3484a39b3eb70ad03aca7bd0d 100644
--- a/src/fix_nh.cpp
+++ b/src/fix_nh.cpp
@@ -51,11 +51,10 @@ enum{ISO,ANISO,TRICLINIC};
---------------------------------------------------------------------- */
FixNH::FixNH(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
-tstat_flag(0), pstat_flag(0),
rfix(NULL), id_dilate(NULL), irregular(NULL), id_temp(NULL), id_press(NULL),
-tcomputeflag(0), pcomputeflag(0), eta(NULL), eta_dot(NULL), eta_dotdot(NULL),
+eta(NULL), eta_dot(NULL), eta_dotdot(NULL),
eta_mass(NULL), etap(NULL), etap_dot(NULL), etap_dotdot(NULL),
-etap_mass(NULL), mpchain(0)
+etap_mass(NULL)
{
if (narg < 4) error->all(FLERR,"Illegal fix nvt/npt/nph command");
@@ -478,8 +477,10 @@ etap_mass(NULL), mpchain(0)
// pre_exchange only required if flips can occur due to shape changes
- if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5])) pre_exchange_flag = 1;
- if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 || domain->xy != 0.0))
+ if (flipflag && (p_flag[3] || p_flag[4] || p_flag[5]))
+ pre_exchange_flag = 1;
+ if (flipflag && (domain->yz != 0.0 || domain->xz != 0.0 ||
+ domain->xy != 0.0))
pre_exchange_flag = 1;
}
diff --git a/src/input.cpp b/src/input.cpp
index 905e7f0590438ad068f347b4b30cd87d5a2102b0..fdb551d3da455f987daf03813377986c699e26e3 100644
--- a/src/input.cpp
+++ b/src/input.cpp
@@ -278,7 +278,8 @@ void Input::file(const char *filename)
}
/* ----------------------------------------------------------------------
- copy command in single to line, parse and execute it
+ invoke one command in single
+ first copy to line, then parse, then execute it
return command name to caller
------------------------------------------------------------------------- */
diff --git a/src/library.cpp b/src/library.cpp
index 44e2c407825b51cc0b523173475a23eb302dcbb8..90c7cfa9d5162d120131385b3bed0cdd858fcebe 100644
--- a/src/library.cpp
+++ b/src/library.cpp
@@ -18,9 +18,11 @@
#include <string.h>
#include <stdlib.h>
#include "library.h"
+#include "lmptype.h"
#include "lammps.h"
#include "universe.h"
#include "input.h"
+#include "atom_vec.h"
#include "atom.h"
#include "domain.h"
#include "update.h"
@@ -38,8 +40,12 @@
using namespace LAMMPS_NS;
+// ----------------------------------------------------------------------
+// utility macros
+// ----------------------------------------------------------------------
+
/* ----------------------------------------------------------------------
- Utility macros for optional code path which captures all exceptions
+ macros for optional code path which captures all exceptions
and stores the last error message. These assume there is a variable lmp
which is a pointer to the current LAMMPS instance.
@@ -76,6 +82,37 @@ using namespace LAMMPS_NS;
#define END_CAPTURE
#endif
+// ----------------------------------------------------------------------
+// helper functions, not in library API
+// ----------------------------------------------------------------------
+
+/* ----------------------------------------------------------------------
+ concatenate one or more LAMMPS input lines starting at ptr
+ removes NULL terminator when last printable char of line = '&'
+ by replacing both NULL and '&' with space character
+ repeat as many times as needed
+ on return, ptr now points to longer line
+------------------------------------------------------------------------- */
+
+void concatenate_lines(char *ptr)
+{
+ int nend = strlen(ptr);
+ int n = nend-1;
+ while (n && isspace(ptr[n])) n--;
+ while (ptr[n] == '&') {
+ ptr[nend] = ' ';
+ ptr[n] = ' ';
+ strtok(ptr,"\n");
+ nend = strlen(ptr);
+ n = nend-1;
+ while (n && isspace(ptr[n])) n--;
+ }
+}
+
+// ----------------------------------------------------------------------
+// library API functions to create/destroy an instance of LAMMPS
+// and communicate commands to it
+// ----------------------------------------------------------------------
/* ----------------------------------------------------------------------
create an instance of LAMMPS and return pointer to it
@@ -172,12 +209,14 @@ void lammps_file(void *ptr, char *str)
/* ----------------------------------------------------------------------
process a single input command in str
+ does not matter if str ends in newline
+ return command name to caller
------------------------------------------------------------------------- */
char *lammps_command(void *ptr, char *str)
{
LAMMPS *lmp = (LAMMPS *) ptr;
- char * result = NULL;
+ char *result = NULL;
BEGIN_CAPTURE
{
@@ -188,6 +227,69 @@ char *lammps_command(void *ptr, char *str)
return result;
}
+/* ----------------------------------------------------------------------
+ process multiple input commands in cmds = list of strings
+ does not matter if each string ends in newline
+ create long contatentated string for processing by commands_string()
+ insert newlines in concatenated string as needed
+------------------------------------------------------------------------- */
+
+void lammps_commands_list(void *ptr, int ncmd, char **cmds)
+{
+ LAMMPS *lmp = (LAMMPS *) ptr;
+
+ int n = ncmd+1;
+ for (int i = 0; i < ncmd; i++) n += strlen(cmds[i]);
+
+ char *str = (char *) lmp->memory->smalloc(n,"lib/commands/list:str");
+ str[0] = '\0';
+ n = 0;
+
+ for (int i = 0; i < ncmd; i++) {
+ strcpy(&str[n],cmds[i]);
+ n += strlen(cmds[i]);
+ if (str[n-1] != '\n') {
+ str[n] = '\n';
+ str[n+1] = '\0';
+ n++;
+ }
+ }
+
+ lammps_commands_string(ptr,str);
+ lmp->memory->sfree(str);
+}
+
+/* ----------------------------------------------------------------------
+ process multiple input commands in single long str, separated by newlines
+ single command can span multiple lines via continuation characters
+ multi-line commands enabled by triple quotes will not work
+------------------------------------------------------------------------- */
+
+void lammps_commands_string(void *ptr, char *str)
+{
+ LAMMPS *lmp = (LAMMPS *) ptr;
+
+ // make copy of str so can strtok() it
+
+ int n = strlen(str) + 1;
+ char *copy = new char[n];
+ strcpy(copy,str);
+
+ BEGIN_CAPTURE
+ {
+ char *ptr = strtok(copy,"\n");
+ if (ptr) concatenate_lines(ptr);
+ while (ptr) {
+ lmp->input->one(ptr);
+ ptr = strtok(NULL,"\n");
+ if (ptr) concatenate_lines(ptr);
+ }
+ }
+ END_CAPTURE
+
+ delete [] copy;
+}
+
/* ----------------------------------------------------------------------
clean-up function to free memory allocated by lib and returned to caller
------------------------------------------------------------------------- */
@@ -197,6 +299,10 @@ void lammps_free(void *ptr)
free(ptr);
}
+// ----------------------------------------------------------------------
+// library API functions to extract info from LAMMPS or set info in LAMMPS
+// ----------------------------------------------------------------------
+
/* ----------------------------------------------------------------------
add LAMMPS-specific library functions
all must receive LAMMPS pointer as argument
@@ -209,6 +315,8 @@ void lammps_free(void *ptr)
returns a void pointer to the entity
which the caller can cast to the proper data type
returns a NULL if name not listed below
+ this function need only be invoked once
+ the returned pointer is a permanent valid reference to the quantity
customize by adding names
------------------------------------------------------------------------- */
@@ -232,14 +340,16 @@ void *lammps_extract_global(void *ptr, char *name)
if (strcmp(name,"ndihedrals") == 0) return (void *) &lmp->atom->ndihedrals;
if (strcmp(name,"nimpropers") == 0) return (void *) &lmp->atom->nimpropers;
if (strcmp(name,"nlocal") == 0) return (void *) &lmp->atom->nlocal;
+ if (strcmp(name,"nghost") == 0) return (void *) &lmp->atom->nghost;
+ if (strcmp(name,"nmax") == 0) return (void *) &lmp->atom->nmax;
if (strcmp(name,"ntimestep") == 0) return (void *) &lmp->update->ntimestep;
- // NOTE: we cannot give access to the thermo "time" data by reference,
- // as that is a recomputed property. Only "atime" can be provided as pointer.
- // please use lammps_get_thermo() defined below to access all supported
- // thermo keywords by value.
+ // update atime can be referenced as a pointer
+ // thermo "timer" data cannot be, since it is computed on request
+ // lammps_get_thermo() can access all thermo keywords by value
if (strcmp(name,"atime") == 0) return (void *) &lmp->update->atime;
+ if (strcmp(name,"atimestep") == 0) return (void *) &lmp->update->atimestep;
return NULL;
}
@@ -250,6 +360,8 @@ void *lammps_extract_global(void *ptr, char *name)
returns a void pointer to the entity
which the caller can cast to the proper data type
returns a NULL if Atom::extract() does not recognize the name
+ the returned pointer is not a permanent valid reference to the
+ per-atom quantity, since LAMMPS may reallocate per-atom data
customize by adding names to Atom::extract()
------------------------------------------------------------------------- */
@@ -261,6 +373,7 @@ void *lammps_extract_atom(void *ptr, char *name)
/* ----------------------------------------------------------------------
extract a pointer to an internal LAMMPS compute-based entity
+ the compute is invoked if its value(s) is not current
id = compute ID
style = 0 for global data, 1 for per-atom data, 2 for local data
type = 0 for scalar, 1 for vector, 2 for array
@@ -275,6 +388,8 @@ void *lammps_extract_atom(void *ptr, char *name)
returns a void pointer to the compute's internal data structure
for the entity which the caller can cast to the proper data type
returns a NULL if id is not recognized or style/type not supported
+ the returned pointer is not a permanent valid reference to the
+ compute data, this function should be re-invoked
IMPORTANT: if the compute is not current it will be invoked
LAMMPS cannot easily check here if it is valid to invoke the compute,
so caller must insure that it is OK
@@ -492,11 +607,11 @@ int lammps_set_variable(void *ptr, char *name, char *str)
}
/* ----------------------------------------------------------------------
- return the current value of a thermo keyword as double.
+ return the current value of a thermo keyword as a double
unlike lammps_extract_global() this does not give access to the
- storage of the data in question, and thus needs to be called
- again to retrieve an updated value. The upshot is that it allows
- accessing information that is only computed on-the-fly.
+ storage of the data in question
+ instead it triggers the Thermo class to compute the current value
+ and returns it
------------------------------------------------------------------------- */
double lammps_get_thermo(void *ptr, char *name)
@@ -529,12 +644,13 @@ int lammps_get_natoms(void *ptr)
/* ----------------------------------------------------------------------
gather the named atom-based entity across all processors
+ atom IDs must be consecutive from 1 to N
name = desired quantity, e.g. x or charge
type = 0 for integer values, 1 for double values
count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
return atom-based values in data, ordered by count, then by atom ID
e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
- data must be pre-allocated by caller to correct length
+ data must be pre-allocated by caller to correct length
------------------------------------------------------------------------- */
void lammps_gather_atoms(void *ptr, char *name,
@@ -547,7 +663,8 @@ void lammps_gather_atoms(void *ptr, char *name,
// error if tags are not defined or not consecutive
int flag = 0;
- if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;
+ if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0)
+ flag = 1;
if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
if (flag) {
if (lmp->comm->me == 0)
@@ -586,7 +703,7 @@ void lammps_gather_atoms(void *ptr, char *name,
for (j = 0; j < count; j++)
copy[offset++] = array[i][0];
}
-
+
MPI_Allreduce(copy,data,count*natoms,MPI_INT,MPI_SUM,lmp->world);
lmp->memory->destroy(copy);
@@ -623,6 +740,7 @@ void lammps_gather_atoms(void *ptr, char *name,
/* ----------------------------------------------------------------------
scatter the named atom-based entity across all processors
+ atom IDs must be consecutive from 1 to N
name = desired quantity, e.g. x or charge
type = 0 for integer values, 1 for double values
count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f
@@ -640,7 +758,8 @@ void lammps_scatter_atoms(void *ptr, char *name,
// error if tags are not defined or not consecutive or no atom map
int flag = 0;
- if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0) flag = 1;
+ if (lmp->atom->tag_enable == 0 || lmp->atom->tag_consecutive() == 0)
+ flag = 1;
if (lmp->atom->natoms > MAXSMALLINT) flag = 1;
if (lmp->atom->map_style == 0) flag = 1;
if (flag) {
@@ -702,9 +821,91 @@ void lammps_scatter_atoms(void *ptr, char *name,
END_CAPTURE
}
+/* ----------------------------------------------------------------------
+ create N atoms and assign them to procs based on coords
+ id = atom IDs (optional, NULL if just use 1 to N)
+ type = N-length vector of atom types (required)
+ x = 3N-length vector of atom coords (required)
+ v = 3N-length vector of atom velocities (optional, NULL if just 0.0)
+ x,v = ordered by xyz, then by atom
+ e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],...
+------------------------------------------------------------------------- */
+
+void lammps_create_atoms(void *ptr, int n, tagint *id, int *type,
+ double *x, double *v)
+{
+ LAMMPS *lmp = (LAMMPS *) ptr;
+
+ BEGIN_CAPTURE
+ {
+ // error if box does not exist or tags not defined
+
+ int flag = 0;
+ if (lmp->domain->box_exist == 0) flag = 1;
+ if (lmp->atom->tag_enable == 0) flag = 1;
+ if (flag) {
+ if (lmp->comm->me == 0)
+ lmp->error->warning(FLERR,"Library error in lammps_create_atoms");
+ return;
+ }
+
+ // loop over N atoms of entire system
+ // if this proc owns it based on coords, invoke create_atom()
+ // optionally set atom tags and velocities
+
+ Atom *atom = lmp->atom;
+ Domain *domain = lmp->domain;
+ int nlocal = atom->nlocal;
+
+ int nprev = nlocal;
+ double xdata[3];
+
+ for (int i = 0; i < n; i++) {
+ xdata[0] = x[3*i];
+ xdata[1] = x[3*i+1];
+ xdata[2] = x[3*i+2];
+ if (!domain->ownatom(xdata)) continue;
+
+ atom->avec->create_atom(type[i],xdata);
+ if (id) atom->tag[nlocal] = id[i];
+ else atom->tag[nlocal] = i+1;
+ if (v) {
+ atom->v[nlocal][0] = v[3*i];
+ atom->v[nlocal][1] = v[3*i+1];
+ atom->v[nlocal][2] = v[3*i+2];
+ }
+ nlocal++;
+ }
+
+ // need to reset atom->natoms inside LAMMPS
+
+ bigint ncurrent = nlocal;
+ MPI_Allreduce(&ncurrent,&lmp->atom->natoms,1,MPI_LMP_BIGINT,
+ MPI_SUM,lmp->world);
+
+ // init per-atom fix/compute/variable values for created atoms
+
+ atom->data_fix_compute_variable(nprev,nlocal);
+
+ // if global map exists, reset it
+ // invoke map_init() b/c atom count has grown
+
+ if (lmp->atom->map_style) {
+ lmp->atom->map_init();
+ lmp->atom->map_set();
+ }
+ }
+ END_CAPTURE
+}
+
+// ----------------------------------------------------------------------
+// library API functions for error handling
+// ----------------------------------------------------------------------
+
#ifdef LAMMPS_EXCEPTIONS
+
/* ----------------------------------------------------------------------
- Check if a new error message
+ check if a new error message
------------------------------------------------------------------------- */
int lammps_has_error(void *ptr) {
@@ -714,8 +915,8 @@ int lammps_has_error(void *ptr) {
}
/* ----------------------------------------------------------------------
- Copy the last error message of LAMMPS into a character buffer
- The return value encodes which type of error it is.
+ copy the last error message of LAMMPS into a character buffer
+ return value encodes which type of error:
1 = normal error (recoverable)
2 = abort error (non-recoverable)
------------------------------------------------------------------------- */
@@ -732,4 +933,5 @@ int lammps_get_last_error_message(void *ptr, char * buffer, int buffer_size) {
}
return 0;
}
+
#endif
diff --git a/src/library.h b/src/library.h
index e4dffc63d1ccfe2cc2d456cc0ca207ea78397af5..a05b7b4eb8a79b1dd599e76fa9807c2f84207535 100644
--- a/src/library.h
+++ b/src/library.h
@@ -30,6 +30,8 @@ void lammps_close(void *);
int lammps_version(void *);
void lammps_file(void *, char *);
char *lammps_command(void *, char *);
+void lammps_commands_list(void *, int, char **);
+void lammps_commands_string(void *, char *);
void lammps_free(void *);
void *lammps_extract_global(void *, char *);
@@ -40,10 +42,11 @@ void *lammps_extract_variable(void *, char *, char *);
int lammps_set_variable(void *, char *, char *);
double lammps_get_thermo(void *, char *);
-int lammps_get_natoms(void *);
+int lammps_get_natoms(void *);
void lammps_gather_atoms(void *, char *, int, int, void *);
void lammps_scatter_atoms(void *, char *, int, int, void *);
+void lammps_create_atoms(void *, int, int *, int *, double *, double *);
#ifdef LAMMPS_EXCEPTIONS
int lammps_has_error(void *);
diff --git a/src/molecule.cpp b/src/molecule.cpp
index 5e9aa292d0072bca21f8bbccc36e29b22a520217..0febb784c84566c0298c00d78d3d2bca2e5ccc08 100644
--- a/src/molecule.cpp
+++ b/src/molecule.cpp
@@ -595,6 +595,10 @@ void Molecule::read(int flag)
// set maxspecial on first pass, so allocate() has a size
if (bondflag && specialflag == 0) {
+ if (domain->box_exist == 0)
+ error->all(FLERR,"Cannot auto-generate special bonds before "
+ "simulation box is defined");
+
maxspecial = atom->maxspecial;
if (flag) {
special_generate();
diff --git a/src/read_data.cpp b/src/read_data.cpp
index e170d909cf3986b262e1f126400187954503f86f..6ce3e434851b011e0e4e133ea43617ac719966fe 100644
--- a/src/read_data.cpp
+++ b/src/read_data.cpp
@@ -703,7 +703,11 @@ void ReadData::command(int narg, char **arg)
if (addflag == NONE) atom->deallocate_topology();
atom->avec->grow(atom->nmax);
}
+
+ // init per-atom fix/compute/variable values for created atoms
+ atom->data_fix_compute_variable(nlocal_previous,atom->nlocal);
+
// assign atoms added by this data file to specified group
if (groupbit) {
@@ -830,6 +834,7 @@ void ReadData::command(int narg, char **arg)
}
// restore old styles, when reading with nocoeff flag given
+
if (coeffflag == 0) {
if (force->pair) delete force->pair;
force->pair = saved_pair;
diff --git a/src/read_dump.cpp b/src/read_dump.cpp
index 0d8d91a32fc2469a46d6e9662f633ae6f46de5e6..8251a01cee02bc2a15d606b2c9e4ce20e71b1b32 100644
--- a/src/read_dump.cpp
+++ b/src/read_dump.cpp
@@ -908,27 +908,9 @@ void ReadDump::process_atoms(int n)
}
}
- // invoke set_arrays() for fixes/computes/variables
- // that need initialization of attributes of new atoms
- // same as in CreateAtoms
- // don't use modify->create_attributes() since would be inefficient
- // for large number of atoms
-
- nlocal = atom->nlocal;
- for (int m = 0; m < modify->nfix; m++) {
- Fix *fix = modify->fix[m];
- if (fix->create_attribute)
- for (i = nlocal_previous; i < nlocal; i++)
- fix->set_arrays(i);
- }
- for (int m = 0; m < modify->ncompute; m++) {
- Compute *compute = modify->compute[m];
- if (compute->create_attribute)
- for (i = nlocal_previous; i < nlocal; i++)
- compute->set_arrays(i);
- }
- for (int i = nlocal_previous; i < nlocal; i++)
- input->variable->set_arrays(i);
+ // init per-atom fix/compute/variable values for created atoms
+
+ atom->data_fix_compute_variable(nlocal_previous,atom->nlocal);
}
/* ----------------------------------------------------------------------
diff --git a/src/region.cpp b/src/region.cpp
index 1244b5736599ec3d21c0c33da80fc5259f9762f5..f42b4e6ad0d92dc0e322c6d625c68246a876b771 100644
--- a/src/region.cpp
+++ b/src/region.cpp
@@ -486,7 +486,7 @@ void Region::set_velocity()
else v[0] = v[1] = v[2] = 0.0;
prev[0] = dx;
prev[1] = dy;
- prev[2] = dz;
+ prev[2] = dz;
}
if (rotateflag) {
@@ -546,14 +546,13 @@ void Region::velocity_contact(double *vwall, double *x, int ic)
void Region::length_restart_string(int &n)
{
- n += sizeof(int) + strlen(id)+1 +
+ n += sizeof(int) + strlen(id)+1 +
sizeof(int) + strlen(style)+1 + sizeof(int) +
size_restart*sizeof(double);
}
/* ----------------------------------------------------------------------
- region writes its current style, id, number of sub-regions
- and position/angle
+ region writes its current style, id, number of sub-regions, position/angle
needed by fix/wall/gran/region to compute velocity by differencing scheme
------------------------------------------------------------------------- */
@@ -562,50 +561,36 @@ void Region::write_restart(FILE *fp)
int sizeid = (strlen(id)+1);
int sizestyle = (strlen(style)+1);
fwrite(&sizeid, sizeof(int), 1, fp);
- fwrite(id, 1, sizeid, fp);
- fwrite(&sizestyle, sizeof(int), 1, fp);
- fwrite(style, 1, sizestyle, fp);
+ fwrite(id,1,sizeid,fp);
+ fwrite(&sizestyle,sizeof(int),1,fp);
+ fwrite(style,1,sizestyle,fp);
fwrite(&nregion,sizeof(int),1,fp);
-
- fwrite(prev, sizeof(double), size_restart, fp);
+ fwrite(prev,sizeof(double),size_restart,fp);
}
/* ----------------------------------------------------------------------
region reads style, id, number of sub-regions from restart file
- if they match current region, also read previous position/angle
+ if they match current region, also read previous position/angle
needed by fix/wall/gran/region to compute velocity by differencing scheme
------------------------------------------------------------------------- */
int Region::restart(char *buf, int &n)
{
- int sizeid = buf[n];
+ int size = *((int *) (&buf[n]));
n += sizeof(int);
- char *restart_id = new char[sizeid];
- for (int i = 0; i < sizeid; i++)
- restart_id[i] = buf[n++];
- if (strcmp(restart_id,id) != 0) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],id) != 0)) return 0;
+ n += size;
- int sizestyle = buf[n];
+ size = *((int *) (&buf[n]));
n += sizeof(int);
- char *restart_style = new char[sizestyle];
- for (int i = 0; i < sizestyle; i++)
- restart_style[i] = buf[n++];
- if (strcmp(restart_style,style) != 0) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],style) != 0)) return 0;
+ n += size;
- int restart_nreg = buf[n];
+ int restart_nreg = *((int *) (&buf[n]));
n += sizeof(int);
if (restart_nreg != nregion) return 0;
- char *rlist = new char[size_restart*sizeof(double)];
- for (int i = 0; i < size_restart*sizeof(double); i++)
- rlist[i] = buf[n++];
- for (int i = 0; i < size_restart; i++){
- prev[i] = ((double *)rlist)[i];
- }
-
- delete [] rlist;
- delete [] restart_id;
- delete [] restart_style;
+ memcpy(prev,&buf[n],size_restart*sizeof(double));
return 1;
}
diff --git a/src/region.h b/src/region.h
index c813772a82c4b9b7ae07d9f46370136484471ddd..9c693bfcd5cc954ef2d5276c8b6fe9324284931d 100644
--- a/src/region.h
+++ b/src/region.h
@@ -60,7 +60,7 @@ class Region : protected Pointers {
double omega[3]; // angular velocity
double rprev; // speed of time-dependent radius, if applicable
double xcenter[3]; // translated/rotated center of cylinder/sphere (only used if varshape)
- double prev[5]; // stores displacement (X3), angle and if
+ double prev[5]; // stores displacement (X3), angle and if
// necessary, region variable size (e.g. radius)
// at previous time step
int vel_timestep; // store timestep at which set_velocity was called
@@ -106,11 +106,12 @@ class Region : protected Pointers {
void options(int, char **);
void point_on_line_segment(double *, double *, double *, double *);
void forward_transform(double &, double &, double &);
+ double point[3],runit[3];
private:
char *xstr,*ystr,*zstr,*tstr;
int xvar,yvar,zvar,tvar;
- double axis[3],point[3],runit[3];
+ double axis[3];
void inverse_transform(double &, double &, double &);
void rotate(double &, double &, double &, double);
diff --git a/src/region_block.cpp b/src/region_block.cpp
index f48aa5af89bf27b3ca3dc9254906504b0046a259..f1787d2ca2f2ee403ff880120315461bb799f5ad 100644
--- a/src/region_block.cpp
+++ b/src/region_block.cpp
@@ -344,7 +344,7 @@ int RegBlock::surface_exterior(double *x, double cutoff)
store closest point in xc,yc,zc
--------------------------------------------------------------------------*/
-double RegBlock::find_closest_point(int i, double *x,
+double RegBlock::find_closest_point(int i, double *x,
double &xc, double &yc, double &zc)
{
double dot,d2,d2min;
@@ -372,7 +372,7 @@ double RegBlock::find_closest_point(int i, double *x,
} else {
point_on_line_segment(corners[i][0],corners[i][1],x,p);
- d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
+ d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
(p[2]-x[2])*(p[2]-x[2]);
if (d2 < d2min) {
d2min = d2;
@@ -382,7 +382,7 @@ double RegBlock::find_closest_point(int i, double *x,
}
point_on_line_segment(corners[i][1],corners[i][2],x,p);
- d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
+ d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
(p[2]-x[2])*(p[2]-x[2]);
if (d2 < d2min) {
d2min = d2;
@@ -392,7 +392,7 @@ double RegBlock::find_closest_point(int i, double *x,
}
point_on_line_segment(corners[i][2],corners[i][3],x,p);
- d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
+ d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
(p[2]-x[2])*(p[2]-x[2]);
if (d2 < d2min) {
d2min = d2;
@@ -402,7 +402,7 @@ double RegBlock::find_closest_point(int i, double *x,
}
point_on_line_segment(corners[i][3],corners[i][4],x,p);
- d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
+ d2 = (p[0]-x[0])*(p[0]-x[0]) + (p[1]-x[1])*(p[1]-x[1]) +
(p[2]-x[2])*(p[2]-x[2]);
if (d2 < d2min) {
d2min = d2;
diff --git a/src/region_cone.cpp b/src/region_cone.cpp
index dce0e64b0fd39bc7b057a168a7aa705523b8cbc9..2af39f93a933f210acd712bda94e93fccebb328f 100644
--- a/src/region_cone.cpp
+++ b/src/region_cone.cpp
@@ -416,7 +416,7 @@ int RegCone::surface_exterior(double *x, double cutoff)
// x is far enough from cone that there is no contact
// x is interior to cone
- if (r >= maxradius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff)
+ if (r >= maxradius+cutoff || x[0] <= lo-cutoff || x[0] >= hi+cutoff)
return 0;
if (r < currentradius && x[0] > lo && x[0] < hi) return 0;
diff --git a/src/region_cylinder.cpp b/src/region_cylinder.cpp
index 6ad70dc8ed93863f4791d69fef582010a72b2fc9..5e666e472f610b6fa052e5f4ecac343775947891 100644
--- a/src/region_cylinder.cpp
+++ b/src/region_cylinder.cpp
@@ -36,7 +36,7 @@ RegCylinder::RegCylinder(LAMMPS *lmp, int narg, char **arg) :
// check open face settings
- if (openflag && (open_faces[2] || open_faces[3] ||
+ if (openflag && (open_faces[2] || open_faces[3] ||
open_faces[4] || open_faces[5]))
error->all(FLERR,"Invalid region cylinder open setting");
@@ -691,16 +691,16 @@ void RegCylinder::variable_check()
/* ----------------------------------------------------------------------
- Set values needed to calculate velocity due to shape changes.
+ Set values needed to calculate velocity due to shape changes.
These values do not depend on the contact, so this function is
called once per timestep by fix/wall/gran/region.
-
+
------------------------------------------------------------------------- */
void RegCylinder::set_velocity_shape()
{
if (axis == 'x'){
- xcenter[0] = 0;
+ xcenter[0] = 0;
xcenter[1] = c1;
xcenter[2] = c2;
}
@@ -715,15 +715,15 @@ void RegCylinder::set_velocity_shape()
xcenter[2] = 0;
}
forward_transform(xcenter[0], xcenter[1], xcenter[2]);
- if (update->ntimestep > 0) rprev = prev[4];
+ if (update->ntimestep > 0) rprev = prev[4];
else rprev = radius;
- prev[4] = radius;
+ prev[4] = radius;
}
/* ----------------------------------------------------------------------
- add velocity due to shape change to wall velocity
+ add velocity due to shape change to wall velocity
------------------------------------------------------------------------- */
void RegCylinder::velocity_contact_shape(double *vwall, double *xc)
@@ -746,7 +746,7 @@ void RegCylinder::velocity_contact_shape(double *vwall, double *xc)
}
vwall[0] += delx/update->dt;
vwall[1] += dely/update->dt;
- vwall[2] += delz/update->dt;
+ vwall[2] += delz/update->dt;
//printf ("R is %g, prev %g, velocity of wall at %g %g %g is %g %g %g\n",radius,rprev,xc[0],xc[1],xc[2],vwall[0],vwall[1],vwall[2]);
}
diff --git a/src/region_intersect.cpp b/src/region_intersect.cpp
index 02e3769ca683d662e8ffb9d29c3fa1a5f8427c06..468bd2867041d7c480c2546ee508c5c1d918d8a7 100644
--- a/src/region_intersect.cpp
+++ b/src/region_intersect.cpp
@@ -43,7 +43,7 @@ RegIntersect::RegIntersect(LAMMPS *lmp, int narg, char **arg) :
idsub[nregion] = new char[m];
strcpy(idsub[nregion],arg[iarg+3]);
iregion = domain->find_region(idsub[nregion]);
- if (iregion == -1)
+ if (iregion == -1)
error->all(FLERR,"Region intersect region ID does not exist");
list[nregion++] = iregion;
}
@@ -124,7 +124,7 @@ void RegIntersect::init()
int iregion;
for (int ilist = 0; ilist < nregion; ilist++) {
iregion = domain->find_region(idsub[ilist]);
- if (iregion == -1)
+ if (iregion == -1)
error->all(FLERR,"Region union region ID does not exist");
list[ilist] = iregion;
}
@@ -290,7 +290,7 @@ void RegIntersect::set_velocity()
void RegIntersect::length_restart_string(int& n)
{
- n += sizeof(int) + strlen(id)+1 +
+ n += sizeof(int) + strlen(id)+1 +
sizeof(int) + strlen(style)+1 + sizeof(int);
for (int ilist = 0; ilist < nregion; ilist++)
domain->regions[list[ilist]]->length_restart_string(n);
@@ -308,7 +308,7 @@ void RegIntersect::write_restart(FILE *fp)
fwrite(&sizeid, sizeof(int), 1, fp);
fwrite(id, 1, sizeid, fp);
fwrite(&sizestyle, sizeof(int), 1, fp);
- fwrite(style, 1, sizestyle, fp);
+ fwrite(style, 1, sizestyle, fp);
fwrite(&nregion,sizeof(int),1,fp);
for (int ilist = 0; ilist < nregion; ilist++){
@@ -321,33 +321,26 @@ void RegIntersect::write_restart(FILE *fp)
needed by fix/wall/gran/region to compute velocity by differencing scheme
------------------------------------------------------------------------- */
-int RegIntersect::restart(char *buf, int& n)
+int RegIntersect::restart(char *buf, int &n)
{
- int sizeid = buf[n];
+ int size = *((int *) (&buf[n]));
n += sizeof(int);
- char *restart_id = new char[sizeid];
- for (int i = 0; i < sizeid; i++)
- restart_id[i] = buf[n++];
- if (strcmp(restart_id,id) != 0) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],id) != 0)) return 0;
+ n += size;
- int sizestyle = buf[n];
+ size = *((int *) (&buf[n]));
n += sizeof(int);
+ if ((size <= 0) || (strcmp(&buf[n],style) != 0)) return 0;
+ n += size;
- char *restart_style = new char[sizestyle];
- for (int i = 0; i < sizestyle; i++)
- restart_style[i] = buf[n++];
- if (strcmp(restart_style,style) != 0) return 0;
-
- int restart_nreg = buf[n];
+ int restart_nreg = *((int *) (&buf[n]));
n += sizeof(int);
if (restart_nreg != nregion) return 0;
- for (int ilist = 0; ilist < nregion; ilist++){
- if (!domain->regions[list[ilist]]->restart(buf, n)){
- return 0;
- }
- }
- return 1;
+ for (int ilist = 0; ilist < nregion; ilist++)
+ if (!domain->regions[list[ilist]]->restart(buf,n)) return 0;
+
+ return 1;
}
/* ----------------------------------------------------------------------
diff --git a/src/region_sphere.cpp b/src/region_sphere.cpp
index 2f51baacc4c1f63d0e680fac5cf87415c97df32c..8aef087d38eab60b8d46b910951797d1405c36bc 100644
--- a/src/region_sphere.cpp
+++ b/src/region_sphere.cpp
@@ -188,10 +188,10 @@ void RegSphere::variable_check()
/* ----------------------------------------------------------------------
- Set values needed to calculate velocity due to shape changes.
+ Set values needed to calculate velocity due to shape changes.
These values do not depend on the contact, so this function is
called once per timestep by fix/wall/gran/region.
-
+
------------------------------------------------------------------------- */
void RegSphere::set_velocity_shape()
@@ -200,27 +200,27 @@ void RegSphere::set_velocity_shape()
xcenter[1] = yc;
xcenter[2] = zc;
forward_transform(xcenter[0], xcenter[1], xcenter[2]);
- if (update->ntimestep > 0) rprev = prev[4];
+ if (update->ntimestep > 0) rprev = prev[4];
else rprev = radius;
- prev[4] = radius;
+ prev[4] = radius;
}
/* ----------------------------------------------------------------------
- add velocity due to shape change to wall velocity
+ add velocity due to shape change to wall velocity
------------------------------------------------------------------------- */
void RegSphere::velocity_contact_shape(double *vwall, double *xc)
{
double delx, dely, delz; // Displacement of contact point in x,y,z
-
+
delx = (xc[0] - xcenter[0])*(1 - rprev/radius);
dely = (xc[1] - xcenter[1])*(1 - rprev/radius);
delz = (xc[2] - xcenter[2])*(1 - rprev/radius);
vwall[0] += delx/update->dt;
vwall[1] += dely/update->dt;
- vwall[2] += delz/update->dt;
+ vwall[2] += delz/update->dt;
}
diff --git a/src/region_union.cpp b/src/region_union.cpp
index c506c811c300751266ab0b94dcb74a59ecd750b5..51bb1b339de541104d2f46a6b30a51f005dd1e92 100644
--- a/src/region_union.cpp
+++ b/src/region_union.cpp
@@ -44,7 +44,7 @@ RegUnion::RegUnion(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
idsub[nregion] = new char[m];
strcpy(idsub[nregion],arg[iarg+3]);
iregion = domain->find_region(idsub[nregion]);
- if (iregion == -1)
+ if (iregion == -1)
error->all(FLERR,"Region union region ID does not exist");
list[nregion++] = iregion;
}
@@ -118,7 +118,7 @@ void RegUnion::init()
int iregion;
for (int ilist = 0; ilist < nregion; ilist++) {
iregion = domain->find_region(idsub[ilist]);
- if (iregion == -1)
+ if (iregion == -1)
error->all(FLERR,"Region union region ID does not exist");
list[ilist] = iregion;
}
@@ -171,7 +171,7 @@ int RegUnion::surface_interior(double *x, double cutoff)
for (jlist = 0; jlist < nregion; jlist++) {
if (jlist == ilist) continue;
jregion = list[jlist];
- if (regions[jregion]->match(xs,ys,zs) &&
+ if (regions[jregion]->match(xs,ys,zs) &&
!regions[jregion]->openflag) break;
}
if (jlist == nregion) {
@@ -284,7 +284,7 @@ void RegUnion::set_velocity()
void RegUnion::length_restart_string(int& n)
{
- n += sizeof(int) + strlen(id)+1 +
+ n += sizeof(int) + strlen(id)+1 +
sizeof(int) + strlen(style)+1 + sizeof(int);
for (int ilist = 0; ilist < nregion; ilist++)
domain->regions[list[ilist]]->length_restart_string(n);
@@ -302,10 +302,10 @@ void RegUnion::write_restart(FILE *fp)
fwrite(&sizeid, sizeof(int), 1, fp);
fwrite(id, 1, sizeid, fp);
fwrite(&sizestyle, sizeof(int), 1, fp);
- fwrite(style, 1, sizestyle, fp);
+ fwrite(style, 1, sizestyle, fp);
fwrite(&nregion,sizeof(int),1,fp);
for (int ilist = 0; ilist < nregion; ilist++)
- domain->regions[list[ilist]]->write_restart(fp);
+ domain->regions[list[ilist]]->write_restart(fp);
}
/* ----------------------------------------------------------------------
@@ -315,27 +315,23 @@ void RegUnion::write_restart(FILE *fp)
int RegUnion::restart(char *buf, int &n)
{
- int sizeid = buf[n];
+ int size = *((int *) (&buf[n]));
n += sizeof(int);
- char *restart_id = new char[sizeid];
- for (int i = 0; i < sizeid; i++)
- restart_id[i] = buf[n++];
- if (strcmp(restart_id,id) != 0) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],id) != 0)) return 0;
+ n += size;
- int sizestyle = buf[n];
+ size = *((int *) (&buf[n]));
n += sizeof(int);
- char *restart_style = new char[sizestyle];
- for (int i = 0; i < sizestyle; i++)
- restart_style[i] = buf[n++];
- if (strcmp(restart_style,style) != 0) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],style) != 0)) return 0;
+ n += size;
- int restart_nreg = buf[n];
+ int restart_nreg = *((int *) (&buf[n]));
n += sizeof(int);
if (restart_nreg != nregion) return 0;
- for (int ilist = 0; ilist < nregion; ilist++){
- if (!domain->regions[list[ilist]]->restart(buf, n)) return 0;
- }
+ for (int ilist = 0; ilist < nregion; ilist++)
+ if (!domain->regions[list[ilist]]->restart(buf,n)) return 0;
+
return 1;
}
diff --git a/src/verlet.cpp b/src/verlet.cpp
index ce06a1c45688c03fa6281b195065061eab0d60ba..035cab79ef6f124d5189f544ffeb17ee9c011159 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -95,6 +95,9 @@ void Verlet::setup()
timer->print_timeout(screen);
}
+ if (lmp->kokkos)
+ error->all(FLERR,"KOKKOS package requires run_style verlet/kk");
+
update->setupflag = 1;
// setup domain, communication and neighboring
diff --git a/src/verlet.h b/src/verlet.h
index 80e46f623668b379b095a6f4a6a1479da44f465a..0e2a333fabc76ad1b32e794db741af67e8dab566 100644
--- a/src/verlet.h
+++ b/src/verlet.h
@@ -53,4 +53,9 @@ W: No fixes defined, atoms won't move
If you are not using a fix like nve, nvt, npt then atom velocities and
coordinates will not be updated during timestepping.
+E: KOKKOS package requires run_style verlet/kk
+
+The KOKKOS package requires the Kokkos version of run_style verlet; the
+regular version cannot be used.
+
*/
diff --git a/src/version.h b/src/version.h
index 3cdddc4efdb039cccbe0fa7cad7126ac19df4c98..c6e869ccdedad1656b6dc5e20bba32742220b266 100644
--- a/src/version.h
+++ b/src/version.h
@@ -1 +1 @@
-#define LAMMPS_VERSION "19 Oct 2016"
+#define LAMMPS_VERSION "27 Oct 2016"