jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

new library functions

This commit is contained in:
Steve Plimpton 2016-10-27 09:34:04 -06:00
parent efaa8feab5
commit d9891abdf4
44 changed files with 1341 additions and 211 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
doc/src/Manual.txt
View File

@ -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

86
doc/src/Section_howto.txt
View File

@ -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

8
doc/src/Section_python.txt
View File

@ -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

2
examples/COUPLE/simple/in.lj
View File

@ -20,4 +20,6 @@ neigh_modify delay 0 every 20 check no
fix 1 all nve
variable fx atom fx
run 10

73
examples/COUPLE/simple/simple.c
View File

@ -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();
}

61
examples/COUPLE/simple/simple.cpp
View File

@ -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();
}

24
python/examples/simple.py
View File

@ -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()

60
python/lammps.py
View File

@ -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:

4
src/KOKKOS/Install.sh
View File

@ -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

2
src/KOKKOS/fix_qeq_reax_kokkos.cpp
View File

@ -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);

124
src/KOKKOS/fix_reaxc_bonds_kokkos.cpp Normal file
View File

@ -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;
}

42
src/KOKKOS/fix_reaxc_bonds_kokkos.h Normal file
View File

@ -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

157
src/KOKKOS/fix_reaxc_species_kokkos.cpp Normal file
View File

@ -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;
}
}

41
src/KOKKOS/fix_reaxc_species_kokkos.h Normal file
View File

@ -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

3
src/KOKKOS/neighbor_kokkos.cpp
View File

@ -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;

4
src/KOKKOS/neighbor_kokkos.h
View File

@ -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

8
src/KOKKOS/pair_kokkos.h
View File

@ -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 >

242
src/KOKKOS/pair_reax_c_kokkos.cpp
View File

@ -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>;

72
src/KOKKOS/pair_reax_c_kokkos.h
View File

@ -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);
}
};
}

2
src/USER-REAXC/compute_spec_atom.cpp
View File

@ -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];

7
src/USER-REAXC/fix_reaxc_bonds.cpp
View File

@ -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");

10
src/USER-REAXC/fix_reaxc_bonds.h
View File

@ -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;

9
src/USER-REAXC/fix_reaxc_species.cpp
View File

@ -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;

9
src/USER-REAXC/fix_reaxc_species.h
View File

@ -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);

30
src/atom.cpp
View File

@ -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

5
src/atom.h
View File

@ -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);

27
src/create_atoms.cpp
View File

@ -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

57
src/delete_atoms.cpp
View File

@ -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;

2
src/delete_atoms.h
View File

@ -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 **);

24
src/domain.cpp
View File

@ -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]);
}

2
src/domain.h
View File

@ -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 **);

11
src/fix_nh.cpp
View File

@ -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;
}

3
src/input.cpp
View File

@ -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
------------------------------------------------------------------------- */

238
src/library.cpp
View File

@ -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
}
#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

5
src/library.h
View File

@ -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 *);

5
src/read_data.cpp
View File

@ -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;

24
src/read_dump.cpp
View File

@ -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);
}
/* ----------------------------------------------------------------------

26
src/region.cpp
View File

@ -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;
}

3
src/region.h
View File

@ -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);

14
src/region_intersect.cpp
View File

@ -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;
}

12
src/region_union.cpp
View File

@ -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;
}

3
src/verlet.cpp
View File

@ -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

5
src/verlet.h
View File

@ -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.
*/

2
src/version.h
View File

@ -1 +1 @@
#define LAMMPS_VERSION "19 Oct 2016"
#define LAMMPS_VERSION "27 Oct 2016"