diff --git a/doc/src/Manual.txt b/doc/src/Manual.txt
index 5677d09a8..25319570a 100644
--- a/doc/src/Manual.txt
+++ b/doc/src/Manual.txt
@@ -1,340 +1,340 @@
LAMMPS Users Manual
-
+
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
LAMMPS Documentation :c,h3
-20 Jan 2017 version :c,h4
+26 Jan 2017 version :c,h4
Version info: :h4
The LAMMPS "version" is the date when it was released, such as 1 May
2010. LAMMPS is updated continuously. Whenever we fix a bug or add a
feature, we release it immediately, and post a notice on "this page of
the WWW site"_bug. Every 2-4 months one of the incremental releases
is subjected to more thorough testing and labeled as a {stable} version.
Each dated copy of LAMMPS contains all the
features and bug-fixes up to and including that version date. The
version date is printed to the screen and logfile every time you run
LAMMPS. It is also in the file src/version.h and in the LAMMPS
directory name created when you unpack a tarball, and at the top of
the first page of the manual (this page).
If you browse the HTML doc pages on the LAMMPS WWW site, they always
describe the most current version of LAMMPS. :ulb,l
If you browse the HTML doc pages included in your tarball, they
describe the version you have. :l
The "PDF file"_Manual.pdf on the WWW site or in the tarball is updated
about once per month. This is because it is large, and we don't want
it to be part of every patch. :l
There is also a "Developer.pdf"_Developer.pdf file in the doc
directory, which describes the internal structure and algorithms of
LAMMPS. :l
:ule
LAMMPS stands for Large-scale Atomic/Molecular Massively Parallel
Simulator.
LAMMPS is a classical molecular dynamics simulation code designed to
run efficiently on parallel computers. It was developed at Sandia
National Laboratories, a US Department of Energy facility, with
funding from the DOE. It is an open-source code, distributed freely
under the terms of the GNU Public License (GPL).
The current core group of LAMMPS developers is at Sandia National
Labs and Temple University:
"Steve Plimpton"_sjp, sjplimp at sandia.gov :ulb,l
Aidan Thompson, athomps at sandia.gov :l
Stan Moore, stamoore at sandia.gov :l
"Axel Kohlmeyer"_ako, akohlmey at gmail.com :l
:ule
Past core developers include Paul Crozier, Ray Shan and Mark Stevens,
all at Sandia. The [LAMMPS home page] at
"http://lammps.sandia.gov"_http://lammps.sandia.gov has more information
about the code and its uses. Interaction with external LAMMPS developers,
bug reports and feature requests are mainly coordinated through the
"LAMMPS project on GitHub."_https://github.com/lammps/lammps
The lammps.org domain, currently hosting "public continuous integration
testing"_https://ci.lammps.org/job/lammps/ and "precompiled Linux
RPM and Windows installer packages"_http://rpm.lammps.org is located
at Temple University and managed by Richard Berger,
richard.berger at temple.edu.
:link(bug,http://lammps.sandia.gov/bug.html)
:link(sjp,http://www.sandia.gov/~sjplimp)
:link(ako,http://goo.gl/1wk0)
:line
The LAMMPS documentation is organized into the following sections. If
you find errors or omissions in this manual or have suggestions for
useful information to add, please send an email to the developers so
we can improve the LAMMPS documentation.
Once you are familiar with LAMMPS, you may want to bookmark "this
page"_Section_commands.html#comm at Section_commands.html#comm since
it gives quick access to documentation for all LAMMPS commands.
"PDF file"_Manual.pdf of the entire manual, generated by
"htmldoc"_http://freecode.com/projects/htmldoc
"Introduction"_Section_intro.html :olb,l
1.1 "What is LAMMPS"_intro_1 :ulb,b
1.2 "LAMMPS features"_intro_2 :b
1.3 "LAMMPS non-features"_intro_3 :b
1.4 "Open source distribution"_intro_4 :b
1.5 "Acknowledgments and citations"_intro_5 :ule,b
"Getting started"_Section_start.html :l
2.1 "What's in the LAMMPS distribution"_start_1 :ulb,b
2.2 "Making LAMMPS"_start_2 :b
2.3 "Making LAMMPS with optional packages"_start_3 :b
2.4 "Building LAMMPS via the Make.py script"_start_4 :b
2.5 "Building LAMMPS as a library"_start_5 :b
2.6 "Running LAMMPS"_start_6 :b
2.7 "Command-line options"_start_7 :b
2.8 "Screen output"_start_8 :b
2.9 "Tips for users of previous versions"_start_9 :ule,b
"Commands"_Section_commands.html :l
3.1 "LAMMPS input script"_cmd_1 :ulb,b
3.2 "Parsing rules"_cmd_2 :b
3.3 "Input script structure"_cmd_3 :b
3.4 "Commands listed by category"_cmd_4 :b
3.5 "Commands listed alphabetically"_cmd_5 :ule,b
"Packages"_Section_packages.html :l
4.1 "Standard packages"_pkg_1 :ulb,b
4.2 "User packages"_pkg_2 :ule,b
"Accelerating LAMMPS performance"_Section_accelerate.html :l
5.1 "Measuring performance"_acc_1 :ulb,b
5.2 "Algorithms and code options to boost performace"_acc_2 :b
5.3 "Accelerator packages with optimized styles"_acc_3 :b
5.3.1 "GPU package"_accelerate_gpu.html :ulb,b
5.3.2 "USER-INTEL package"_accelerate_intel.html :b
5.3.3 "KOKKOS package"_accelerate_kokkos.html :b
5.3.4 "USER-OMP package"_accelerate_omp.html :b
5.3.5 "OPT package"_accelerate_opt.html :ule,b
5.4 "Comparison of various accelerator packages"_acc_4 :ule,b
"How-to discussions"_Section_howto.html :l
6.1 "Restarting a simulation"_howto_1 :ulb,b
6.2 "2d simulations"_howto_2 :b
6.3 "CHARMM and AMBER force fields"_howto_3 :b
6.4 "Running multiple simulations from one input script"_howto_4 :b
6.5 "Multi-replica simulations"_howto_5 :b
6.6 "Granular models"_howto_6 :b
6.7 "TIP3P water model"_howto_7 :b
6.8 "TIP4P water model"_howto_8 :b
6.9 "SPC water model"_howto_9 :b
6.10 "Coupling LAMMPS to other codes"_howto_10 :b
6.11 "Visualizing LAMMPS snapshots"_howto_11 :b
6.12 "Triclinic (non-orthogonal) simulation boxes"_howto_12 :b
6.13 "NEMD simulations"_howto_13 :b
6.14 "Finite-size spherical and aspherical particles"_howto_14 :b
6.15 "Output from LAMMPS (thermo, dumps, computes, fixes, variables)"_howto_15 :b
6.16 "Thermostatting, barostatting, and compute temperature"_howto_16 :b
6.17 "Walls"_howto_17 :b
6.18 "Elastic constants"_howto_18 :b
6.19 "Library interface to LAMMPS"_howto_19 :b
6.20 "Calculating thermal conductivity"_howto_20 :b
6.21 "Calculating viscosity"_howto_21 :b
6.22 "Calculating a diffusion coefficient"_howto_22 :b
6.23 "Using chunks to calculate system properties"_howto_23 :b
6.24 "Setting parameters for pppm/disp"_howto_24 :b
6.25 "Polarizable models"_howto_25 :b
6.26 "Adiabatic core/shell model"_howto_26 :b
6.27 "Drude induced dipoles"_howto_27 :ule,b
"Example problems"_Section_example.html :l
"Performance & scalability"_Section_perf.html :l
"Additional tools"_Section_tools.html :l
"Modifying & extending LAMMPS"_Section_modify.html :l
10.1 "Atom styles"_mod_1 :ulb,b
10.2 "Bond, angle, dihedral, improper potentials"_mod_2 :b
10.3 "Compute styles"_mod_3 :b
10.4 "Dump styles"_mod_4 :b
10.5 "Dump custom output options"_mod_5 :b
10.6 "Fix styles"_mod_6 :b
10.7 "Input script commands"_mod_7 :b
10.8 "Kspace computations"_mod_8 :b
10.9 "Minimization styles"_mod_9 :b
10.10 "Pairwise potentials"_mod_10 :b
10.11 "Region styles"_mod_11 :b
10.12 "Body styles"_mod_12 :b
10.13 "Thermodynamic output options"_mod_13 :b
10.14 "Variable options"_mod_14 :b
10.15 "Submitting new features for inclusion in LAMMPS"_mod_15 :ule,b
"Python interface"_Section_python.html :l
11.1 "Overview of running LAMMPS from Python"_py_1 :ulb,b
11.2 "Overview of using Python from a LAMMPS script"_py_2 :b
11.3 "Building LAMMPS as a shared library"_py_3 :b
11.4 "Installing the Python wrapper into Python"_py_4 :b
11.5 "Extending Python with MPI to run in parallel"_py_5 :b
11.6 "Testing the Python-LAMMPS interface"_py_6 :b
11.7 "Using LAMMPS from Python"_py_7 :b
11.8 "Example Python scripts that use LAMMPS"_py_8 :ule,b
"Errors"_Section_errors.html :l
12.1 "Common problems"_err_1 :ulb,b
12.2 "Reporting bugs"_err_2 :b
12.3 "Error & warning messages"_err_3 :ule,b
"Future and history"_Section_history.html :l
13.1 "Coming attractions"_hist_1 :ulb,b
13.2 "Past versions"_hist_2 :ule,b
:ole
:link(intro_1,Section_intro.html#intro_1)
:link(intro_2,Section_intro.html#intro_2)
:link(intro_3,Section_intro.html#intro_3)
:link(intro_4,Section_intro.html#intro_4)
:link(intro_5,Section_intro.html#intro_5)
:link(start_1,Section_start.html#start_1)
:link(start_2,Section_start.html#start_2)
:link(start_3,Section_start.html#start_3)
:link(start_4,Section_start.html#start_4)
:link(start_5,Section_start.html#start_5)
:link(start_6,Section_start.html#start_6)
:link(start_7,Section_start.html#start_7)
:link(start_8,Section_start.html#start_8)
:link(start_9,Section_start.html#start_9)
:link(cmd_1,Section_commands.html#cmd_1)
:link(cmd_2,Section_commands.html#cmd_2)
:link(cmd_3,Section_commands.html#cmd_3)
:link(cmd_4,Section_commands.html#cmd_4)
:link(cmd_5,Section_commands.html#cmd_5)
:link(pkg_1,Section_packages.html#pkg_1)
:link(pkg_2,Section_packages.html#pkg_2)
:link(acc_1,Section_accelerate.html#acc_1)
:link(acc_2,Section_accelerate.html#acc_2)
:link(acc_3,Section_accelerate.html#acc_3)
:link(acc_4,Section_accelerate.html#acc_4)
:link(howto_1,Section_howto.html#howto_1)
:link(howto_2,Section_howto.html#howto_2)
:link(howto_3,Section_howto.html#howto_3)
:link(howto_4,Section_howto.html#howto_4)
:link(howto_5,Section_howto.html#howto_5)
:link(howto_6,Section_howto.html#howto_6)
:link(howto_7,Section_howto.html#howto_7)
:link(howto_8,Section_howto.html#howto_8)
:link(howto_9,Section_howto.html#howto_9)
:link(howto_10,Section_howto.html#howto_10)
:link(howto_11,Section_howto.html#howto_11)
:link(howto_12,Section_howto.html#howto_12)
:link(howto_13,Section_howto.html#howto_13)
:link(howto_14,Section_howto.html#howto_14)
:link(howto_15,Section_howto.html#howto_15)
:link(howto_16,Section_howto.html#howto_16)
:link(howto_17,Section_howto.html#howto_17)
:link(howto_18,Section_howto.html#howto_18)
:link(howto_19,Section_howto.html#howto_19)
:link(howto_20,Section_howto.html#howto_20)
:link(howto_21,Section_howto.html#howto_21)
:link(howto_22,Section_howto.html#howto_22)
:link(howto_23,Section_howto.html#howto_23)
:link(howto_24,Section_howto.html#howto_24)
:link(howto_25,Section_howto.html#howto_25)
:link(howto_26,Section_howto.html#howto_26)
:link(howto_27,Section_howto.html#howto_27)
:link(mod_1,Section_modify.html#mod_1)
:link(mod_2,Section_modify.html#mod_2)
:link(mod_3,Section_modify.html#mod_3)
:link(mod_4,Section_modify.html#mod_4)
:link(mod_5,Section_modify.html#mod_5)
:link(mod_6,Section_modify.html#mod_6)
:link(mod_7,Section_modify.html#mod_7)
:link(mod_8,Section_modify.html#mod_8)
:link(mod_9,Section_modify.html#mod_9)
:link(mod_10,Section_modify.html#mod_10)
:link(mod_11,Section_modify.html#mod_11)
:link(mod_12,Section_modify.html#mod_12)
:link(mod_13,Section_modify.html#mod_13)
:link(mod_14,Section_modify.html#mod_14)
:link(mod_15,Section_modify.html#mod_15)
:link(py_1,Section_python.html#py_1)
:link(py_2,Section_python.html#py_2)
:link(py_3,Section_python.html#py_3)
:link(py_4,Section_python.html#py_4)
:link(py_5,Section_python.html#py_5)
:link(py_6,Section_python.html#py_6)
:link(err_1,Section_errors.html#err_1)
:link(err_2,Section_errors.html#err_2)
:link(err_3,Section_errors.html#err_3)
:link(hist_1,Section_history.html#hist_1)
:link(hist_2,Section_history.html#hist_2)
diff --git a/potentials/WL.meam b/potentials/WL.meam
new file mode 100644
index 000000000..43eecb422
--- /dev/null
+++ b/potentials/WL.meam
@@ -0,0 +1,13 @@
+# DATE: 2017-01-25 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov, CITATION: Lee, Baskes, Kim, Cho. Phys. Rev. B, 64, 184102 (2001)
+rc = 3.8
+delr = 0.1
+augt1 = 0
+erose_form = 2
+zbl(1,1) = 0
+nn2(1,1) = 1
+Ec(1,1) = 8.66
+re(1,1) = 2.74
+attrac(1,1) = 0
+repuls(1,1) = 0
+Cmin(1,1,1) = 0.49
+Cmax(1,1,1) = 2.8
diff --git a/src/angle_hybrid.cpp b/src/angle_hybrid.cpp
index a477d7f8f..9b3af1856 100644
--- a/src/angle_hybrid.cpp
+++ b/src/angle_hybrid.cpp
@@ -1,376 +1,376 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include
#include
#include
#include "angle_hybrid.h"
#include "atom.h"
#include "neighbor.h"
#include "domain.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define EXTRA 1000
/* ---------------------------------------------------------------------- */
AngleHybrid::AngleHybrid(LAMMPS *lmp) : Angle(lmp)
{
writedata = 0;
nstyles = 0;
}
/* ---------------------------------------------------------------------- */
AngleHybrid::~AngleHybrid()
{
if (nstyles) {
for (int i = 0; i < nstyles; i++) delete styles[i];
delete [] styles;
for (int i = 0; i < nstyles; i++) delete [] keywords[i];
delete [] keywords;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(map);
delete [] nanglelist;
delete [] maxangle;
for (int i = 0; i < nstyles; i++)
memory->destroy(anglelist[i]);
delete [] anglelist;
}
}
/* ---------------------------------------------------------------------- */
void AngleHybrid::compute(int eflag, int vflag)
{
int i,j,m,n;
// save ptrs to original anglelist
int nanglelist_orig = neighbor->nanglelist;
int **anglelist_orig = neighbor->anglelist;
// if this is re-neighbor step, create sub-style anglelists
// nanglelist[] = length of each sub-style list
// realloc sub-style anglelist if necessary
// load sub-style anglelist with 4 values from original anglelist
if (neighbor->ago == 0) {
for (m = 0; m < nstyles; m++) nanglelist[m] = 0;
for (i = 0; i < nanglelist_orig; i++) {
m = map[anglelist_orig[i][3]];
if (m >= 0) nanglelist[m]++;
}
for (m = 0; m < nstyles; m++) {
if (nanglelist[m] > maxangle[m]) {
memory->destroy(anglelist[m]);
maxangle[m] = nanglelist[m] + EXTRA;
memory->create(anglelist[m],maxangle[m],4,"angle_hybrid:anglelist");
}
nanglelist[m] = 0;
}
for (i = 0; i < nanglelist_orig; i++) {
m = map[anglelist_orig[i][3]];
if (m < 0) continue;
n = nanglelist[m];
anglelist[m][n][0] = anglelist_orig[i][0];
anglelist[m][n][1] = anglelist_orig[i][1];
anglelist[m][n][2] = anglelist_orig[i][2];
anglelist[m][n][3] = anglelist_orig[i][3];
nanglelist[m]++;
}
}
// call each sub-style's compute function
// set neighbor->anglelist to sub-style anglelist before call
// accumulate sub-style global/peratom energy/virial in hybrid
if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = 0;
+ else evflag = eflag_global = vflag_global = eflag_atom = vflag_atom = 0;
for (m = 0; m < nstyles; m++) {
neighbor->nanglelist = nanglelist[m];
neighbor->anglelist = anglelist[m];
styles[m]->compute(eflag,vflag);
if (eflag_global) energy += styles[m]->energy;
if (vflag_global)
for (n = 0; n < 6; n++) virial[n] += styles[m]->virial[n];
if (eflag_atom) {
n = atom->nlocal;
if (force->newton_bond) n += atom->nghost;
double *eatom_substyle = styles[m]->eatom;
for (i = 0; i < n; i++) eatom[i] += eatom_substyle[i];
}
if (vflag_atom) {
n = atom->nlocal;
if (force->newton_bond) n += atom->nghost;
double **vatom_substyle = styles[m]->vatom;
for (i = 0; i < n; i++)
for (j = 0; j < 6; j++)
vatom[i][j] += vatom_substyle[i][j];
}
}
// restore ptrs to original anglelist
neighbor->nanglelist = nanglelist_orig;
neighbor->anglelist = anglelist_orig;
}
/* ---------------------------------------------------------------------- */
void AngleHybrid::allocate()
{
allocated = 1;
int n = atom->nangletypes;
memory->create(map,n+1,"angle:map");
memory->create(setflag,n+1,"angle:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
nanglelist = new int[nstyles];
maxangle = new int[nstyles];
anglelist = new int**[nstyles];
for (int m = 0; m < nstyles; m++) maxangle[m] = 0;
for (int m = 0; m < nstyles; m++) anglelist[m] = NULL;
}
/* ----------------------------------------------------------------------
create one angle style for each arg in list
------------------------------------------------------------------------- */
void AngleHybrid::settings(int narg, char **arg)
{
int i,m,istyle;
if (narg < 1) error->all(FLERR,"Illegal angle_style command");
// delete old lists, since cannot just change settings
if (nstyles) {
for (int i = 0; i < nstyles; i++) delete styles[i];
delete [] styles;
for (int i = 0; i < nstyles; i++) delete [] keywords[i];
delete [] keywords;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(map);
delete [] nanglelist;
delete [] maxangle;
for (int i = 0; i < nstyles; i++)
memory->destroy(anglelist[i]);
delete [] anglelist;
}
allocated = 0;
// count sub-styles by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric word
// need a better way to skip these exceptions
nstyles = 0;
i = 0;
while (i < narg) {
if (strcmp(arg[i],"table") == 0) i++;
i++;
while (i < narg && !isalpha(arg[i][0])) i++;
nstyles++;
}
// allocate list of sub-styles
styles = new Angle*[nstyles];
keywords = new char*[nstyles];
// allocate each sub-style and call its settings() with subset of args
// allocate uses suffix, but don't store suffix version in keywords,
// else syntax in coeff() will not match
// define subset of args for a sub-style by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric
// need a better way to skip these exceptions
int dummy;
nstyles = 0;
i = 0;
while (i < narg) {
for (m = 0; m < nstyles; m++)
if (strcmp(arg[i],keywords[m]) == 0)
error->all(FLERR,"Angle style hybrid cannot use "
"same angle style twice");
if (strcmp(arg[i],"hybrid") == 0)
error->all(FLERR,"Angle style hybrid cannot have hybrid as an argument");
if (strcmp(arg[i],"none") == 0)
error->all(FLERR,"Angle style hybrid cannot have none as an argument");
styles[nstyles] = force->new_angle(arg[i],1,dummy);
force->store_style(keywords[nstyles],arg[i],0);
istyle = i;
if (strcmp(arg[i],"table") == 0) i++;
i++;
while (i < narg && !isalpha(arg[i][0])) i++;
styles[nstyles]->settings(i-istyle-1,&arg[istyle+1]);
nstyles++;
}
}
/* ----------------------------------------------------------------------
set coeffs for one type
---------------------------------------------------------------------- */
void AngleHybrid::coeff(int narg, char **arg)
{
if (!allocated) allocate();
int ilo,ihi;
force->bounds(FLERR,arg[0],atom->nangletypes,ilo,ihi);
// 2nd arg = angle sub-style name
// allow for "none" or "skip" as valid sub-style name
int m;
for (m = 0; m < nstyles; m++)
if (strcmp(arg[1],keywords[m]) == 0) break;
int none = 0;
int skip = 0;
if (m == nstyles) {
if (strcmp(arg[1],"none") == 0) none = 1;
else if (strcmp(arg[1],"skip") == 0) none = skip = 1;
else if (strcmp(arg[1],"ba") == 0)
error->all(FLERR,"BondAngle coeff for hybrid angle has invalid format");
else if (strcmp(arg[1],"bb") == 0)
error->all(FLERR,"BondBond coeff for hybrid angle has invalid format");
else error->all(FLERR,"Angle coeff for hybrid has invalid style");
}
// move 1st arg to 2nd arg
// just copy ptrs, since arg[] points into original input line
arg[1] = arg[0];
// invoke sub-style coeff() starting with 1st arg
if (!none) styles[m]->coeff(narg-1,&arg[1]);
// set setflag and which type maps to which sub-style
// if sub-style is skip: auxiliary class2 setting in data file so ignore
// if sub-style is none: set hybrid setflag, wipe out map
for (int i = ilo; i <= ihi; i++) {
if (skip) continue;
else if (none) {
setflag[i] = 1;
map[i] = -1;
} else {
setflag[i] = styles[m]->setflag[i];
map[i] = m;
}
}
}
/* ----------------------------------------------------------------------
run angle style specific initialization
------------------------------------------------------------------------- */
void AngleHybrid::init_style()
{
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();
}
/* ----------------------------------------------------------------------
return an equilbrium angle length
------------------------------------------------------------------------- */
double AngleHybrid::equilibrium_angle(int i)
{
if (map[i] < 0)
error->one(FLERR,"Invoked angle equil angle on angle style none");
return styles[map[i]]->equilibrium_angle(i);
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void AngleHybrid::write_restart(FILE *fp)
{
fwrite(&nstyles,sizeof(int),1,fp);
int n;
for (int m = 0; m < nstyles; m++) {
n = strlen(keywords[m]) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(keywords[m],sizeof(char),n,fp);
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void AngleHybrid::read_restart(FILE *fp)
{
int me = comm->me;
if (me == 0) fread(&nstyles,sizeof(int),1,fp);
MPI_Bcast(&nstyles,1,MPI_INT,0,world);
styles = new Angle*[nstyles];
keywords = new char*[nstyles];
allocate();
int n,dummy;
for (int m = 0; m < nstyles; m++) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
keywords[m] = new char[n];
if (me == 0) fread(keywords[m],sizeof(char),n,fp);
MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
styles[m] = force->new_angle(keywords[m],0,dummy);
}
}
/* ---------------------------------------------------------------------- */
double AngleHybrid::single(int type, int i1, int i2, int i3)
{
if (map[type] < 0) error->one(FLERR,"Invoked angle single on angle style none");
return styles[map[type]]->single(type,i1,i2,i3);
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double AngleHybrid::memory_usage()
{
double bytes = maxeatom * sizeof(double);
bytes += maxvatom*6 * sizeof(double);
for (int m = 0; m < nstyles; m++) bytes += maxangle[m]*4 * sizeof(int);
for (int m = 0; m < nstyles; m++)
if (styles[m]) bytes += styles[m]->memory_usage();
return bytes;
}
diff --git a/src/bond_hybrid.cpp b/src/bond_hybrid.cpp
index 1244c2822..9a16d0e1f 100644
--- a/src/bond_hybrid.cpp
+++ b/src/bond_hybrid.cpp
@@ -1,362 +1,362 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include
#include
#include
#include "bond_hybrid.h"
#include "atom.h"
#include "neighbor.h"
#include "domain.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define EXTRA 1000
/* ---------------------------------------------------------------------- */
BondHybrid::BondHybrid(LAMMPS *lmp) : Bond(lmp)
{
writedata = 0;
nstyles = 0;
}
/* ---------------------------------------------------------------------- */
BondHybrid::~BondHybrid()
{
if (nstyles) {
for (int i = 0; i < nstyles; i++) delete styles[i];
delete [] styles;
for (int i = 0; i < nstyles; i++) delete [] keywords[i];
delete [] keywords;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(map);
delete [] nbondlist;
delete [] maxbond;
for (int i = 0; i < nstyles; i++)
memory->destroy(bondlist[i]);
delete [] bondlist;
}
}
/* ---------------------------------------------------------------------- */
void BondHybrid::compute(int eflag, int vflag)
{
int i,j,m,n;
// save ptrs to original bondlist
int nbondlist_orig = neighbor->nbondlist;
int **bondlist_orig = neighbor->bondlist;
// if this is re-neighbor step, create sub-style bondlists
// nbondlist[] = length of each sub-style list
// realloc sub-style bondlist if necessary
// load sub-style bondlist with 3 values from original bondlist
if (neighbor->ago == 0) {
for (m = 0; m < nstyles; m++) nbondlist[m] = 0;
for (i = 0; i < nbondlist_orig; i++) {
m = map[bondlist_orig[i][2]];
if (m >= 0) nbondlist[m]++;
}
for (m = 0; m < nstyles; m++) {
if (nbondlist[m] > maxbond[m]) {
memory->destroy(bondlist[m]);
maxbond[m] = nbondlist[m] + EXTRA;
memory->create(bondlist[m],maxbond[m],3,"bond_hybrid:bondlist");
}
nbondlist[m] = 0;
}
for (i = 0; i < nbondlist_orig; i++) {
m = map[bondlist_orig[i][2]];
if (m < 0) continue;
n = nbondlist[m];
bondlist[m][n][0] = bondlist_orig[i][0];
bondlist[m][n][1] = bondlist_orig[i][1];
bondlist[m][n][2] = bondlist_orig[i][2];
nbondlist[m]++;
}
}
// call each sub-style's compute function
// set neighbor->bondlist to sub-style bondlist before call
// accumulate sub-style global/peratom energy/virial in hybrid
if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = 0;
+ else evflag = eflag_global = vflag_global = eflag_atom = vflag_atom = 0;
for (m = 0; m < nstyles; m++) {
neighbor->nbondlist = nbondlist[m];
neighbor->bondlist = bondlist[m];
styles[m]->compute(eflag,vflag);
if (eflag_global) energy += styles[m]->energy;
if (vflag_global)
for (n = 0; n < 6; n++) virial[n] += styles[m]->virial[n];
if (eflag_atom) {
n = atom->nlocal;
if (force->newton_bond) n += atom->nghost;
double *eatom_substyle = styles[m]->eatom;
for (i = 0; i < n; i++) eatom[i] += eatom_substyle[i];
}
if (vflag_atom) {
n = atom->nlocal;
if (force->newton_bond) n += atom->nghost;
double **vatom_substyle = styles[m]->vatom;
for (i = 0; i < n; i++)
for (j = 0; j < 6; j++)
vatom[i][j] += vatom_substyle[i][j];
}
}
// restore ptrs to original bondlist
neighbor->nbondlist = nbondlist_orig;
neighbor->bondlist = bondlist_orig;
}
/* ---------------------------------------------------------------------- */
void BondHybrid::allocate()
{
allocated = 1;
int n = atom->nbondtypes;
memory->create(map,n+1,"bond:map");
memory->create(setflag,n+1,"bond:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
nbondlist = new int[nstyles];
maxbond = new int[nstyles];
bondlist = new int**[nstyles];
for (int m = 0; m < nstyles; m++) maxbond[m] = 0;
for (int m = 0; m < nstyles; m++) bondlist[m] = NULL;
}
/* ----------------------------------------------------------------------
create one bond style for each arg in list
------------------------------------------------------------------------- */
void BondHybrid::settings(int narg, char **arg)
{
int i,m,istyle;
if (narg < 1) error->all(FLERR,"Illegal bond_style command");
// delete old lists, since cannot just change settings
if (nstyles) {
for (int i = 0; i < nstyles; i++) delete styles[i];
delete [] styles;
for (int i = 0; i < nstyles; i++) delete [] keywords[i];
delete [] keywords;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(map);
delete [] nbondlist;
delete [] maxbond;
for (int i = 0; i < nstyles; i++)
memory->destroy(bondlist[i]);
delete [] bondlist;
}
allocated = 0;
// count sub-styles by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric word
// need a better way to skip these exceptions
nstyles = 0;
i = 0;
while (i < narg) {
if (strcmp(arg[i],"table") == 0) i++;
i++;
while (i < narg && !isalpha(arg[i][0])) i++;
nstyles++;
}
// allocate list of sub-styles
styles = new Bond*[nstyles];
keywords = new char*[nstyles];
// allocate each sub-style and call its settings() with subset of args
// allocate uses suffix, but don't store suffix version in keywords,
// else syntax in coeff() will not match
// define subset of args for a sub-style by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric
// need a better way to skip these exceptions
int dummy;
nstyles = 0;
i = 0;
while (i < narg) {
for (m = 0; m < nstyles; m++)
if (strcmp(arg[i],keywords[m]) == 0)
error->all(FLERR,"Bond style hybrid cannot use same bond style twice");
if (strcmp(arg[i],"hybrid") == 0)
error->all(FLERR,"Bond style hybrid cannot have hybrid as an argument");
if (strcmp(arg[i],"none") == 0)
error->all(FLERR,"Bond style hybrid cannot have none as an argument");
styles[nstyles] = force->new_bond(arg[i],1,dummy);
force->store_style(keywords[nstyles],arg[i],0);
istyle = i;
if (strcmp(arg[i],"table") == 0) i++;
i++;
while (i < narg && !isalpha(arg[i][0])) i++;
styles[nstyles]->settings(i-istyle-1,&arg[istyle+1]);
nstyles++;
}
}
/* ----------------------------------------------------------------------
set coeffs for one type
---------------------------------------------------------------------- */
void BondHybrid::coeff(int narg, char **arg)
{
if (!allocated) allocate();
int ilo,ihi;
force->bounds(FLERR,arg[0],atom->nbondtypes,ilo,ihi);
// 2nd arg = bond sub-style name
// allow for "none" as valid sub-style name
int m;
for (m = 0; m < nstyles; m++)
if (strcmp(arg[1],keywords[m]) == 0) break;
int none = 0;
if (m == nstyles) {
if (strcmp(arg[1],"none") == 0) none = 1;
else error->all(FLERR,"Bond coeff for hybrid has invalid style");
}
// move 1st arg to 2nd arg
// just copy ptrs, since arg[] points into original input line
arg[1] = arg[0];
// invoke sub-style coeff() starting with 1st arg
if (!none) styles[m]->coeff(narg-1,&arg[1]);
// set setflag and which type maps to which sub-style
// if sub-style is none: set hybrid setflag, wipe out map
for (int i = ilo; i <= ihi; i++) {
setflag[i] = 1;
if (none) map[i] = -1;
else map[i] = m;
}
}
/* ---------------------------------------------------------------------- */
void BondHybrid::init_style()
{
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();
}
/* ----------------------------------------------------------------------
return an equilbrium bond length
------------------------------------------------------------------------- */
double BondHybrid::equilibrium_distance(int i)
{
if (map[i] < 0)
error->one(FLERR,"Invoked bond equil distance on bond style none");
return styles[map[i]]->equilibrium_distance(i);
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void BondHybrid::write_restart(FILE *fp)
{
fwrite(&nstyles,sizeof(int),1,fp);
int n;
for (int m = 0; m < nstyles; m++) {
n = strlen(keywords[m]) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(keywords[m],sizeof(char),n,fp);
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void BondHybrid::read_restart(FILE *fp)
{
int me = comm->me;
if (me == 0) fread(&nstyles,sizeof(int),1,fp);
MPI_Bcast(&nstyles,1,MPI_INT,0,world);
styles = new Bond*[nstyles];
keywords = new char*[nstyles];
allocate();
int n,dummy;
for (int m = 0; m < nstyles; m++) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
keywords[m] = new char[n];
if (me == 0) fread(keywords[m],sizeof(char),n,fp);
MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
styles[m] = force->new_bond(keywords[m],0,dummy);
}
}
/* ---------------------------------------------------------------------- */
double BondHybrid::single(int type, double rsq, int i, int j,
double &fforce)
{
if (map[type] < 0) error->one(FLERR,"Invoked bond single on bond style none");
return styles[map[type]]->single(type,rsq,i,j,fforce);
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double BondHybrid::memory_usage()
{
double bytes = maxeatom * sizeof(double);
bytes += maxvatom*6 * sizeof(double);
for (int m = 0; m < nstyles; m++) bytes += maxbond[m]*3 * sizeof(int);
for (int m = 0; m < nstyles; m++)
if (styles[m]) bytes += styles[m]->memory_usage();
return bytes;
}
diff --git a/src/dihedral_hybrid.cpp b/src/dihedral_hybrid.cpp
index 0ae396b88..372a858d0 100644
--- a/src/dihedral_hybrid.cpp
+++ b/src/dihedral_hybrid.cpp
@@ -1,353 +1,353 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include
#include
#include
#include "dihedral_hybrid.h"
#include "atom.h"
#include "neighbor.h"
#include "domain.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define EXTRA 1000
/* ---------------------------------------------------------------------- */
DihedralHybrid::DihedralHybrid(LAMMPS *lmp) : Dihedral(lmp)
{
nstyles = 0;
}
/* ---------------------------------------------------------------------- */
DihedralHybrid::~DihedralHybrid()
{
if (nstyles) {
for (int i = 0; i < nstyles; i++) delete styles[i];
delete [] styles;
for (int i = 0; i < nstyles; i++) delete [] keywords[i];
delete [] keywords;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(map);
delete [] ndihedrallist;
delete [] maxdihedral;
for (int i = 0; i < nstyles; i++)
memory->destroy(dihedrallist[i]);
delete [] dihedrallist;
}
}
/* ---------------------------------------------------------------------- */
void DihedralHybrid::compute(int eflag, int vflag)
{
int i,j,m,n;
// save ptrs to original dihedrallist
int ndihedrallist_orig = neighbor->ndihedrallist;
int **dihedrallist_orig = neighbor->dihedrallist;
// if this is re-neighbor step, create sub-style dihedrallists
// ndihedrallist[] = length of each sub-style list
// realloc sub-style dihedrallist if necessary
// load sub-style dihedrallist with 5 values from original dihedrallist
if (neighbor->ago == 0) {
for (m = 0; m < nstyles; m++) ndihedrallist[m] = 0;
for (i = 0; i < ndihedrallist_orig; i++) {
m = map[dihedrallist_orig[i][4]];
if (m >= 0) ndihedrallist[m]++;
}
for (m = 0; m < nstyles; m++) {
if (ndihedrallist[m] > maxdihedral[m]) {
memory->destroy(dihedrallist[m]);
maxdihedral[m] = ndihedrallist[m] + EXTRA;
memory->create(dihedrallist[m],maxdihedral[m],5,
"dihedral_hybrid:dihedrallist");
}
ndihedrallist[m] = 0;
}
for (i = 0; i < ndihedrallist_orig; i++) {
m = map[dihedrallist_orig[i][4]];
if (m < 0) continue;
n = ndihedrallist[m];
dihedrallist[m][n][0] = dihedrallist_orig[i][0];
dihedrallist[m][n][1] = dihedrallist_orig[i][1];
dihedrallist[m][n][2] = dihedrallist_orig[i][2];
dihedrallist[m][n][3] = dihedrallist_orig[i][3];
dihedrallist[m][n][4] = dihedrallist_orig[i][4];
ndihedrallist[m]++;
}
}
// call each sub-style's compute function
// set neighbor->dihedrallist to sub-style dihedrallist before call
// accumulate sub-style global/peratom energy/virial in hybrid
if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = 0;
+ else evflag = eflag_global = vflag_global = eflag_atom = vflag_atom = 0;
for (m = 0; m < nstyles; m++) {
neighbor->ndihedrallist = ndihedrallist[m];
neighbor->dihedrallist = dihedrallist[m];
styles[m]->compute(eflag,vflag);
if (eflag_global) energy += styles[m]->energy;
if (vflag_global)
for (n = 0; n < 6; n++) virial[n] += styles[m]->virial[n];
if (eflag_atom) {
n = atom->nlocal;
if (force->newton_bond) n += atom->nghost;
double *eatom_substyle = styles[m]->eatom;
for (i = 0; i < n; i++) eatom[i] += eatom_substyle[i];
}
if (vflag_atom) {
n = atom->nlocal;
if (force->newton_bond) n += atom->nghost;
double **vatom_substyle = styles[m]->vatom;
for (i = 0; i < n; i++)
for (j = 0; j < 6; j++)
vatom[i][j] += vatom_substyle[i][j];
}
}
// restore ptrs to original dihedrallist
neighbor->ndihedrallist = ndihedrallist_orig;
neighbor->dihedrallist = dihedrallist_orig;
}
/* ---------------------------------------------------------------------- */
void DihedralHybrid::allocate()
{
allocated = 1;
int n = atom->ndihedraltypes;
memory->create(map,n+1,"dihedral:map");
memory->create(setflag,n+1,"dihedral:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
ndihedrallist = new int[nstyles];
maxdihedral = new int[nstyles];
dihedrallist = new int**[nstyles];
for (int m = 0; m < nstyles; m++) maxdihedral[m] = 0;
for (int m = 0; m < nstyles; m++) dihedrallist[m] = NULL;
}
/* ----------------------------------------------------------------------
create one dihedral style for each arg in list
------------------------------------------------------------------------- */
void DihedralHybrid::settings(int narg, char **arg)
{
int i,m,istyle;
if (narg < 1) error->all(FLERR,"Illegal dihedral_style command");
// delete old lists, since cannot just change settings
if (nstyles) {
for (int i = 0; i < nstyles; i++) delete styles[i];
delete [] styles;
for (int i = 0; i < nstyles; i++) delete [] keywords[i];
delete [] keywords;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(map);
delete [] ndihedrallist;
delete [] maxdihedral;
for (int i = 0; i < nstyles; i++)
memory->destroy(dihedrallist[i]);
delete [] dihedrallist;
}
allocated = 0;
// count sub-styles by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric word
// need a better way to skip these exceptions
nstyles = 0;
i = 0;
while (i < narg) {
if (strcmp(arg[i],"table") == 0) i++;
i++;
while (i < narg && !isalpha(arg[i][0])) i++;
nstyles++;
}
// allocate list of sub-styles
styles = new Dihedral*[nstyles];
keywords = new char*[nstyles];
// allocate each sub-style and call its settings() with subset of args
// allocate uses suffix, but don't store suffix version in keywords,
// else syntax in coeff() will not match
// define subset of args for a sub-style by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric
// need a better way to skip these exceptions
int dummy;
nstyles = 0;
i = 0;
while (i < narg) {
for (m = 0; m < nstyles; m++)
if (strcmp(arg[i],keywords[m]) == 0)
error->all(FLERR,"Dihedral style hybrid cannot use "
"same dihedral style twice");
if (strcmp(arg[i],"hybrid") == 0)
error->all(FLERR,
"Dihedral style hybrid cannot have hybrid as an argument");
if (strcmp(arg[i],"none") == 0)
error->all(FLERR,"Dihedral style hybrid cannot have none as an argument");
styles[nstyles] = force->new_dihedral(arg[i],1,dummy);
force->store_style(keywords[nstyles],arg[i],0);
istyle = i;
if (strcmp(arg[i],"table") == 0) i++;
i++;
while (i < narg && !isalpha(arg[i][0])) i++;
styles[nstyles]->settings(i-istyle-1,&arg[istyle+1]);
nstyles++;
}
}
/* ----------------------------------------------------------------------
set coeffs for one type
---------------------------------------------------------------------- */
void DihedralHybrid::coeff(int narg, char **arg)
{
if (!allocated) allocate();
int ilo,ihi;
force->bounds(FLERR,arg[0],atom->ndihedraltypes,ilo,ihi);
// 2nd arg = dihedral sub-style name
// allow for "none" or "skip" as valid sub-style name
int m;
for (m = 0; m < nstyles; m++)
if (strcmp(arg[1],keywords[m]) == 0) break;
int none = 0;
int skip = 0;
if (m == nstyles) {
if (strcmp(arg[1],"none") == 0) none = 1;
else if (strcmp(arg[1],"skip") == 0) none = skip = 1;
else error->all(FLERR,"Dihedral coeff for hybrid has invalid style");
}
// move 1st arg to 2nd arg
// just copy ptrs, since arg[] points into original input line
arg[1] = arg[0];
// invoke sub-style coeff() starting with 1st arg
if (!none) styles[m]->coeff(narg-1,&arg[1]);
// set setflag and which type maps to which sub-style
// if sub-style is skip: auxiliary class2 setting in data file so ignore
// if sub-style is none and not skip: set hybrid setflag, wipe out map
for (int i = ilo; i <= ihi; i++) {
if (skip) continue;
else if (none) {
setflag[i] = 1;
map[i] = -1;
} else {
setflag[i] = styles[m]->setflag[i];
map[i] = m;
}
}
}
/* ---------------------------------------------------------------------- */
void DihedralHybrid::init_style()
{
for (int m = 0; m < nstyles; m++)
if (styles[m]) styles[m]->init_style();
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void DihedralHybrid::write_restart(FILE *fp)
{
fwrite(&nstyles,sizeof(int),1,fp);
int n;
for (int m = 0; m < nstyles; m++) {
n = strlen(keywords[m]) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(keywords[m],sizeof(char),n,fp);
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void DihedralHybrid::read_restart(FILE *fp)
{
int me = comm->me;
if (me == 0) fread(&nstyles,sizeof(int),1,fp);
MPI_Bcast(&nstyles,1,MPI_INT,0,world);
styles = new Dihedral*[nstyles];
keywords = new char*[nstyles];
allocate();
int n,dummy;
for (int m = 0; m < nstyles; m++) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
keywords[m] = new char[n];
if (me == 0) fread(keywords[m],sizeof(char),n,fp);
MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
styles[m] = force->new_dihedral(keywords[m],0,dummy);
}
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double DihedralHybrid::memory_usage()
{
double bytes = maxeatom * sizeof(double);
bytes += maxvatom*6 * sizeof(double);
for (int m = 0; m < nstyles; m++) bytes += maxdihedral[m]*5 * sizeof(int);
for (int m = 0; m < nstyles; m++)
if (styles[m]) bytes += styles[m]->memory_usage();
return bytes;
}
diff --git a/src/improper_hybrid.cpp b/src/improper_hybrid.cpp
index 78ee69569..abaaae02d 100644
--- a/src/improper_hybrid.cpp
+++ b/src/improper_hybrid.cpp
@@ -1,350 +1,350 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include
#include
#include
#include "improper_hybrid.h"
#include "atom.h"
#include "neighbor.h"
#include "domain.h"
#include "comm.h"
#include "force.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
#define EXTRA 1000
/* ---------------------------------------------------------------------- */
ImproperHybrid::ImproperHybrid(LAMMPS *lmp) : Improper(lmp)
{
nstyles = 0;
}
/* ---------------------------------------------------------------------- */
ImproperHybrid::~ImproperHybrid()
{
if (nstyles) {
for (int i = 0; i < nstyles; i++) delete styles[i];
delete [] styles;
for (int i = 0; i < nstyles; i++) delete [] keywords[i];
delete [] keywords;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(map);
delete [] nimproperlist;
delete [] maximproper;
for (int i = 0; i < nstyles; i++)
memory->destroy(improperlist[i]);
delete [] improperlist;
}
}
/* ---------------------------------------------------------------------- */
void ImproperHybrid::compute(int eflag, int vflag)
{
int i,j,m,n;
// save ptrs to original improperlist
int nimproperlist_orig = neighbor->nimproperlist;
int **improperlist_orig = neighbor->improperlist;
// if this is re-neighbor step, create sub-style improperlists
// nimproperlist[] = length of each sub-style list
// realloc sub-style improperlist if necessary
// load sub-style improperlist with 5 values from original improperlist
if (neighbor->ago == 0) {
for (m = 0; m < nstyles; m++) nimproperlist[m] = 0;
for (i = 0; i < nimproperlist_orig; i++) {
m = map[improperlist_orig[i][4]];
if (m >= 0) nimproperlist[m]++;
}
for (m = 0; m < nstyles; m++) {
if (nimproperlist[m] > maximproper[m]) {
memory->destroy(improperlist[m]);
maximproper[m] = nimproperlist[m] + EXTRA;
memory->create(improperlist[m],maximproper[m],5,
"improper_hybrid:improperlist");
}
nimproperlist[m] = 0;
}
for (i = 0; i < nimproperlist_orig; i++) {
m = map[improperlist_orig[i][4]];
if (m < 0) continue;
n = nimproperlist[m];
improperlist[m][n][0] = improperlist_orig[i][0];
improperlist[m][n][1] = improperlist_orig[i][1];
improperlist[m][n][2] = improperlist_orig[i][2];
improperlist[m][n][3] = improperlist_orig[i][3];
improperlist[m][n][4] = improperlist_orig[i][4];
nimproperlist[m]++;
}
}
// call each sub-style's compute function
// set neighbor->improperlist to sub-style improperlist before call
// accumulate sub-style global/peratom energy/virial in hybrid
if (eflag || vflag) ev_setup(eflag,vflag);
- else evflag = 0;
+ else evflag = eflag_global = vflag_global = eflag_atom = vflag_atom = 0;
for (m = 0; m < nstyles; m++) {
neighbor->nimproperlist = nimproperlist[m];
neighbor->improperlist = improperlist[m];
styles[m]->compute(eflag,vflag);
if (eflag_global) energy += styles[m]->energy;
if (vflag_global)
for (n = 0; n < 6; n++) virial[n] += styles[m]->virial[n];
if (eflag_atom) {
n = atom->nlocal;
if (force->newton_bond) n += atom->nghost;
double *eatom_substyle = styles[m]->eatom;
for (i = 0; i < n; i++) eatom[i] += eatom_substyle[i];
}
if (vflag_atom) {
n = atom->nlocal;
if (force->newton_bond) n += atom->nghost;
double **vatom_substyle = styles[m]->vatom;
for (i = 0; i < n; i++)
for (j = 0; j < 6; j++)
vatom[i][j] += vatom_substyle[i][j];
}
}
// restore ptrs to original improperlist
neighbor->nimproperlist = nimproperlist_orig;
neighbor->improperlist = improperlist_orig;
}
/* ---------------------------------------------------------------------- */
void ImproperHybrid::allocate()
{
allocated = 1;
int n = atom->nimpropertypes;
memory->create(map,n+1,"improper:map");
memory->create(setflag,n+1,"improper:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
nimproperlist = new int[nstyles];
maximproper = new int[nstyles];
improperlist = new int**[nstyles];
for (int m = 0; m < nstyles; m++) maximproper[m] = 0;
for (int m = 0; m < nstyles; m++) improperlist[m] = NULL;
}
/* ---------------------------------------------------------------------- */
void ImproperHybrid::init_style()
{
for (int i = 0; i < nstyles; i++)
styles[i]->init_style();
}
/* ----------------------------------------------------------------------
create one improper style for each arg in list
------------------------------------------------------------------------- */
void ImproperHybrid::settings(int narg, char **arg)
{
int i,m,istyle;
if (narg < 1) error->all(FLERR,"Illegal improper_style command");
// delete old lists, since cannot just change settings
if (nstyles) {
for (int i = 0; i < nstyles; i++) delete styles[i];
delete [] styles;
for (int i = 0; i < nstyles; i++) delete [] keywords[i];
delete [] keywords;
}
if (allocated) {
memory->destroy(setflag);
memory->destroy(map);
delete [] nimproperlist;
delete [] maximproper;
for (int i = 0; i < nstyles; i++)
memory->destroy(improperlist[i]);
delete [] improperlist;
}
allocated = 0;
// count sub-styles by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric word
// need a better way to skip these exceptions
nstyles = 0;
i = 0;
while (i < narg) {
if (strcmp(arg[i],"table") == 0) i++;
i++;
while (i < narg && !isalpha(arg[i][0])) i++;
nstyles++;
}
// allocate list of sub-styles
styles = new Improper*[nstyles];
keywords = new char*[nstyles];
// allocate each sub-style and call its settings() with subset of args
// allocate uses suffix, but don't store suffix version in keywords,
// else syntax in coeff() will not match
// define subset of args for a sub-style by skipping numeric args
// one exception is 1st arg of style "table", which is non-numeric
// need a better way to skip these exceptions
int dummy;
nstyles = 0;
i = 0;
while (i < narg) {
for (m = 0; m < nstyles; m++)
if (strcmp(arg[i],keywords[m]) == 0)
error->all(FLERR,"Improper style hybrid cannot use "
"same improper style twice");
if (strcmp(arg[i],"hybrid") == 0)
error->all(FLERR,
"Improper style hybrid cannot have hybrid as an argument");
if (strcmp(arg[i],"none") == 0)
error->all(FLERR,"Improper style hybrid cannot have none as an argument");
styles[nstyles] = force->new_improper(arg[i],1,dummy);
force->store_style(keywords[nstyles],arg[i],0);
istyle = i;
if (strcmp(arg[i],"table") == 0) i++;
i++;
while (i < narg && !isalpha(arg[i][0])) i++;
styles[nstyles]->settings(i-istyle-1,&arg[istyle+1]);
nstyles++;
}
}
/* ----------------------------------------------------------------------
set coeffs for one type
---------------------------------------------------------------------- */
void ImproperHybrid::coeff(int narg, char **arg)
{
if (!allocated) allocate();
int ilo,ihi;
force->bounds(FLERR,arg[0],atom->nimpropertypes,ilo,ihi);
// 2nd arg = improper sub-style name
// allow for "none" as valid sub-style name
int m;
for (m = 0; m < nstyles; m++)
if (strcmp(arg[1],keywords[m]) == 0) break;
int none = 0;
if (m == nstyles) {
if (strcmp(arg[1],"none") == 0) none = 1;
else error->all(FLERR,"Improper coeff for hybrid has invalid style");
}
// move 1st arg to 2nd arg
// just copy ptrs, since arg[] points into original input line
arg[1] = arg[0];
// invoke sub-style coeff() starting with 1st arg
if (!none) styles[m]->coeff(narg-1,&arg[1]);
// set setflag and which type maps to which sub-style
// if sub-style is none: set hybrid setflag, wipe out map
for (int i = ilo; i <= ihi; i++) {
if (none) {
setflag[i] = 1;
map[i] = -1;
} else {
setflag[i] = styles[m]->setflag[i];
map[i] = m;
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void ImproperHybrid::write_restart(FILE *fp)
{
fwrite(&nstyles,sizeof(int),1,fp);
int n;
for (int m = 0; m < nstyles; m++) {
n = strlen(keywords[m]) + 1;
fwrite(&n,sizeof(int),1,fp);
fwrite(keywords[m],sizeof(char),n,fp);
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void ImproperHybrid::read_restart(FILE *fp)
{
int me = comm->me;
if (me == 0) fread(&nstyles,sizeof(int),1,fp);
MPI_Bcast(&nstyles,1,MPI_INT,0,world);
styles = new Improper*[nstyles];
keywords = new char*[nstyles];
allocate();
int n,dummy;
for (int m = 0; m < nstyles; m++) {
if (me == 0) fread(&n,sizeof(int),1,fp);
MPI_Bcast(&n,1,MPI_INT,0,world);
keywords[m] = new char[n];
if (me == 0) fread(keywords[m],sizeof(char),n,fp);
MPI_Bcast(keywords[m],n,MPI_CHAR,0,world);
styles[m] = force->new_improper(keywords[m],0,dummy);
}
}
/* ----------------------------------------------------------------------
memory usage
------------------------------------------------------------------------- */
double ImproperHybrid::memory_usage()
{
double bytes = maxeatom * sizeof(double);
bytes += maxvatom*6 * sizeof(double);
for (int m = 0; m < nstyles; m++) bytes += maximproper[m]*5 * sizeof(int);
for (int m = 0; m < nstyles; m++)
if (styles[m]) bytes += styles[m]->memory_usage();
return bytes;
}
diff --git a/src/neigh_request.cpp b/src/neigh_request.cpp
index 6354af4d3..4007daadc 100644
--- a/src/neigh_request.cpp
+++ b/src/neigh_request.cpp
@@ -1,240 +1,245 @@
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include "neigh_request.h"
#include "atom.h"
#include "memory.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
NeighRequest::NeighRequest(LAMMPS *lmp) : Pointers(lmp)
{
// default ID = 0
id = 0;
// class user of list: default is pair request
// only one is set to 1
pair = 1;
fix = compute = command = 0;
// kind of list: default is half neighbor list
// only one is set to 1
half = 1;
full = 0;
gran = granhistory = 0;
respainner = respamiddle = respaouter = 0;
half_from_full = 0;
full_cluster = 0;
// only set when command = 1;
command_style = NULL;
// combination of settings, mutiple can be set to 1
- // default is every reneighboring
+ // default is every reneighboring, not occasional
// default is use newton_pair setting in force
// default is no size history (when gran = 1)
// default is no one-sided sphere/surface interactions (when gran = 1)
// default is no auxiliary floating point values
// default is no neighbors of ghosts
// default is no multi-threaded neighbor list build
// default is no Kokkos neighbor list build
// default is no Shardlow Splitting Algorithm (SSA) neighbor list build
+ // default is neighbors of atoms, not bonds
occasional = 0;
newton = 0;
granonesided = 0;
dnum = 0;
ghost = 0;
omp = 0;
intel = 0;
kokkos_host = kokkos_device = 0;
ssa = 0;
+ bond = 0;
// copy/skip/derive info, default is no copy or skip
// none or only one option is set
copy = 0;
skip = 0;
iskip = NULL;
ijskip = NULL;
off2on = 0;
otherlist = -1;
}
/* ---------------------------------------------------------------------- */
NeighRequest::~NeighRequest()
{
delete [] iskip;
memory->destroy(ijskip);
}
/* ----------------------------------------------------------------------
archive request params that Neighbor may change after call to identical()
------------------------------------------------------------------------- */
void NeighRequest::archive()
{
half_original = half;
half_from_full_original = half_from_full;
copy_original = copy;
otherlist_original = otherlist;
}
/* ----------------------------------------------------------------------
compare this request to other request
identical means all params set by requestor are the same
compare to original values in other if Neighbor may have changed them
return 1 if identical, 0 if not
------------------------------------------------------------------------- */
int NeighRequest::identical(NeighRequest *other)
{
int same = 1;
// set same = 0 if old list was never processed
// use of requestor_instance and instance counter
// prevents an old fix from being unfix/refix in same memory location
// and appearing to be old, when it is really new
// only needed for classes with persistent neigh lists: Fix, Compute, Pair
if (requestor != other->requestor) same = 0;
if (requestor_instance != other->requestor_instance) same = 0;
if (id != other->id) same = 0;
if (pair != other->pair) same = 0;
if (fix != other->fix) same = 0;
if (compute != other->compute) same = 0;
if (command != other->command) same = 0;
if (half != other->half_original) same = 0;
if (full != other->full) same = 0;
if (gran != other->gran) same = 0;
if (respainner != other->respainner) same = 0;
if (respamiddle != other->respamiddle) same = 0;
if (respaouter != other->respaouter) same = 0;
if (half_from_full != other->half_from_full_original) same = 0;
if (newton != other->newton) same = 0;
if (occasional != other->occasional) same = 0;
if (granhistory != other->granhistory) same = 0;
if (granonesided != other->granonesided) same = 0;
if (dnum != other->dnum) same = 0;
if (ghost != other->ghost) same = 0;
if (omp != other->omp) same = 0;
if (intel != other->intel) same = 0;
if (kokkos_host != other->kokkos_host) same = 0;
if (kokkos_device != other->kokkos_device) same = 0;
if (ssa != other->ssa) same = 0;
+ if (bond != other->bond) same = 0;
if (copy != other->copy_original) same = 0;
if (same_skip(other) == 0) same = 0;
if (otherlist != other->otherlist_original) same = 0;
return same;
}
/* ----------------------------------------------------------------------
compare kind of this request to other request
return 1 if same, 0 if different
------------------------------------------------------------------------- */
int NeighRequest::same_kind(NeighRequest *other)
{
int same = 1;
// class caller will be same (pair), b/c called by pair hybrid styles
// kind must be the the same
if (half != other->half) same = 0;
if (full != other->full) same = 0;
if (gran != other->gran) same = 0;
if (granhistory != other->granhistory) same = 0;
if (respainner != other->respainner) same = 0;
if (respamiddle != other->respamiddle) same = 0;
if (respaouter != other->respaouter) same = 0;
if (half_from_full != other->half_from_full) same = 0;
// settings that must be the same
// other settings do not need to be the same (e.g. granonesided)
if (newton != other->newton) same = 0;
if (ghost != other->ghost) same = 0;
if (omp != other->omp) same = 0;
if (intel != other->intel) same = 0;
if (kokkos_host != other->kokkos_host) same = 0;
if (kokkos_device != other->kokkos_device) same = 0;
if (ssa != other->ssa) same = 0;
+ if (bond != other->bond) same = 0;
// copy/skip/derive info does not need to be the same
return same;
}
/* ----------------------------------------------------------------------
compare skip attributes of this request to other request
return 1 if same, 0 if different
------------------------------------------------------------------------- */
int NeighRequest::same_skip(NeighRequest *other)
{
int i,j;
int same = 1;
if (skip != other->skip) same = 0;
if (skip && other->skip) {
int ntypes = atom->ntypes;
for (i = 1; i <= ntypes; i++)
if (iskip[i] != other->iskip[i]) same = 0;
for (i = 1; i <= ntypes; i++)
for (j = 1; j <= ntypes; j++)
if (ijskip[i][j] != other->ijskip[i][j]) same = 0;
}
return same;
}
/* ----------------------------------------------------------------------
set kind and optional values in this request to those of other request
not copy/skip/derive values since this may be non-skip parent of child
------------------------------------------------------------------------- */
void NeighRequest::copy_request(NeighRequest *other)
{
half = other->half;
full = other->full;
gran = other->gran;
granhistory = other->granhistory;
respainner = other->respainner;
respamiddle = other->respamiddle;
respaouter = other->respaouter;
half_from_full = other->half_from_full;
newton = other->newton;
granonesided = other->granonesided;
dnum = other->dnum;
ghost = other->ghost;
omp = other->omp;
intel = other->intel;
kokkos_host = other->kokkos_host;
kokkos_device = other->kokkos_device;
ssa = other->ssa;
+ bond = other->bond;
}
diff --git a/src/neigh_request.h b/src/neigh_request.h
index 0b561710e..2738f4c03 100644
--- a/src/neigh_request.h
+++ b/src/neigh_request.h
@@ -1,139 +1,143 @@
/* -*- 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.
------------------------------------------------------------------------- */
#ifndef LMP_NEIGH_REQUEST_H
#define LMP_NEIGH_REQUEST_H
#include "pointers.h"
namespace LAMMPS_NS {
class NeighRequest : protected Pointers {
public:
int index; // index of which neigh request this is
void *requestor; // class that made request
int requestor_instance; // instance of that class (only Fix, Compute, Pair)
int id; // ID of request as stored by requestor
// used to track multiple requests from one class
// which class style requests the list, one flag is 1, others are 0
int pair; // set by default
int fix;
int compute;
int command;
// kind of list requested, one flag is 1, others are 0
// NOTE: should make only the first 3 settings be unique,
// allow others as add-on flags
// this will require changed flags in pair requestors
// this will lead to simpler logic in Neighbor::choose_build()
int half; // 1 if half neigh list (set by default)
int full; // 1 if full neigh list
int half_from_full; // 1 if half list computed from previous full list
int gran; // 1 if granular list
int granhistory; // 1 if history info for granular contact pairs
int respainner; // 1 if a rRESPA inner list
int respamiddle; // 1 if a rRESPA middle list
int respaouter; // 1 if a rRESPA outer list
int full_cluster; // only used by Kokkos pair styles
// command_style only set if command = 1
// allows print_pair_info() to access command name
const char *command_style;
// -----------------
// optional settings, set by caller, all are 0 by default
// -----------------
// 0 if needed every reneighboring during run
// 1 if occasionally needed by a fix, compute, etc
int occasional;
// 0 if use force::newton_pair setting
// 1 if override with pair newton on
// 2 if override with pair newton off
int newton;
// 1 if one-sided granular list for sphere/surf interactions (gran = 1)
int granonesided;
// number of auxiliary floating point values to store, 0 if none set
int dnum;
// 1 if also need neighbors of ghosts
int ghost;
// 1 if using multi-threaded neighbor list build for USER-OMP or USER-INTEL
int omp;
int intel;
// 1 if using Kokkos neighbor build
int kokkos_host;
int kokkos_device;
// 1 if using Shardlow Splitting Algorithm (SSA) neighbor list build
int ssa;
+ // 1 if bond neighbors, not atom neighbors
+
+ int bond;
+
// -----------------
// end of optional settings
// -----------------
// set by pair_hybrid and neighbor after all requests are made
// these settings do not change kind value or optional settings
int copy; // 1 if this list copied from another list
int skip; // 1 if this list skips atom types from another list
int *iskip; // iskip[i] if atoms of type I are not in list
int **ijskip; // ijskip[i][j] if pairs of type I,J are not in list
int off2on; // 1 if this is newton on list, but skips from off list
int otherlist; // index of other list to copy or skip from
// original params by requestor
// stored to compare against in identical() in case Neighbor changes them
int half_original;
int half_from_full_original;
int copy_original;
int otherlist_original;
// pointer to FSH class, set by caller
class FixShearHistory *fix_history; // fix that stores history info
// methods
NeighRequest(class LAMMPS *);
~NeighRequest();
void archive();
int identical(NeighRequest *);
int same_kind(NeighRequest *);
int same_skip(NeighRequest *);
void copy_request(NeighRequest *);
};
}
#endif
diff --git a/src/neighbor.cpp b/src/neighbor.cpp
index cbcd5d640..79958a29f 100644
--- a/src/neighbor.cpp
+++ b/src/neighbor.cpp
@@ -1,2084 +1,2086 @@
/* ----------------------------------------------------------------------
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 (triclinic and multi-neigh) : Pieter in 't Veld (SNL)
------------------------------------------------------------------------- */
#include
#include
#include
#include
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "style_nbin.h"
#include "style_nstencil.h"
#include "style_npair.h"
#include "style_ntopo.h"
#include "atom.h"
#include "atom_vec.h"
#include "comm.h"
#include "force.h"
#include "pair.h"
#include "domain.h"
#include "group.h"
#include "modify.h"
#include "fix.h"
#include "compute.h"
#include "update.h"
#include "respa.h"
#include "output.h"
#include "citeme.h"
#include "memory.h"
#include "error.h"
#include