diff --git a/doc/src/Manual.txt b/doc/src/Manual.txt index a957470735..8590e33ace 100644 --- a/doc/src/Manual.txt +++ b/doc/src/Manual.txt @@ -1,7 +1,7 @@ LAMMPS Users Manual - + @@ -21,7 +21,7 @@

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 7b8ed2d9a6..33cbfa9586 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. +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. -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_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 -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 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). + +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 352c1fa592..b3b5171e74 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 5f6ebdc4f2..3dc80bdfea 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 c50f289b7a..cc813d78fe 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 */ - if (lammps == 1) { - lammps_command(ptr,"run 10"); + double *x = NULL; + double *v = NULL; - int natoms = lammps_get_natoms(ptr); - double *x = (double *) malloc(3*natoms*sizeof(double)); - lammps_gather_atoms(ptr,"x",1,3,x); + if (lammps == 1) { + lammps_command(lmp,"run 10"); + + 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(ptr,"run 1"); + lammps_command(lmp,"run 1"); } - if (lammps == 1) lammps_close(ptr); + // 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 = 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 (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 a440bcb1ae..b279ae704c 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 (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 (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 9315a555fc..6336c89de3 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 8c57928064..7344a16be3 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 c45be8c973..41f9304f8b 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 e54be74124..d39d5cc84a 100644 --- a/src/KOKKOS/fix_qeq_reax_kokkos.cpp +++ b/src/KOKKOS/fix_qeq_reax_kokkos.cpp @@ -431,8 +431,8 @@ void FixQEqReaxKokkos::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 0000000000..6e96c49196 --- /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 +#include +#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 (atom->natoms); + + numbonds = 0; + if (reaxc->execution_space == Device) + ((PairReaxCKokkos*) reaxc)->FindBond(numbonds); + else + ((PairReaxCKokkos*) 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*) reaxc)->PackBondBuffer(k_buf,nbuf_local); + else + ((PairReaxCKokkos*) 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 0000000000..f27e558fdd --- /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 0000000000..d3de0c998c --- /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 +#include +#include "atom.h" +#include +#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 0000000000..64de060006 --- /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 3bb8b0dce4..9622129e70 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 3693ebd9a4..8c097139a7 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 b880654fbf..3710c460c0 100644 --- a/src/KOKKOS/pair_kokkos.h +++ b/src/KOKKOS/pair_kokkos.h @@ -611,7 +611,7 @@ struct PairComputeFunctor { // pair_compute_neighlist and pair_compute_fullcluster will match - either the dummy version // or the real one further below. template -EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if*>::type list) { +EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if*>::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 -EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if*>::type list) { +EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if*>::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 -EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if*>::type list) { +EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable_if<(NEIGHFLAG&PairStyle::EnabledNeighFlags) != 0, NeighListKokkos*>::type list) { EV_FLOAT ev; if(fpair->atom->ntypes > MAX_TYPES_STACKPARAMS) { PairComputeFunctor ff(fpair,list); @@ -646,7 +646,7 @@ EV_FLOAT pair_compute_neighlist (PairStyle* fpair, typename Kokkos::Impl::enable // Submit ParallelFor for NEIGHFLAG=FULLCLUSTER template -EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if*>::type list) { +EV_FLOAT pair_compute_fullcluster (PairStyle* fpair, typename Kokkos::Impl::enable_if<(FULLCLUSTER&PairStyle::EnabledNeighFlags) != 0, NeighListKokkos*>::type list) { EV_FLOAT ev; if(fpair->atom->ntypes > MAX_TYPES_STACKPARAMS) { typedef PairComputeFunctor diff --git a/src/KOKKOS/pair_reax_c_kokkos.cpp b/src/KOKKOS/pair_reax_c_kokkos.cpp index 4fcb3652e3..894c3ab53c 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::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::PairReaxCKokkos(LAMMPS *lmp) : PairReaxC(lmp) template PairReaxCKokkos::~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::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::compute(int eflag_in, int vflag_in) k_vatom.template sync(); } + if (fixspecies_flag) + FindBondSpecies(); + copymode = 0; } @@ -1346,6 +1362,19 @@ void PairReaxCKokkos::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::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 +void PairReaxCKokkos::FindBond(int &numbonds) +{ + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,nmax),*this); + DeviceType::fence(); + + bo_cut_bond = control->bg_cut; + + atomKK->sync(execution_space,TAG_MASK); + tag = atomKK->k_tag.view(); + + const int inum = list->inum; + NeighListKokkos* k_list = static_cast*>(list); + d_ilist = k_list->d_ilist; + k_list->clean_copy(); + + numbonds = 0; + PairReaxCKokkosFindBondFunctor find_bond_functor(this); + Kokkos::parallel_reduce(inum,find_bond_functor,numbonds); + DeviceType::fence(); + copymode = 0; +} + +template +KOKKOS_INLINE_FUNCTION +void PairReaxCKokkos::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 +KOKKOS_INLINE_FUNCTION +void PairReaxCKokkos::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 +void PairReaxCKokkos::PackBondBuffer(DAT::tdual_ffloat_1d k_buf, int &nbuf_local) +{ + d_buf = k_buf.view(); + k_params_sing.template sync(); + atomKK->sync(execution_space,TAG_MASK|TYPE_MASK|Q_MASK|MOLECULE_MASK); + + tag = atomKK->k_tag.view(); + type = atomKK->k_type.view(); + q = atomKK->k_q.view(); + if (atom->molecule) + molecule = atomKK->k_molecule.view(); + + copymode = 1; + nlocal = atomKK->nlocal; + PairReaxCKokkosPackBondBufferFunctor pack_bond_buffer_functor(this); + Kokkos::parallel_scan(nlocal,pack_bond_buffer_functor); + DeviceType::fence(); + copymode = 0; + + k_buf.modify(); + k_nbuf_local.modify(); + + k_buf.sync(); + k_nbuf_local.sync(); + nbuf_local = k_nbuf_local.h_view(); +} + +template +KOKKOS_INLINE_FUNCTION +void PairReaxCKokkos::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()() = j - 1; +} + +/* ---------------------------------------------------------------------- */ + +template +void PairReaxCKokkos::FindBondSpecies() +{ + copymode = 1; + Kokkos::parallel_for(Kokkos::RangePolicy(0,nmax),*this); + DeviceType::fence(); + + nlocal = atomKK->nlocal; + Kokkos::parallel_for(Kokkos::RangePolicy(0,nlocal),*this); + DeviceType::fence(); + copymode = 0; + + // NOTE: Could improve performance if a Kokkos version of ComputeSpecAtom is added + + k_tmpbo.modify(); + k_tmpid.modify(); + k_error_flag.modify(); + + k_tmpbo.sync(); + k_tmpid.sync(); + k_error_flag.sync(); + + if (k_error_flag.h_view()) + error->all(FLERR,"Increase MAXSPECBOND in reaxc_defs.h"); +} + +template +KOKKOS_INLINE_FUNCTION +void PairReaxCKokkos::operator()(PairReaxFindBondSpeciesZero, const int &i) const { + for (int j = 0; j < MAXSPECBOND; j++) { + k_tmpbo.view()(i,j) = 0.0; + k_tmpid.view()(i,j) = 0; + } +} + +template +KOKKOS_INLINE_FUNCTION +void PairReaxCKokkos::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()(i,nj) = j; + k_tmpbo.view()(i,nj) = bo_tmp; + nj++; + if (nj > MAXSPECBOND) k_error_flag.view()() = 1; + } + } +} + template class PairReaxCKokkos; #ifdef KOKKOS_HAVE_CUDA template class PairReaxCKokkos; diff --git a/src/KOKKOS/pair_reax_c_kokkos.h b/src/KOKKOS/pair_reax_c_kokkos.h index 8c07ee2a03..204ad7732f 100644 --- a/src/KOKKOS/pair_reax_c_kokkos.h +++ b/src/KOKKOS/pair_reax_c_kokkos.h @@ -115,6 +115,13 @@ struct PairReaxComputeTorsion{}; template struct PairReaxComputeHydrogen{}; +struct PairReaxFindBondZero{}; + +struct PairReaxFindBondSpeciesZero{}; + +struct PairReaxFindBondSpecies{}; + + template 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 KOKKOS_INLINE_FUNCTION @@ -245,6 +255,21 @@ class PairReaxCKokkos : public PairReaxC { KOKKOS_INLINE_FUNCTION void operator()(PairReaxComputeHydrogen, 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 +struct PairReaxCKokkosFindBondFunctor { + typedef DeviceType device_type; + typedef int value_type; + PairReaxCKokkos c; + PairReaxCKokkosFindBondFunctor(PairReaxCKokkos* 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 +struct PairReaxCKokkosPackBondBufferFunctor { + typedef DeviceType device_type; + typedef int value_type; + PairReaxCKokkos c; + PairReaxCKokkosPackBondBufferFunctor(PairReaxCKokkos* 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 38cf9ad6e6..4af8efcae7 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 5f653557a5..543669de76 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 b91085163b..e1dcd82c0e 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 c8acc5f0e0..ead73f02a1 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 a3a47b29dd..872ea2528f 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 762341e002..053a18430b 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 00149c5dae..9abbb49569 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, const char *, int); void set_mass(const char *, int, int, double); diff --git a/src/create_atoms.cpp b/src/create_atoms.cpp index 82315e7dc5..776a94c8b9 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 4ac7c02622..116c2a2f1e 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 + + 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); } - - 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 6fc2ea1f90..62ba47d715 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 *hash; void delete_group(int, char **); diff --git a/src/domain.cpp b/src/domain.cpp index 2c3f61401d..c348167dfb 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 ef3de42071..e2425f4585 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 01ba81e554..1799a00558 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 7a46929445..c1279a967d 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 44e2c40782..90c7cfa9d5 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -18,9 +18,11 @@ #include #include #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 } -#ifdef LAMMPS_EXCEPTIONS /* ---------------------------------------------------------------------- - Check if a new error message + 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 ------------------------------------------------------------------------- */ 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 e4dffc63d1..a05b7b4eb8 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 7759b9b03d..d6a33d6e9d 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 0d8d91a32f..8251a01cee 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 a37b81dbd5..f42b4e6ad0 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 ecc285d6c6..9c693bfcd5 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 03c9ccc360..468bd28670 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 e1b86e159e..51bb1b339d 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 ce06a1c456..035cab79ef 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 80e46f6236..0e2a333fab 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 3cdddc4efd..c6e869ccde 100644 --- a/src/version.h +++ b/src/version.h @@ -1 +1 @@ -#define LAMMPS_VERSION "19 Oct 2016" +#define LAMMPS_VERSION "27 Oct 2016"