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/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_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/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/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/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/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_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/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 a37b81dbd5f16cc910295ccbe1a35a33de1763a9..f42b4e6ad0d92dc0e322c6d625c68246a876b771 100644
--- a/src/region.cpp
+++ b/src/region.cpp
@@ -552,8 +552,7 @@ void Region::length_restart_string(int &n)
}
/* ----------------------------------------------------------------------
- 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,37 +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 size = *((int *)(buf+n));
+ int size = *((int *) (&buf[n]));
n += sizeof(int);
- if ((size <= 0) || (strcmp(buf+n,id) != 0)) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],id) != 0)) return 0;
n += size;
- size = *((int *)(buf+n));
+ size = *((int *) (&buf[n]));
n += sizeof(int);
- if ((size <= 0) || (strcmp(buf+n,style) != 0)) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],style) != 0)) return 0;
n += size;
- int restart_nreg = *((int *)(buf+n));
+ int restart_nreg = *((int *) (&buf[n]));
n += sizeof(int);
if (restart_nreg != nregion) return 0;
- memcpy(prev,buf+n,size_restart*sizeof(double));
+ memcpy(prev,&buf[n],size_restart*sizeof(double));
return 1;
}
diff --git a/src/region.h b/src/region.h
index ecc285d6c6a12a6ebe0669fe86e7ecf0e3aedbac..9c693bfcd5cc954ef2d5276c8b6fe9324284931d 100644
--- a/src/region.h
+++ b/src/region.h
@@ -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_intersect.cpp b/src/region_intersect.cpp
index 03c9ccc360ef41b2cb0a05ec216cb59a643a4c96..468bd2867041d7c480c2546ee508c5c1d918d8a7 100644
--- a/src/region_intersect.cpp
+++ b/src/region_intersect.cpp
@@ -321,24 +321,24 @@ 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 size = *((int *)(buf+n));
+ int size = *((int *) (&buf[n]));
n += sizeof(int);
- if ((size <= 0) || (strcmp(buf+n,id) != 0)) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],id) != 0)) return 0;
n += size;
- size = *((int *)(buf+n));
+ size = *((int *) (&buf[n]));
n += sizeof(int);
- if ((size <= 0) || (strcmp(buf+n,style) != 0)) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],style) != 0)) return 0;
n += size;
- int restart_nreg = *((int *)(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;
+ if (!domain->regions[list[ilist]]->restart(buf,n)) return 0;
return 1;
}
diff --git a/src/region_union.cpp b/src/region_union.cpp
index e1b86e159e4204a7e4ee833c7ab882bc01b526dc..51bb1b339de541104d2f46a6b30a51f005dd1e92 100644
--- a/src/region_union.cpp
+++ b/src/region_union.cpp
@@ -315,22 +315,22 @@ void RegUnion::write_restart(FILE *fp)
int RegUnion::restart(char *buf, int &n)
{
- int size = *((int *)(buf+n));
+ int size = *((int *) (&buf[n]));
n += sizeof(int);
- if ((size <= 0) || (strcmp(buf+n,id) != 0)) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],id) != 0)) return 0;
n += size;
- size = *((int *)(buf+n));
+ size = *((int *) (&buf[n]));
n += sizeof(int);
- if ((size <= 0) || (strcmp(buf+n,style) != 0)) return 0;
+ if ((size <= 0) || (strcmp(&buf[n],style) != 0)) return 0;
n += size;
- int restart_nreg = *((int *)(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;
+ 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"