tx258/lammps-gran-kokkos
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #2919 from akohlmey/collected-small-changes

Browse Source 
Collected small changes and fixes
This commit is contained in:
Axel Kohlmeyer
2021-09-03 19:03:29 -04:00
committed by GitHub
parent 9085662e2b 6290054e52
commit 801cd647c3
9 changed files with 325 additions and 302 deletions
Download Patch File Download Diff File Expand all files Collapse all files

86
doc/src/Commands_input.rst
View File

@ -1,55 +1,75 @@
LAMMPS input scripts LAMMPS input scripts
==================== ====================
LAMMPS executes by reading commands from a input script (text file), LAMMPS executes calculations by reading commands from a input script (text file), one
one line at a time. When the input script ends, LAMMPS exits. Each line at a time. When the input script ends, LAMMPS exits. This is different
command causes LAMMPS to take some action. It may set an internal from programs that read and process the entire input before starting a calculation.
variable, read in a file, or run a simulation. Most commands have
default settings, which means you only need to use the command if you Each command causes LAMMPS to take some immediate action without regard
wish to change the default. for any commands that may be processed later. Commands may set an
internal variable, read in a file, or run a simulation. These actions
can be grouped into three categories:
a) commands that change a global setting (examples: timestep, newton,
echo, log, thermo, restart),
b) commands that add, modify, remove, or replace "styles" that are
executed during a "run" (examples: pair_style, fix, compute, dump,
thermo_style, pair_modify), and
c) commands that execute a "run" or perform some other computation or
operation (examples: print, run, minimize, temper, write_dump, rerun,
read_data, read_restart)
Commands in category a) have default settings, which means you only
need to use the command if you wish to change the defaults.
In many cases, the ordering of commands in an input script is not In many cases, the ordering of commands in an input script is not
important. However the following rules apply: important, but can have consequences when the global state is changed
between commands in the c) category. The following rules apply:
(1) LAMMPS does not read your entire input script and then perform a (1) LAMMPS does not read your entire input script and then perform a
simulation with all the settings. Rather, the input script is read simulation with all the settings. Rather, the input script is read
one line at a time and each command takes effect when it is read. one line at a time and each command takes effect when it is read.
Thus this sequence of commands: Thus this sequence of commands:
.. code-block:: LAMMPS .. code-block:: LAMMPS
timestep 0.5 timestep 0.5
run 100 run 100
run 100 run 100
does something different than this sequence: does something different than this sequence:
.. code-block:: LAMMPS .. code-block:: LAMMPS
run 100 run 100
timestep 0.5 timestep 0.5
run 100 run 100
In the first case, the specified timestep (0.5 fs) is used for two In the first case, the specified timestep (0.5 fs) is used for two
simulations of 100 timesteps each. In the second case, the default simulations of 100 timesteps each. In the second case, the default
timestep (1.0 fs) is used for the first 100 step simulation and a 0.5 fs timestep (1.0 fs) is used for the first 100 step simulation and a
timestep is used for the second one. 0.5 fs timestep is used for the second one.
(2) Some commands are only valid when they follow other commands. For (2) Some commands are only valid when they follow other commands. For
example you cannot set the temperature of a group of atoms until atoms example you cannot set the temperature of a group of atoms until
have been defined and a group command is used to define which atoms atoms have been defined and a group command is used to define which
belong to the group. atoms belong to the group.
(3) Sometimes command B will use values that can be set by command A. (3) Sometimes command B will use values that can be set by command A.
This means command A must precede command B in the input script if it This means command A must precede command B in the input script if
is to have the desired effect. For example, the it is to have the desired effect. For example, the :doc:`read_data
:doc:`read_data <read_data>` command initializes the system by setting <read_data>` command initializes the system by setting up the
up the simulation box and assigning atoms to processors. If default simulation box and assigning atoms to processors. If default values
values are not desired, the :doc:`processors <processors>` and are not desired, the :doc:`processors <processors>` and
:doc:`boundary <boundary>` commands need to be used before read_data to :doc:`boundary <boundary>` commands need to be used before read_data
tell LAMMPS how to map processors to the simulation box. to tell LAMMPS how to map processors to the simulation box.
Many input script errors are detected by LAMMPS and an ERROR or Many input script errors are detected by LAMMPS and an ERROR or
WARNING message is printed. The :doc:`Errors <Errors>` page gives WARNING message is printed. The :doc:`Errors <Errors>` page gives
more information on what errors mean. The documentation for each more information on what errors mean. The documentation for each
command lists restrictions on how the command can be used. command lists restrictions on how the command can be used.
You can use the :ref:`-skiprun <skiprun>` command line flag
to have LAMMPS skip the execution of any "run", "minimize", or similar
commands to check the entire input for correct syntax to avoid crashes
on typos or syntax errors in long runs.

17
doc/src/Run_options.rst
View File

@ -2,7 +2,7 @@ Command-line options
==================== ====================
At run time, LAMMPS recognizes several optional command-line switches At run time, LAMMPS recognizes several optional command-line switches
which may be used in any order. Either the full word or a one-or-two which may be used in any order. Either the full word or a one or two
letter abbreviation can be used: letter abbreviation can be used:
* :ref:`-e or -echo <echo>` * :ref:`-e or -echo <echo>`
@ -22,6 +22,7 @@ letter abbreviation can be used:
* :ref:`-r2data or -restart2data <restart2data>` * :ref:`-r2data or -restart2data <restart2data>`
* :ref:`-r2dump or -restart2dump <restart2dump>` * :ref:`-r2dump or -restart2dump <restart2dump>`
* :ref:`-sc or -screen <screen>` * :ref:`-sc or -screen <screen>`
* :ref:`-sr or skiprun <skiprun>`
* :ref:`-sf or -suffix <suffix>` * :ref:`-sf or -suffix <suffix>`
* :ref:`-v or -var <var>` * :ref:`-v or -var <var>`
@ -532,6 +533,20 @@ partition screen files file.N.
---------- ----------
.. _skiprun:
**-skiprun**
Insert the command :doc:`timer timerout 0 every 1 <timer>` at the
beginning of an input file or after a :doc:`clear <clear>` command.
This has the effect that the entire LAMMPS input script is processed
without executing actual runs or minimizations (their main loops are
skipped). This can be helpful and convenient to test input scripts for
long running calculations to avoid having them crash after a long time
due to a typo or syntax error.
----------
.. _suffix: .. _suffix:
**-suffix style args** **-suffix style args**

18
doc/src/fix_bond_swap.rst
View File

@ -88,18 +88,18 @@ and bond partners B2 of B1 a is performed. For each pair of A1-A2 and
B1-B2 bonds to be eligible for swapping, the following 4 criteria must B1-B2 bonds to be eligible for swapping, the following 4 criteria must
be met: be met:
(1) All 4 monomers must be in the fix group. 1. All 4 monomers must be in the fix group.
(2) All 4 monomers must be owned by the processor (not ghost atoms). 2. All 4 monomers must be owned by the processor (not ghost atoms).
This insures that another processor does not attempt to swap bonds This insures that another processor does not attempt to swap bonds
involving the same atoms on the same timestep. Note that this also involving the same atoms on the same timestep. Note that this also
means that bond pairs which straddle processor boundaries are not means that bond pairs which straddle processor boundaries are not
eligible for swapping on this step. eligible for swapping on this step.
(3) The distances between 4 pairs of atoms -- (A1,A2), (B1,B2), 3. The distances between 4 pairs of atoms -- (A1,A2), (B1,B2), (A1,B2),
(A1,B2), (B1,A2) -- must all be less than the specified *cutoff*\ . (B1,A2) -- must all be less than the specified *cutoff*.
(4) The molecule IDs of A1 and B1 must be the same (see below). 4. The molecule IDs of A1 and B1 must be the same (see below).
If an eligible B1 partner is found, the energy change due to swapping If an eligible B1 partner is found, the energy change due to swapping
the 2 bonds is computed. This includes changes in pairwise, bond, and the 2 bonds is computed. This includes changes in pairwise, bond, and

2
src/REPLICA/neb.cpp
View File

@ -51,7 +51,7 @@ NEB::NEB(LAMMPS *lmp) : Command(lmp), all(nullptr), rdist(nullptr) {}
NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in, NEB::NEB(LAMMPS *lmp, double etol_in, double ftol_in, int n1steps_in,
int n2steps_in, int nevery_in, double *buf_init, double *buf_final) int n2steps_in, int nevery_in, double *buf_init, double *buf_final)
: Command(lmp), all(nullptr), rdist(nullptr) : Command(lmp), fp(nullptr), all(nullptr), rdist(nullptr)
{ {
double delx,dely,delz; double delx,dely,delz;

2
src/SPIN/neb_spin.cpp
View File

@ -68,7 +68,7 @@ static const char cite_neb_spin[] =
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
NEBSpin::NEBSpin(LAMMPS *lmp) : Command(lmp) { NEBSpin::NEBSpin(LAMMPS *lmp) : Command(lmp), fp(nullptr) {
if (lmp->citeme) lmp->citeme->add(cite_neb_spin); if (lmp->citeme) lmp->citeme->add(cite_neb_spin);
} }

49
src/lammps.cpp
View File

@ -126,6 +126,8 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
cslib = nullptr; cslib = nullptr;
cscomm = 0; cscomm = 0;
skiprunflag = 0;
screen = nullptr; screen = nullptr;
logfile = nullptr; logfile = nullptr;
infile = nullptr; infile = nullptr;
@ -391,6 +393,11 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
screenflag = iarg + 1; screenflag = iarg + 1;
iarg += 2; iarg += 2;
} else if (strcmp(arg[iarg],"-skiprun") == 0 ||
strcmp(arg[iarg],"-sr") == 0) {
skiprunflag = 1;
++iarg;
} else if (strcmp(arg[iarg],"-suffix") == 0 || } else if (strcmp(arg[iarg],"-suffix") == 0 ||
strcmp(arg[iarg],"-sf") == 0) { strcmp(arg[iarg],"-sf") == 0) {
if (iarg+2 > narg) if (iarg+2 > narg)
@ -462,9 +469,8 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
else { else {
universe->uscreen = fopen(arg[screenflag],"w"); universe->uscreen = fopen(arg[screenflag],"w");
if (universe->uscreen == nullptr) if (universe->uscreen == nullptr)
error->universe_one(FLERR,fmt::format("Cannot open universe screen " error->universe_one(FLERR,fmt::format("Cannot open universe screen file {}: {}",
"file {}: {}",arg[screenflag], arg[screenflag],utils::getsyserror()));
utils::getsyserror()));
} }
if (logflag == 0) { if (logflag == 0) {
if (helpflag == 0) { if (helpflag == 0) {
@ -478,9 +484,8 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
else { else {
universe->ulogfile = fopen(arg[logflag],"w"); universe->ulogfile = fopen(arg[logflag],"w");
if (universe->ulogfile == nullptr) if (universe->ulogfile == nullptr)
error->universe_one(FLERR,fmt::format("Cannot open universe log " error->universe_one(FLERR,fmt::format("Cannot open universe log file {}: {}",
"file {}: {}",arg[logflag], arg[logflag],utils::getsyserror()));
utils::getsyserror()));
} }
} }
@ -505,8 +510,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
else if (strcmp(arg[inflag], "none") == 0) infile = stdin; else if (strcmp(arg[inflag], "none") == 0) infile = stdin;
else infile = fopen(arg[inflag],"r"); else infile = fopen(arg[inflag],"r");
if (infile == nullptr) if (infile == nullptr)
error->one(FLERR,"Cannot open input script {}: {}", error->one(FLERR,"Cannot open input script {}: {}",arg[inflag], utils::getsyserror());
arg[inflag], utils::getsyserror());
} }
if ((universe->me == 0) && !helpflag) if ((universe->me == 0) && !helpflag)
@ -530,16 +534,14 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
str = fmt::format("screen.{}",universe->iworld); str = fmt::format("screen.{}",universe->iworld);
screen = fopen(str.c_str(),"w"); screen = fopen(str.c_str(),"w");
if (screen == nullptr) if (screen == nullptr)
error->one(FLERR,"Cannot open screen file {}: {}", error->one(FLERR,"Cannot open screen file {}: {}",str,utils::getsyserror());
str,utils::getsyserror());
} else if (strcmp(arg[screenflag],"none") == 0) { } else if (strcmp(arg[screenflag],"none") == 0) {
screen = nullptr; screen = nullptr;
} else { } else {
str = fmt::format("{}.{}",arg[screenflag],universe->iworld); str = fmt::format("{}.{}",arg[screenflag],universe->iworld);
screen = fopen(str.c_str(),"w"); screen = fopen(str.c_str(),"w");
if (screen == nullptr) if (screen == nullptr)
error->one(FLERR,"Cannot open screen file {}: {}", error->one(FLERR,"Cannot open screen file {}: {}",arg[screenflag],utils::getsyserror());
arg[screenflag],utils::getsyserror());
} }
} else if (strcmp(arg[partscreenflag],"none") == 0) { } else if (strcmp(arg[partscreenflag],"none") == 0) {
screen = nullptr; screen = nullptr;
@ -547,8 +549,7 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
str = fmt::format("{}.{}",arg[partscreenflag],universe->iworld); str = fmt::format("{}.{}",arg[partscreenflag],universe->iworld);
screen = fopen(str.c_str(),"w"); screen = fopen(str.c_str(),"w");
if (screen == nullptr) if (screen == nullptr)
error->one(FLERR,"Cannot open screen file {}: {}", error->one(FLERR,"Cannot open screen file {}: {}",str,utils::getsyserror());
str,utils::getsyserror());
} }
if (partlogflag == 0) { if (partlogflag == 0) {
@ -556,16 +557,14 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
str = fmt::format("log.lammps.{}",universe->iworld); str = fmt::format("log.lammps.{}",universe->iworld);
logfile = fopen(str.c_str(),"w"); logfile = fopen(str.c_str(),"w");
if (logfile == nullptr) if (logfile == nullptr)
error->one(FLERR,"Cannot open logfile {}: {}", error->one(FLERR,"Cannot open logfile {}: {}",str, utils::getsyserror());
str, utils::getsyserror());
} else if (strcmp(arg[logflag],"none") == 0) { } else if (strcmp(arg[logflag],"none") == 0) {
logfile = nullptr; logfile = nullptr;
} else { } else {
str = fmt::format("{}.{}",arg[logflag],universe->iworld); str = fmt::format("{}.{}",arg[logflag],universe->iworld);
logfile = fopen(str.c_str(),"w"); logfile = fopen(str.c_str(),"w");
if (logfile == nullptr) if (logfile == nullptr)
error->one(FLERR,"Cannot open logfile {}: {}", error->one(FLERR,"Cannot open logfile {}: {}",str, utils::getsyserror());
str, utils::getsyserror());
} }
} else if (strcmp(arg[partlogflag],"none") == 0) { } else if (strcmp(arg[partlogflag],"none") == 0) {
logfile = nullptr; logfile = nullptr;
@ -573,15 +572,13 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
str = fmt::format("{}.{}",arg[partlogflag],universe->iworld); str = fmt::format("{}.{}",arg[partlogflag],universe->iworld);
logfile = fopen(str.c_str(),"w"); logfile = fopen(str.c_str(),"w");
if (logfile == nullptr) if (logfile == nullptr)
error->one(FLERR,"Cannot open logfile {}: {}", error->one(FLERR,"Cannot open logfile {}: {}",str, utils::getsyserror());
str, utils::getsyserror());
} }
if (strcmp(arg[inflag], "none") != 0) { if (strcmp(arg[inflag], "none") != 0) {
infile = fopen(arg[inflag],"r"); infile = fopen(arg[inflag],"r");
if (infile == nullptr) if (infile == nullptr)
error->one(FLERR,"Cannot open input script {}: {}", error->one(FLERR,"Cannot open input script {}: {}",arg[inflag], utils::getsyserror());
arg[inflag], utils::getsyserror());
} }
} }
@ -615,12 +612,10 @@ LAMMPS::LAMMPS(int narg, char **arg, MPI_Comm communicator) :
int mpisize; int mpisize;
MPI_Type_size(MPI_LMP_TAGINT,&mpisize); MPI_Type_size(MPI_LMP_TAGINT,&mpisize);
if (mpisize != sizeof(tagint)) if (mpisize != sizeof(tagint))
error->all(FLERR,"MPI_LMP_TAGINT and tagint in " error->all(FLERR,"MPI_LMP_TAGINT and tagint in lmptype.h are not compatible");
"lmptype.h are not compatible");
MPI_Type_size(MPI_LMP_BIGINT,&mpisize); MPI_Type_size(MPI_LMP_BIGINT,&mpisize);
if (mpisize != sizeof(bigint)) if (mpisize != sizeof(bigint))
error->all(FLERR,"MPI_LMP_BIGINT and bigint in " error->all(FLERR,"MPI_LMP_BIGINT and bigint in lmptype.h are not compatible");
"lmptype.h are not compatible");
#ifdef LAMMPS_SMALLBIG #ifdef LAMMPS_SMALLBIG
if (sizeof(smallint) != 4 || sizeof(imageint) != 4 || if (sizeof(smallint) != 4 || sizeof(imageint) != 4 ||
@ -837,6 +832,8 @@ void LAMMPS::create()
void LAMMPS::post_create() void LAMMPS::post_create()
{ {
if (skiprunflag) input->one("timer timeout 0 every 1");
// default package command triggered by "-k on" // default package command triggered by "-k on"
if (kokkos && kokkos->kokkos_exists) input->one("package kokkos"); if (kokkos && kokkos->kokkos_exists) input->one("package kokkos");

68
src/lammps.h
View File

@ -22,52 +22,52 @@ namespace LAMMPS_NS {
class LAMMPS { class LAMMPS {
public: public:
// ptrs to fundamental LAMMPS classes // ptrs to fundamental LAMMPS classes
class Memory *memory; // memory allocation functions class Memory *memory; // memory allocation functions
class Error *error; // error handling class Error *error; // error handling
class Universe *universe; // universe of processors class Universe *universe; // universe of processors
class Input *input; // input script processing class Input *input; // input script processing
// ptrs to top-level LAMMPS-specific classes // ptrs to top-level LAMMPS-specific classes
class Atom *atom; // atom-based quantities class Atom *atom; // atom-based quantities
class Update *update; // integrators/minimizers class Update *update; // integrators/minimizers
class Neighbor *neighbor; // neighbor lists class Neighbor *neighbor; // neighbor lists
class Comm *comm; // inter-processor communication class Comm *comm; // inter-processor communication
class Domain *domain; // simulation box class Domain *domain; // simulation box
class Force *force; // inter-particle forces class Force *force; // inter-particle forces
class Modify *modify; // fixes and computes class Modify *modify; // fixes and computes
class Group *group; // groups of atoms class Group *group; // groups of atoms
class Output *output; // thermo/dump/restart class Output *output; // thermo/dump/restart
class Timer *timer; // CPU timing info class Timer *timer; // CPU timing info
//
class KokkosLMP *kokkos; // KOKKOS accelerator class
class AtomKokkos *atomKK; // KOKKOS version of Atom class
class MemoryKokkos *memoryKK; // KOKKOS version of Memory class
class Python *python; // Python interface
class CiteMe *citeme; // handle citation info
const char *version; // LAMMPS version string = date const char *version; // LAMMPS version string = date
int num_ver; // numeric version id derived from *version* int num_ver; // numeric version id derived from *version*
// that is constructed so that will be greater // that is constructed so that will be greater
// for newer versions in numeric or string // for newer versions in numeric or string
// value comparisons // value comparisons
//
MPI_Comm world; // MPI communicator MPI_Comm world; // MPI communicator
FILE *infile; // infile FILE *infile; // infile
FILE *screen; // screen output FILE *screen; // screen output
FILE *logfile; // logfile FILE *logfile; // logfile
//
double initclock; // wall clock at instantiation double initclock; // wall clock at instantiation
int skiprunflag; // 1 inserts timer command to skip run and minimize loops
char *suffix, *suffix2, *suffixp; // suffixes to add to input script style names char *suffix, *suffix2, *suffixp; // suffixes to add to input script style names
int suffix_enable; // 1 if suffixes are enabled, 0 if disabled int suffix_enable; // 1 if suffixes are enabled, 0 if disabled
char *exename; // pointer to argv[0] char *exename; // pointer to argv[0]
//
char ***packargs; // arguments for cmdline package commands char ***packargs; // arguments for cmdline package commands
int num_package; // number of cmdline package commands int num_package; // number of cmdline package commands
//
int clientserver; // 0 = neither, 1 = client, 2 = server int clientserver; // 0 = neither, 1 = client, 2 = server
void *cslib; // client/server messaging via CSlib void *cslib; // client/server messaging via CSlib
MPI_Comm cscomm; // MPI comm for client+server in mpi/one mode MPI_Comm cscomm; // MPI comm for client+server in mpi/one mode
class KokkosLMP *kokkos; // KOKKOS accelerator class
class AtomKokkos *atomKK; // KOKKOS version of Atom class
class MemoryKokkos *memoryKK; // KOKKOS version of Memory class
class Python *python; // Python interface
class CiteMe *citeme; // handle citation info
const char *match_style(const char *style, const char *name); const char *match_style(const char *style, const char *name);
static const char *installed_packages[]; static const char *installed_packages[];

383
src/pair_lj_cut.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -18,19 +17,20 @@
#include "pair_lj_cut.h" #include "pair_lj_cut.h"
#include <cmath>
#include <cstring>
#include "atom.h" #include "atom.h"
#include "comm.h" #include "comm.h"
#include "error.h"
#include "force.h" #include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "respa.h"
#include "math_const.h" #include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "respa.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst; using namespace MathConst;
@ -68,13 +68,13 @@ PairLJCut::~PairLJCut()
void PairLJCut::compute(int eflag, int vflag) void PairLJCut::compute(int eflag, int vflag)
{ {
int i,j,ii,jj,inum,jnum,itype,jtype; int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double rsq,r2inv,r6inv,forcelj,factor_lj; double rsq, r2inv, r6inv, forcelj, factor_lj;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist, *jlist, *numneigh, **firstneigh;
evdwl = 0.0; evdwl = 0.0;
ev_init(eflag,vflag); ev_init(eflag, vflag);
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
@ -107,32 +107,30 @@ void PairLJCut::compute(int eflag, int vflag)
delx = xtmp - x[j][0]; delx = xtmp - x[j][0];
dely = ytmp - x[j][1]; dely = ytmp - x[j][1];
delz = ztmp - x[j][2]; delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz; rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j]; jtype = type[j];
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
r2inv = 1.0/rsq; r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv; r6inv = r2inv * r2inv * r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv; fpair = factor_lj * forcelj * r2inv;
f[i][0] += delx*fpair; f[i][0] += delx * fpair;
f[i][1] += dely*fpair; f[i][1] += dely * fpair;
f[i][2] += delz*fpair; f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) { if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair; f[j][0] -= delx * fpair;
f[j][1] -= dely*fpair; f[j][1] -= dely * fpair;
f[j][2] -= delz*fpair; f[j][2] -= delz * fpair;
} }
if (eflag) { if (eflag) {
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
offset[itype][jtype];
evdwl *= factor_lj; evdwl *= factor_lj;
} }
if (evflag) ev_tally(i,j,nlocal,newton_pair, if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
evdwl,0.0,fpair,delx,dely,delz);
} }
} }
} }
@ -144,10 +142,10 @@ void PairLJCut::compute(int eflag, int vflag)
void PairLJCut::compute_inner() void PairLJCut::compute_inner()
{ {
int i,j,ii,jj,inum,jnum,itype,jtype; int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,fpair; double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
double rsq,r2inv,r6inv,forcelj,factor_lj,rsw; double rsq, r2inv, r6inv, forcelj, factor_lj, rsw;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist, *jlist, *numneigh, **firstneigh;
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
@ -165,8 +163,8 @@ void PairLJCut::compute_inner()
double cut_out_off = cut_respa[1]; double cut_out_off = cut_respa[1];
double cut_out_diff = cut_out_off - cut_out_on; double cut_out_diff = cut_out_off - cut_out_on;
double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_on_sq = cut_out_on * cut_out_on;
double cut_out_off_sq = cut_out_off*cut_out_off; double cut_out_off_sq = cut_out_off * cut_out_off;
// loop over neighbors of my atoms // loop over neighbors of my atoms
@ -187,26 +185,26 @@ void PairLJCut::compute_inner()
delx = xtmp - x[j][0]; delx = xtmp - x[j][0];
dely = ytmp - x[j][1]; dely = ytmp - x[j][1];
delz = ztmp - x[j][2]; delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz; rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cut_out_off_sq) { if (rsq < cut_out_off_sq) {
r2inv = 1.0/rsq; r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv; r6inv = r2inv * r2inv * r2inv;
jtype = type[j]; jtype = type[j];
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv; fpair = factor_lj * forcelj * r2inv;
if (rsq > cut_out_on_sq) { if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff;
fpair *= 1.0 - rsw*rsw*(3.0 - 2.0*rsw); fpair *= 1.0 - rsw * rsw * (3.0 - 2.0 * rsw);
} }
f[i][0] += delx*fpair; f[i][0] += delx * fpair;
f[i][1] += dely*fpair; f[i][1] += dely * fpair;
f[i][2] += delz*fpair; f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) { if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair; f[j][0] -= delx * fpair;
f[j][1] -= dely*fpair; f[j][1] -= dely * fpair;
f[j][2] -= delz*fpair; f[j][2] -= delz * fpair;
} }
} }
} }
@ -217,10 +215,10 @@ void PairLJCut::compute_inner()
void PairLJCut::compute_middle() void PairLJCut::compute_middle()
{ {
int i,j,ii,jj,inum,jnum,itype,jtype; int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,fpair; double xtmp, ytmp, ztmp, delx, dely, delz, fpair;
double rsq,r2inv,r6inv,forcelj,factor_lj,rsw; double rsq, r2inv, r6inv, forcelj, factor_lj, rsw;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist, *jlist, *numneigh, **firstneigh;
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
@ -241,10 +239,10 @@ void PairLJCut::compute_middle()
double cut_in_diff = cut_in_on - cut_in_off; double cut_in_diff = cut_in_on - cut_in_off;
double cut_out_diff = cut_out_off - cut_out_on; double cut_out_diff = cut_out_off - cut_out_on;
double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_off_sq = cut_in_off * cut_in_off;
double cut_in_on_sq = cut_in_on*cut_in_on; double cut_in_on_sq = cut_in_on * cut_in_on;
double cut_out_on_sq = cut_out_on*cut_out_on; double cut_out_on_sq = cut_out_on * cut_out_on;
double cut_out_off_sq = cut_out_off*cut_out_off; double cut_out_off_sq = cut_out_off * cut_out_off;
// loop over neighbors of my atoms // loop over neighbors of my atoms
@ -265,30 +263,30 @@ void PairLJCut::compute_middle()
delx = xtmp - x[j][0]; delx = xtmp - x[j][0];
dely = ytmp - x[j][1]; dely = ytmp - x[j][1];
delz = ztmp - x[j][2]; delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz; rsq = delx * delx + dely * dely + delz * delz;
if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) { if (rsq < cut_out_off_sq && rsq > cut_in_off_sq) {
r2inv = 1.0/rsq; r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv; r6inv = r2inv * r2inv * r2inv;
jtype = type[j]; jtype = type[j];
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv; fpair = factor_lj * forcelj * r2inv;
if (rsq < cut_in_on_sq) { if (rsq < cut_in_on_sq) {
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff;
fpair *= rsw*rsw*(3.0 - 2.0*rsw); fpair *= rsw * rsw * (3.0 - 2.0 * rsw);
} }
if (rsq > cut_out_on_sq) { if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff; rsw = (sqrt(rsq) - cut_out_on) / cut_out_diff;
fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0); fpair *= 1.0 + rsw * rsw * (2.0 * rsw - 3.0);
} }
f[i][0] += delx*fpair; f[i][0] += delx * fpair;
f[i][1] += dely*fpair; f[i][1] += dely * fpair;
f[i][2] += delz*fpair; f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) { if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair; f[j][0] -= delx * fpair;
f[j][1] -= dely*fpair; f[j][1] -= dely * fpair;
f[j][2] -= delz*fpair; f[j][2] -= delz * fpair;
} }
} }
} }
@ -299,13 +297,13 @@ void PairLJCut::compute_middle()
void PairLJCut::compute_outer(int eflag, int vflag) void PairLJCut::compute_outer(int eflag, int vflag)
{ {
int i,j,ii,jj,inum,jnum,itype,jtype; int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double rsq,r2inv,r6inv,forcelj,factor_lj,rsw; double rsq, r2inv, r6inv, forcelj, factor_lj, rsw;
int *ilist,*jlist,*numneigh,**firstneigh; int *ilist, *jlist, *numneigh, **firstneigh;
evdwl = 0.0; evdwl = 0.0;
ev_init(eflag,vflag); ev_init(eflag, vflag);
double **x = atom->x; double **x = atom->x;
double **f = atom->f; double **f = atom->f;
@ -323,8 +321,8 @@ void PairLJCut::compute_outer(int eflag, int vflag)
double cut_in_on = cut_respa[3]; double cut_in_on = cut_respa[3];
double cut_in_diff = cut_in_on - cut_in_off; double cut_in_diff = cut_in_on - cut_in_off;
double cut_in_off_sq = cut_in_off*cut_in_off; double cut_in_off_sq = cut_in_off * cut_in_off;
double cut_in_on_sq = cut_in_on*cut_in_on; double cut_in_on_sq = cut_in_on * cut_in_on;
// loop over neighbors of my atoms // loop over neighbors of my atoms
@ -345,50 +343,48 @@ void PairLJCut::compute_outer(int eflag, int vflag)
delx = xtmp - x[j][0]; delx = xtmp - x[j][0];
dely = ytmp - x[j][1]; dely = ytmp - x[j][1];
delz = ztmp - x[j][2]; delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz; rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j]; jtype = type[j];
if (rsq < cutsq[itype][jtype]) { if (rsq < cutsq[itype][jtype]) {
if (rsq > cut_in_off_sq) { if (rsq > cut_in_off_sq) {
r2inv = 1.0/rsq; r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv; r6inv = r2inv * r2inv * r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv; fpair = factor_lj * forcelj * r2inv;
if (rsq < cut_in_on_sq) { if (rsq < cut_in_on_sq) {
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff; rsw = (sqrt(rsq) - cut_in_off) / cut_in_diff;
fpair *= rsw*rsw*(3.0 - 2.0*rsw); fpair *= rsw * rsw * (3.0 - 2.0 * rsw);
} }
f[i][0] += delx*fpair; f[i][0] += delx * fpair;
f[i][1] += dely*fpair; f[i][1] += dely * fpair;
f[i][2] += delz*fpair; f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) { if (newton_pair || j < nlocal) {
f[j][0] -= delx*fpair; f[j][0] -= delx * fpair;
f[j][1] -= dely*fpair; f[j][1] -= dely * fpair;
f[j][2] -= delz*fpair; f[j][2] -= delz * fpair;
} }
} }
if (eflag) { if (eflag) {
r2inv = 1.0/rsq; r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv; r6inv = r2inv * r2inv * r2inv;
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - evdwl = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
offset[itype][jtype];
evdwl *= factor_lj; evdwl *= factor_lj;
} }
if (vflag) { if (vflag) {
if (rsq <= cut_in_off_sq) { if (rsq <= cut_in_off_sq) {
r2inv = 1.0/rsq; r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv; r6inv = r2inv * r2inv * r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
fpair = factor_lj*forcelj*r2inv; fpair = factor_lj * forcelj * r2inv;
} else if (rsq < cut_in_on_sq) } else if (rsq < cut_in_on_sq)
fpair = factor_lj*forcelj*r2inv; fpair = factor_lj * forcelj * r2inv;
} }
if (evflag) ev_tally(i,j,nlocal,newton_pair, if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
evdwl,0.0,fpair,delx,dely,delz);
} }
} }
} }
@ -401,23 +397,22 @@ void PairLJCut::compute_outer(int eflag, int vflag)
void PairLJCut::allocate() void PairLJCut::allocate()
{ {
allocated = 1; allocated = 1;
int n = atom->ntypes; int n = atom->ntypes + 1;
memory->create(setflag,n+1,n+1,"pair:setflag"); memory->create(setflag, n, n, "pair:setflag");
for (int i = 1; i <= n; i++) for (int i = 1; i < n; i++)
for (int j = i; j <= n; j++) for (int j = i; j < n; j++) setflag[i][j] = 0;
setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq"); memory->create(cutsq, n, n, "pair:cutsq");
memory->create(cut,n+1,n+1,"pair:cut"); memory->create(cut, n, n, "pair:cut");
memory->create(epsilon,n+1,n+1,"pair:epsilon"); memory->create(epsilon, n, n, "pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma"); memory->create(sigma, n, n, "pair:sigma");
memory->create(lj1,n+1,n+1,"pair:lj1"); memory->create(lj1, n, n, "pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2"); memory->create(lj2, n, n, "pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3"); memory->create(lj3, n, n, "pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4"); memory->create(lj4, n, n, "pair:lj4");
memory->create(offset,n+1,n+1,"pair:offset"); memory->create(offset, n, n, "pair:offset");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -426,14 +421,14 @@ void PairLJCut::allocate()
void PairLJCut::settings(int narg, char **arg) void PairLJCut::settings(int narg, char **arg)
{ {
if (narg != 1) error->all(FLERR,"Illegal pair_style command"); if (narg != 1) error->all(FLERR, "Illegal pair_style command");
cut_global = utils::numeric(FLERR,arg[0],false,lmp); cut_global = utils::numeric(FLERR, arg[0], false, lmp);
// reset cutoffs that have been explicitly set // reset cutoffs that have been explicitly set
if (allocated) { if (allocated) {
int i,j; int i, j;
for (i = 1; i <= atom->ntypes; i++) for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) for (j = i; j <= atom->ntypes; j++)
if (setflag[i][j]) cut[i][j] = cut_global; if (setflag[i][j]) cut[i][j] = cut_global;
@ -446,23 +441,22 @@ void PairLJCut::settings(int narg, char **arg)
void PairLJCut::coeff(int narg, char **arg) void PairLJCut::coeff(int narg, char **arg)
{ {
if (narg < 4 || narg > 5) if (narg < 4 || narg > 5) error->all(FLERR, "Incorrect args for pair coefficients");
error->all(FLERR,"Incorrect args for pair coefficients");
if (!allocated) allocate(); if (!allocated) allocate();
int ilo,ihi,jlo,jhi; int ilo, ihi, jlo, jhi;
utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error); utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error); utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
double epsilon_one = utils::numeric(FLERR,arg[2],false,lmp); double epsilon_one = utils::numeric(FLERR, arg[2], false, lmp);
double sigma_one = utils::numeric(FLERR,arg[3],false,lmp); double sigma_one = utils::numeric(FLERR, arg[3], false, lmp);
double cut_one = cut_global; double cut_one = cut_global;
if (narg == 5) cut_one = utils::numeric(FLERR,arg[4],false,lmp); if (narg == 5) cut_one = utils::numeric(FLERR, arg[4], false, lmp);
int count = 0; int count = 0;
for (int i = ilo; i <= ihi; i++) { for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo,i); j <= jhi; j++) { for (int j = MAX(jlo, i); j <= jhi; j++) {
epsilon[i][j] = epsilon_one; epsilon[i][j] = epsilon_one;
sigma[i][j] = sigma_one; sigma[i][j] = sigma_one;
cut[i][j] = cut_one; cut[i][j] = cut_one;
@ -471,7 +465,7 @@ void PairLJCut::coeff(int narg, char **arg)
} }
} }
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -485,12 +479,12 @@ void PairLJCut::init_style()
int irequest; int irequest;
int respa = 0; int respa = 0;
if (update->whichflag == 1 && utils::strmatch(update->integrate_style,"^respa")) { if (update->whichflag == 1 && utils::strmatch(update->integrate_style, "^respa")) {
if (((Respa *) update->integrate)->level_inner >= 0) respa = 1; if (((Respa *) update->integrate)->level_inner >= 0) respa = 1;
if (((Respa *) update->integrate)->level_middle >= 0) respa = 2; if (((Respa *) update->integrate)->level_middle >= 0) respa = 2;
} }
irequest = neighbor->request(this,instance_me); irequest = neighbor->request(this, instance_me);
if (respa >= 1) { if (respa >= 1) {
neighbor->requests[irequest]->respaouter = 1; neighbor->requests[irequest]->respaouter = 1;
@ -500,10 +494,11 @@ void PairLJCut::init_style()
// set rRESPA cutoffs // set rRESPA cutoffs
if (utils::strmatch(update->integrate_style,"^respa") && if (utils::strmatch(update->integrate_style, "^respa") &&
((Respa *) update->integrate)->level_inner >= 0) ((Respa *) update->integrate)->level_inner >= 0)
cut_respa = ((Respa *) update->integrate)->cutoff; cut_respa = ((Respa *) update->integrate)->cutoff;
else cut_respa = nullptr; else
cut_respa = nullptr;
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -513,21 +508,21 @@ void PairLJCut::init_style()
double PairLJCut::init_one(int i, int j) double PairLJCut::init_one(int i, int j)
{ {
if (setflag[i][j] == 0) { if (setflag[i][j] == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], epsilon[i][j] = mix_energy(epsilon[i][i], epsilon[j][j], sigma[i][i], sigma[j][j]);
sigma[i][i],sigma[j][j]); sigma[i][j] = mix_distance(sigma[i][i], sigma[j][j]);
sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
cut[i][j] = mix_distance(cut[i][i],cut[j][j]);
} }
lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j], 12.0);
lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j], 12.0);
lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j], 6.0);
if (offset_flag && (cut[i][j] > 0.0)) { if (offset_flag && (cut[i][j] > 0.0)) {
double ratio = sigma[i][j] / cut[i][j]; double ratio = sigma[i][j] / cut[i][j];
offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio, 12.0) - pow(ratio, 6.0));
} else offset[i][j] = 0.0; } else
offset[i][j] = 0.0;
lj1[j][i] = lj1[i][j]; lj1[j][i] = lj1[i][j];
lj2[j][i] = lj2[i][j]; lj2[j][i] = lj2[i][j];
@ -538,7 +533,7 @@ double PairLJCut::init_one(int i, int j)
// check interior rRESPA cutoff // check interior rRESPA cutoff
if (cut_respa && cut[i][j] < cut_respa[3]) if (cut_respa && cut[i][j] < cut_respa[3])
error->all(FLERR,"Pair cutoff < Respa interior cutoff"); error->all(FLERR, "Pair cutoff < Respa interior cutoff");
// compute I,J contribution to long-range tail correction // compute I,J contribution to long-range tail correction
// count total # of atoms of type I and J via Allreduce // count total # of atoms of type I and J via Allreduce
@ -547,23 +542,22 @@ double PairLJCut::init_one(int i, int j)
int *type = atom->type; int *type = atom->type;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double count[2],all[2]; double count[2], all[2];
count[0] = count[1] = 0.0; count[0] = count[1] = 0.0;
for (int k = 0; k < nlocal; k++) { for (int k = 0; k < nlocal; k++) {
if (type[k] == i) count[0] += 1.0; if (type[k] == i) count[0] += 1.0;
if (type[k] == j) count[1] += 1.0; if (type[k] == j) count[1] += 1.0;
} }
MPI_Allreduce(count,all,2,MPI_DOUBLE,MPI_SUM,world); MPI_Allreduce(count, all, 2, MPI_DOUBLE, MPI_SUM, world);
double sig2 = sigma[i][j]*sigma[i][j]; double sig2 = sigma[i][j] * sigma[i][j];
double sig6 = sig2*sig2*sig2; double sig6 = sig2 * sig2 * sig2;
double rc3 = cut[i][j]*cut[i][j]*cut[i][j]; double rc3 = cut[i][j] * cut[i][j] * cut[i][j];
double rc6 = rc3*rc3; double rc6 = rc3 * rc3;
double rc9 = rc3*rc6; double rc9 = rc3 * rc6;
etail_ij = 8.0*MY_PI*all[0]*all[1]*epsilon[i][j] * double prefactor = 8.0 * MY_PI * all[0] * all[1] * epsilon[i][j] * sig6 / (9.0 * rc9);
sig6 * (sig6 - 3.0*rc6) / (9.0*rc9); etail_ij = prefactor * (sig6 - 3.0 * rc6);
ptail_ij = 16.0*MY_PI*all[0]*all[1]*epsilon[i][j] * ptail_ij = 2.0 * prefactor * (2.0 * sig6 - 3.0 * rc6);
sig6 * (2.0*sig6 - 3.0*rc6) / (9.0*rc9);
} }
return cut[i][j]; return cut[i][j];
@ -577,14 +571,14 @@ void PairLJCut::write_restart(FILE *fp)
{ {
write_restart_settings(fp); write_restart_settings(fp);
int i,j; int i, j;
for (i = 1; i <= atom->ntypes; i++) for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) { for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j],sizeof(int),1,fp); fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) { if (setflag[i][j]) {
fwrite(&epsilon[i][j],sizeof(double),1,fp); fwrite(&epsilon[i][j], sizeof(double), 1, fp);
fwrite(&sigma[i][j],sizeof(double),1,fp); fwrite(&sigma[i][j], sizeof(double), 1, fp);
fwrite(&cut[i][j],sizeof(double),1,fp); fwrite(&cut[i][j], sizeof(double), 1, fp);
} }
} }
} }
@ -598,21 +592,21 @@ void PairLJCut::read_restart(FILE *fp)
read_restart_settings(fp); read_restart_settings(fp);
allocate(); allocate();
int i,j; int i, j;
int me = comm->me; int me = comm->me;
for (i = 1; i <= atom->ntypes; i++) for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) { for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR,&setflag[i][j],sizeof(int),1,fp,nullptr,error); if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) { if (setflag[i][j]) {
if (me == 0) { if (me == 0) {
utils::sfread(FLERR,&epsilon[i][j],sizeof(double),1,fp,nullptr,error); utils::sfread(FLERR, &epsilon[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR,&sigma[i][j],sizeof(double),1,fp,nullptr,error); utils::sfread(FLERR, &sigma[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR,&cut[i][j],sizeof(double),1,fp,nullptr,error); utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
} }
MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&epsilon[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&sigma[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
} }
} }
} }
@ -623,10 +617,10 @@ void PairLJCut::read_restart(FILE *fp)
void PairLJCut::write_restart_settings(FILE *fp) void PairLJCut::write_restart_settings(FILE *fp)
{ {
fwrite(&cut_global,sizeof(double),1,fp); fwrite(&cut_global, sizeof(double), 1, fp);
fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&offset_flag, sizeof(int), 1, fp);
fwrite(&mix_flag,sizeof(int),1,fp); fwrite(&mix_flag, sizeof(int), 1, fp);
fwrite(&tail_flag,sizeof(int),1,fp); fwrite(&tail_flag, sizeof(int), 1, fp);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -637,15 +631,15 @@ void PairLJCut::read_restart_settings(FILE *fp)
{ {
int me = comm->me; int me = comm->me;
if (me == 0) { if (me == 0) {
utils::sfread(FLERR,&cut_global,sizeof(double),1,fp,nullptr,error); utils::sfread(FLERR, &cut_global, sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR,&offset_flag,sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR, &offset_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR,&mix_flag,sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR, &mix_flag, sizeof(int), 1, fp, nullptr, error);
utils::sfread(FLERR,&tail_flag,sizeof(int),1,fp,nullptr,error); utils::sfread(FLERR, &tail_flag, sizeof(int), 1, fp, nullptr, error);
} }
MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_global, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&offset_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&mix_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag, 1, MPI_INT, 0, world);
MPI_Bcast(&tail_flag,1,MPI_INT,0,world); MPI_Bcast(&tail_flag, 1, MPI_INT, 0, world);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -654,8 +648,7 @@ void PairLJCut::read_restart_settings(FILE *fp)
void PairLJCut::write_data(FILE *fp) void PairLJCut::write_data(FILE *fp)
{ {
for (int i = 1; i <= atom->ntypes; i++) for (int i = 1; i <= atom->ntypes; i++) fprintf(fp, "%d %g %g\n", i, epsilon[i][i], sigma[i][i]);
fprintf(fp,"%d %g %g\n",i,epsilon[i][i],sigma[i][i]);
} }
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
@ -666,25 +659,23 @@ void PairLJCut::write_data_all(FILE *fp)
{ {
for (int i = 1; i <= atom->ntypes; i++) for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++) for (int j = i; j <= atom->ntypes; j++)
fprintf(fp,"%d %d %g %g %g\n",i,j,epsilon[i][j],sigma[i][j],cut[i][j]); fprintf(fp, "%d %d %g %g %g\n", i, j, epsilon[i][j], sigma[i][j], cut[i][j]);
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
double PairLJCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq, double PairLJCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
double /*factor_coul*/, double factor_lj, double /*factor_coul*/, double factor_lj, double &fforce)
double &fforce)
{ {
double r2inv,r6inv,forcelj,philj; double r2inv, r6inv, forcelj, philj;
r2inv = 1.0/rsq; r2inv = 1.0 / rsq;
r6inv = r2inv*r2inv*r2inv; r6inv = r2inv * r2inv * r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]); forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
fforce = factor_lj*forcelj*r2inv; fforce = factor_lj * forcelj * r2inv;
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) - philj = r6inv * (lj3[itype][jtype] * r6inv - lj4[itype][jtype]) - offset[itype][jtype];
offset[itype][jtype]; return factor_lj * philj;
return factor_lj*philj;
} }
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
@ -692,7 +683,7 @@ double PairLJCut::single(int /*i*/, int /*j*/, int itype, int jtype, double rsq,
void *PairLJCut::extract(const char *str, int &dim) void *PairLJCut::extract(const char *str, int &dim)
{ {
dim = 2; dim = 2;
if (strcmp(str,"epsilon") == 0) return (void *) epsilon; if (strcmp(str, "epsilon") == 0) return (void *) epsilon;
if (strcmp(str,"sigma") == 0) return (void *) sigma; if (strcmp(str, "sigma") == 0) return (void *) sigma;
return nullptr; return nullptr;
} }

2
src/variable.cpp
View File

@ -4720,7 +4720,7 @@ int Variable::parse_args(char *str, char **args)
while (ptr && narg < MAXFUNCARG) { while (ptr && narg < MAXFUNCARG) {
ptrnext = find_next_comma(ptr); ptrnext = find_next_comma(ptr);
if (ptrnext) *ptrnext = '\0'; if (ptrnext) *ptrnext = '\0';
args[narg] = utils::strdup(ptr); args[narg] = utils::strdup(utils::trim(ptr));
narg++; narg++;
ptr = ptrnext; ptr = ptrnext;
if (ptr) ptr++; if (ptr) ptr++;