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

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

Collected small changes and fixes
This commit is contained in:
Axel Kohlmeyer 2021-03-29 14:05:08 -04:00 committed by GitHub
parent dbc30f9c19 1e3f1c584c
commit 180d2895c4
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
153 changed files with 1275 additions and 3487 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
cmake/CMakeLists.txt
View File

@ -31,7 +31,7 @@ endif()
# If enabled, no need to use LD_LIBRARY_PATH / DYLD_LIBRARY_PATH when installed
option(LAMMPS_INSTALL_RPATH "Set runtime path for shared libraries linked to LAMMPS binaries" OFF)
if (LAMMPS_INSTALL_RPATH)
if(LAMMPS_INSTALL_RPATH)
set(CMAKE_INSTALL_RPATH ${CMAKE_INSTALL_FULL_LIBDIR})
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH ON)
endif()

2
doc/src/Commands_bond.rst
View File

@ -126,7 +126,7 @@ OPT.
* :doc:`quadratic (o) <dihedral_quadratic>`
* :doc:`spherical <dihedral_spherical>`
* :doc:`table (o) <dihedral_table>`
* :doc:`table/cut <dihedral_table_cut>`
* :doc:`table/cut <dihedral_table>`
.. _improper:

1
doc/src/Howto_bioFF.rst
View File

@ -91,6 +91,7 @@ documentation for the formula it computes.
* :doc:`bond_style <bond_harmonic>` harmonic
* :doc:`bond_style <bond_morse>` morse
* :doc:`angle_style <angle_cosine_squared>` cosine/squared
* :doc:`angle_style <angle_harmonic>` harmonic
* :doc:`angle_style <angle_cosine>` cosine
* :doc:`angle_style <angle_cosine_periodic>` cosine/periodic

16
doc/src/angle_cosine_periodic.rst
View File

@ -24,20 +24,17 @@ Examples
Description
"""""""""""
The *cosine/periodic* angle style uses the following potential, which
is commonly used in the :doc:`DREIDING <Howto_bioFF>` force field,
particularly for organometallic systems where :math:`n` = 4 might be used
The *cosine/periodic* angle style uses the following potential, which may be
particularly used for organometallic systems where :math:`n` = 4 might be used
for an octahedral complex and :math:`n` = 3 might be used for a trigonal
center:
.. math::
E = C \left[ 1 - B(-1)^n\cos\left( n\theta\right) \right]
E = \frac{2.0}{n^2} * C \left[ 1 - B(-1)^n\cos\left( n\theta\right) \right]
where :math:`C`, :math:`B` and :math:`n` are coefficients defined for each angle type.
See :ref:`(Mayo) <cosine-Mayo>` for a description of the DREIDING force field
The following coefficients must be defined for each angle type via the
:doc:`angle_coeff <angle_coeff>` command as in the example above, or in
the data file or restart files read by the :doc:`read_data <read_data>`
@ -47,10 +44,9 @@ or :doc:`read_restart <read_restart>` commands:
* :math:`B` = 1 or -1
* :math:`n` = 1, 2, 3, 4, 5 or 6 for periodicity
Note that the prefactor :math:`C` is specified and not the overall force
constant :math:`K = \frac{C}{n^2}`. When :math:`B = 1`, it leads to a minimum for the
linear geometry. When :math:`B = -1`, it leads to a maximum for the linear
geometry.
Note that the prefactor :math:`C` is specified as coefficient and not the overall force
constant :math:`K = \frac{2 C}{n^2}`. When :math:`B = 1`, it leads to a minimum for the
linear geometry. When :math:`B = -1`, it leads to a maximum for the linear geometry.
----------

14
doc/src/angle_cosine_squared.rst
View File

@ -30,8 +30,11 @@ The *cosine/squared* angle style uses the potential
E = K [\cos(\theta) - \cos(\theta_0)]^2
where :math:`\theta_0` is the equilibrium value of the angle, and :math:`K` is a
prefactor. Note that the usual 1/2 factor is included in :math:`K`.
, which is commonly used in the :doc:`DREIDING <Howto_bioFF>` force field,
where :math:`\theta_0` is the equilibrium value of the angle, and :math:`K`
is a prefactor. Note that the usual 1/2 factor is included in :math:`K`.
See :ref:`(Mayo) <cosine-Mayo>` for a description of the DREIDING force field.
The following coefficients must be defined for each angle type via the
:doc:`angle_coeff <angle_coeff>` command as in the example above, or in
@ -66,3 +69,10 @@ Default
"""""""
none
----------
.. _cosine-Mayo:
**(Mayo)** Mayo, Olfason, Goddard III, J Phys Chem, 94, 8897-8909
(1990).

2
doc/src/dihedral_style.rst
View File

@ -113,7 +113,7 @@ more of (g,i,k,o,t) to indicate which accelerated styles exist.
* :doc:`quadratic <dihedral_quadratic>` - dihedral with quadratic term in angle
* :doc:`spherical <dihedral_spherical>` - dihedral which includes angle terms to avoid singularities
* :doc:`table <dihedral_table>` - tabulated dihedral
* :doc:`table/cut <dihedral_table_cut>` - tabulated dihedral with analytic cutoff
* :doc:`table/cut <dihedral_table>` - tabulated dihedral with analytic cutoff
----------

89
doc/src/dihedral_table.rst
View File

@ -1,19 +1,24 @@
.. index:: dihedral_style table
.. index:: dihedral_style table/omp
.. index:: dihedral_style table/cut
dihedral_style table command
============================
Accelerator Variants: *table/omp*
dihedral_style table/cut command
================================
Syntax
""""""
.. code-block:: LAMMPS
dihedral_style table style Ntable
dihedral_style style interp Ntable
* style = *linear* or *spline* = method of interpolation
* style = *table* or *table/cut*
* interp = *linear* or *spline* = method of interpolation
* Ntable = size of the internal lookup table
Examples
@ -26,13 +31,21 @@ Examples
dihedral_coeff 1 file.table DIH_TABLE1
dihedral_coeff 2 file.table DIH_TABLE2
dihedral_style table/cut spline 400
dihedral_style table/cut linear 1000
dihedral_coeff 1 aat 1.0 177 180 file.table DIH_TABLE1
dihedral_coeff 2 aat 0.5 170 180 file.table DIH_TABLE2
Description
"""""""""""
The *table* dihedral style creates interpolation tables of length
*Ntable* from dihedral potential and derivative values listed in a
file(s) as a function of the dihedral angle "phi". The files are read
by the :doc:`dihedral_coeff <dihedral_coeff>` command.
The *table* and *table/cut* dihedral styles create interpolation tables
of length *Ntable* from dihedral potential and derivative values listed
in a file(s) as a function of the dihedral angle "phi". The files are
read by the :doc:`dihedral_coeff <dihedral_coeff>` command. For
dihedral style *table/cut* additionally an analytic cutoff that is
quadratic in the bond-angle (theta) is applied in order to regularize
the dihedral interaction.
The interpolation tables are created by fitting cubic splines to the
file values and interpolating energy and derivative values at each of
@ -51,16 +64,53 @@ interpolated table. For a given dihedral angle (phi), the appropriate
coefficients are chosen from this list, and a cubic polynomial is used
to compute the energy and the derivative at this angle.
The following coefficients must be defined for each dihedral type via
the :doc:`dihedral_coeff <dihedral_coeff>` command as in the example
above.
For dihedral style *table* the following coefficients must be defined
for each dihedral type via the :doc:`dihedral_coeff <dihedral_coeff>`
command as in the example above.
* filename
* keyword
The filename specifies a file containing tabulated energy and
derivative values. The keyword specifies a section of the file. The
format of this file is described below.
The filename specifies a file containing tabulated energy and derivative
values. The keyword specifies which section of the file to read. The
format of this file is the same for both dihedral styles and described
below.
For dihedral style *table/cut* the following coefficients must be
defined for each dihedral type via the :doc:`dihedral_coeff
<dihedral_coeff>` command as in the example above.
* style (aat)
* cutoff prefactor
* cutoff angle1
* cutoff angle2
* filename
* keyword
The cutoff dihedral style uses a tabulated dihedral interaction with a
cutoff function:
.. math::
f(\theta) & = K \qquad\qquad\qquad\qquad\qquad\qquad \theta < \theta_1 \\
f(\theta) & = K \left(1-\frac{(\theta - \theta_1)^2}{(\theta_2 - \theta_1)^2}\right) \qquad \theta_1 < \theta < \theta_2
The cutoff specifies an prefactor to the cutoff function. While this
value would ordinarily equal 1 there may be situations where the value
should change.
The cutoff :math:`\theta_1` specifies the angle (in degrees) below which
the dihedral interaction is unmodified, i.e. the cutoff function is 1.
The cutoff function is applied between :math:`\theta_1` and
:math:`\theta_2`, which is the angle at which the cutoff function drops
to zero. The value of zero effectively "turns off" the dihedral
interaction.
The filename specifies a file containing tabulated energy and derivative
values. The keyword specifies which section of the file to read. The
format of this file is the same for both dihedral styles and described
below.
----------
@ -182,18 +232,19 @@ that matches the specified keyword.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This dihedral style writes the settings for the "dihedral_style table"
command to :doc:`binary restart files <restart>`, so a dihedral_style
command does not need to specified in an input script that reads a
restart file. However, the coefficient information is not stored in
the restart file, since it is tabulated in the potential files. Thus,
These dihedral styles write the settings for the "dihedral_style table"
or "dihedral_style table/cut" command to :doc:`binary restart files
<restart>`, so a dihedral_style command does not need to specified in an
input script that reads a restart file. However, the coefficient
information loaded from the table file(s) is not stored in the restart
file, since it is tabulated in the potential files. Thus, suitable
dihedral_coeff commands do need to be specified in the restart input
script.
script after reading the restart file.
Restrictions
""""""""""""
This dihedral style can only be used if LAMMPS was built with the
These dihedral styles can only be used if LAMMPS was built with the
USER-MISC package. See the :doc:`Build package <Build_package>` doc
page for more info.

229
doc/src/dihedral_table_cut.rst
View File

@ -1,229 +0,0 @@
.. index:: dihedral_style table/cut
dihedral_style table/cut command
================================
Syntax
""""""
.. code-block:: LAMMPS
dihedral_style table/cut style Ntable
* style = *linear* or *spline* = method of interpolation
* Ntable = size of the internal lookup table
Examples
""""""""
.. code-block:: LAMMPS
dihedral_style table/cut spline 400
dihedral_style table/cut linear 1000
dihedral_coeff 1 aat 1.0 177 180 file.table DIH_TABLE1
dihedral_coeff 2 aat 0.5 170 180 file.table DIH_TABLE2
Description
"""""""""""
The *table/cut* dihedral style creates interpolation tables of length
*Ntable* from dihedral potential and derivative values listed in a
file(s) as a function of the dihedral angle "phi". In addition, an
analytic cutoff that is quadratic in the bond-angle (theta) is applied
in order to regularize the dihedral interaction. The dihedral table
files are read by the :doc:`dihedral_coeff <dihedral_coeff>` command.
The interpolation tables are created by fitting cubic splines to the
file values and interpolating energy and derivative values at each of
*Ntable* dihedral angles. During a simulation, these tables are used
to interpolate energy and force values on individual atoms as
needed. The interpolation is done in one of 2 styles: *linear* or
*spline*\ .
For the *linear* style, the dihedral angle (phi) is used to find 2
surrounding table values from which an energy or its derivative is
computed by linear interpolation.
For the *spline* style, cubic spline coefficients are computed and
stored at each of the *Ntable* evenly-spaced values in the
interpolated table. For a given dihedral angle (phi), the appropriate
coefficients are chosen from this list, and a cubic polynomial is used
to compute the energy and the derivative at this angle.
The following coefficients must be defined for each dihedral type via
the :doc:`dihedral_coeff <dihedral_coeff>` command as in the example
above.
* style (aat)
* cutoff prefactor
* cutoff angle1
* cutoff angle2
* filename
* keyword
The cutoff dihedral style uses a tabulated dihedral interaction with a
cutoff function:
.. math::
f(\theta) & = K \qquad\qquad\qquad\qquad\qquad\qquad \theta < \theta_1 \\
f(\theta) & = K \left(1-\frac{(\theta - \theta_1)^2}{(\theta_2 - \theta_1)^2}\right) \qquad \theta_1 < \theta < \theta_2
The cutoff specifies an prefactor to the cutoff function. While this value
would ordinarily equal 1 there may be situations where the value should change.
The cutoff :math:`\theta_1` specifies the angle (in degrees) below which the dihedral
interaction is unmodified, i.e. the cutoff function is 1.
The cutoff function is applied between :math:`\theta_1` and :math:`\theta_2`, which is
the angle at which the cutoff function drops to zero. The value of zero effectively
"turns off" the dihedral interaction.
The filename specifies a file containing tabulated energy and
derivative values. The keyword specifies a section of the file. The
format of this file is described below.
----------
The format of a tabulated file is as follows (without the
parenthesized comments). It can begin with one or more comment
or blank lines.
.. parsed-literal::
# Table of the potential and its negative derivative
DIH_TABLE1 (keyword is the first text on line)
N 30 DEGREES (N, NOF, DEGREES, RADIANS, CHECKU/F)
(blank line)
1 -168.0 -1.40351172223 0.0423346818422
2 -156.0 -1.70447981034 0.00811786522531
3 -144.0 -1.62956100432 -0.0184129719987
...
30 180.0 -0.707106781187 0.0719306095245
# Example 2: table of the potential. Forces omitted
DIH_TABLE2
N 30 NOF CHECKU testU.dat CHECKF testF.dat
1 -168.0 -1.40351172223
2 -156.0 -1.70447981034
3 -144.0 -1.62956100432
...
30 180.0 -0.707106781187
A section begins with a non-blank line whose first character is not a
"#"; blank lines or lines starting with "#" can be used as comments
between sections. The first line begins with a keyword which
identifies the section. The line can contain additional text, but the
initial text must match the argument specified in the
:doc:`dihedral_coeff <dihedral_coeff>` command. The next line lists (in
any order) one or more parameters for the table. Each parameter is a
keyword followed by one or more numeric values.
Following a blank line, the next N lines list the tabulated values. On
each line, the first value is the index from 1 to N, the second value is
the angle value, the third value is the energy (in energy units), and
the fourth is -dE/d(phi) also in energy units). The third term is the
energy of the 4-atom configuration for the specified angle. The fourth
term (when present) is the negative derivative of the energy with
respect to the angle (in degrees, or radians depending on whether the
user selected DEGREES or RADIANS). Thus the units of the last term
are still energy, not force. The dihedral angle values must increase
from one line to the next.
Dihedral table splines are cyclic. There is no discontinuity at 180
degrees (or at any other angle). Although in the examples above, the
angles range from -180 to 180 degrees, in general, the first angle in
the list can have any value (positive, zero, or negative). However
the *range* of angles represented in the table must be *strictly* less
than 360 degrees (2pi radians) to avoid angle overlap. (You may not
supply entries in the table for both 180 and -180, for example.) If
the user's table covers only a narrow range of dihedral angles,
strange numerical behavior can occur in the large remaining gap.
**Parameters:**
The parameter "N" is required and its value is the number of table
entries that follow. Note that this may be different than the N
specified in the :doc:`dihedral_style table <dihedral_style>` command.
Let *Ntable* is the number of table entries requested dihedral_style
command, and let *Nfile* be the parameter following "N" in the
tabulated file ("30" in the sparse example above). What LAMMPS does
is a preliminary interpolation by creating splines using the *Nfile*
tabulated values as nodal points. It uses these to interpolate as
needed to generate energy and derivative values at *Ntable* different
points (which are evenly spaced over a 360 degree range, even if the
angles in the file are not). The resulting tables of length *Ntable*
are then used as described above, when computing energy and force for
individual dihedral angles and their atoms. This means that if you
want the interpolation tables of length *Ntable* to match exactly what
is in the tabulated file (with effectively nopreliminary
interpolation), you should set *Ntable* = *Nfile*\ . To insure the
nodal points in the user's file are aligned with the interpolated
table entries, the angles in the table should be integer multiples of
360/\ *Ntable* degrees, or 2\*PI/\ *Ntable* radians (depending on your
choice of angle units).
The optional "NOF" keyword allows the user to omit the forces
(negative energy derivatives) from the table file (normally located in
the fourth column). In their place, forces will be calculated
automatically by differentiating the potential energy function
indicated by the third column of the table (using either linear or
spline interpolation).
The optional "DEGREES" keyword allows the user to specify angles in
degrees instead of radians (default).
The optional "RADIANS" keyword allows the user to specify angles in
radians instead of degrees. (Note: This changes the way the forces
are scaled in the fourth column of the data file.)
The optional "CHECKU" keyword is followed by a filename. This allows
the user to save all of the *Ntable* different entries in the
interpolated energy table to a file to make sure that the interpolated
function agrees with the user's expectations. (Note: You can
temporarily increase the *Ntable* parameter to a high value for this
purpose. "\ *Ntable*\ " is explained above.)
The optional "CHECKF" keyword is analogous to the "CHECKU" keyword.
It is followed by a filename, and it allows the user to check the
interpolated force table. This option is available even if the user
selected the "NOF" option.
Note that one file can contain many sections, each with a tabulated
potential. LAMMPS reads the file section by section until it finds one
that matches the specified keyword.
Restart, fix_modify, output, run start/stop, minimize info
"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
This dihedral style writes the settings for the "dihedral_style table/cut"
command to :doc:`binary restart files <restart>`, so a dihedral_style
command does not need to specified in an input script that reads a
restart file. However, the coefficient information is not stored in
the restart file, since it is tabulated in the potential files. Thus,
dihedral_coeff commands do need to be specified in the restart input
script.
Restrictions
""""""""""""
This dihedral style can only be used if LAMMPS was built with the
USER-MISC package. See the :doc:`Build package <Build_package>` doc
page for more info.
Related commands
""""""""""""""""
:doc:`dihedral_coeff <dihedral_coeff>`, :doc:`dihedral_style table <dihedral_table>`
Default
"""""""
none
.. _dihedralcut-Salerno:
**(Salerno)** Salerno, Bernstein, J Chem Theory Comput, --, ---- (2018).

2
doc/src/pair_polymorphic.rst
View File

@ -319,7 +319,7 @@ This pair style is part of the MANYBODY package. It is only enabled if
LAMMPS was built with that package. See the :doc:`Build package
<Build_package>` doc page for more info.
This pair potential requires the :doc:`newtion <newton>` setting to be
This pair potential requires the :doc:`newton <newton>` setting to be
"on" for pair interactions.
The potential files provided with LAMMPS (see the potentials directory)

4
src/ASPHERE/compute_temp_asphere.cpp
View File

@ -58,9 +58,7 @@ ComputeTempAsphere::ComputeTempAsphere(LAMMPS *lmp, int narg, char **arg) :
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute temp/asphere command");
tempbias = 1;
int n = strlen(arg[iarg+1]) + 1;
id_bias = new char[n];
strcpy(id_bias,arg[iarg+1]);
id_bias = utils::strdup(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"dof") == 0) {
if (iarg+2 > narg)

17
src/ASPHERE/fix_nph_asphere.cpp
View File

@ -12,10 +12,9 @@
------------------------------------------------------------------------- */
#include "fix_nph_asphere.h"
#include <cstring>
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -36,21 +35,15 @@ FixNPHAsphere::FixNPHAsphere(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp,tcmd.c_str());
modify->add_compute(tcmd + " all temp/asphere");
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} all temp/asphere",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
std::string pcmd = id + std::string("_press");
id_press = new char[pcmd.size()+1];
strcpy(id_press,pcmd.c_str());
modify->add_compute(pcmd + " all pressure " + std::string(id_temp));
id_press = utils::strdup(std::string(id) + "_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
pcomputeflag = 1;
}

17
src/ASPHERE/fix_npt_asphere.cpp
View File

@ -12,10 +12,9 @@
------------------------------------------------------------------------- */
#include "fix_npt_asphere.h"
#include <cstring>
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -35,21 +34,15 @@ FixNPTAsphere::FixNPTAsphere(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp,tcmd.c_str());
modify->add_compute(tcmd + " all temp/asphere");
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} all temp/asphere",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
std::string pcmd = id + std::string("_press");
id_press = new char[pcmd.size()+1];
strcpy(id_press,pcmd.c_str());
modify->add_compute(pcmd + " all pressure " + std::string(id_temp));
id_press = utils::strdup(std::string(id) + "_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
pcomputeflag = 1;
}

3
src/ASPHERE/fix_nve_asphere.cpp
View File

@ -16,10 +16,11 @@
------------------------------------------------------------------------- */
#include "fix_nve_asphere.h"
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "error.h"
#include "math_extra.h"
using namespace LAMMPS_NS;
using namespace FixConst;

13
src/ASPHERE/fix_nvt_asphere.cpp
View File

@ -12,12 +12,10 @@
------------------------------------------------------------------------- */
#include "fix_nvt_asphere.h"
#include <cstring>
#include "error.h"
#include "group.h"
#include "modify.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -35,11 +33,8 @@ FixNVTAsphere::FixNVTAsphere(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp style
// id = fix-ID + temp
std::string cmd = id + std::string("_temp");
id_temp = new char[cmd.size()+1];
strcpy(id_temp,cmd.c_str());
cmd += fmt::format(" {} temp/asphere",group->names[igroup]);
modify->add_compute(cmd);
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp/asphere",
id_temp,group->names[igroup]));
tcomputeflag = 1;
}

4
src/BODY/compute_temp_body.cpp
View File

@ -56,9 +56,7 @@ ComputeTempBody::ComputeTempBody(LAMMPS *lmp, int narg, char **arg) :
if (iarg+2 > narg)
error->all(FLERR,"Illegal compute temp/body command");
tempbias = 1;
int n = strlen(arg[iarg+1]) + 1;
id_bias = new char[n];
strcpy(id_bias,arg[iarg+1]);
id_bias = utils::strdup(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"dof") == 0) {
if (iarg+2 > narg)

19
src/BODY/fix_nph_body.cpp
View File

@ -16,12 +16,12 @@
------------------------------------------------------------------------- */
#include "fix_nph_body.h"
#include <cstring>
#include "error.h"
#include "group.h"
#include "modify.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -41,22 +41,15 @@ FixNPHBody::FixNPHBody(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp,tcmd.c_str());
tcmd += fmt::format(" {} temp/body",group->names[igroup]);
modify->add_compute(tcmd);
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} all temp/body",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
std::string pcmd = id + std::string("_press");
id_press = new char[pcmd.size()+1];
strcpy(id_press,pcmd.c_str());
modify->add_compute(pcmd + " all pressure " + std::string(id_temp));
id_press = utils::strdup(std::string(id) + "_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
pcomputeflag = 1;
}

19
src/BODY/fix_npt_body.cpp
View File

@ -16,12 +16,12 @@
------------------------------------------------------------------------- */
#include "fix_npt_body.h"
#include <cstring>
#include "error.h"
#include "group.h"
#include "modify.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -41,22 +41,15 @@ FixNPTBody::FixNPTBody(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp,tcmd.c_str());
tcmd += fmt::format(" {} temp/body",group->names[igroup]);
modify->add_compute(tcmd);
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} all temp/body",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
std::string pcmd = id + std::string("_press");
id_press = new char[pcmd.size()+1];
strcpy(id_press,pcmd.c_str());
modify->add_compute(pcmd + " all pressure " + std::string(id_temp));
id_press = utils::strdup(std::string(id) + "_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
pcomputeflag = 1;
}

13
src/BODY/fix_nvt_body.cpp
View File

@ -16,12 +16,12 @@
------------------------------------------------------------------------- */
#include "fix_nvt_body.h"
#include <cstring>
#include "error.h"
#include "group.h"
#include "modify.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -39,11 +39,8 @@ FixNVTBody::FixNVTBody(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp style
// id = fix-ID + temp
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp,tcmd.c_str());
tcmd += fmt::format(" {} temp/body",group->names[igroup]);
modify->add_compute(tcmd);
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp/body",
id_temp,group->names[igroup]));
tcomputeflag = 1;
}

36
src/GPU/fix_npt_gpu.cpp
View File

@ -11,10 +11,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <cstring>
#include "fix_npt_gpu.h"
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -25,44 +25,24 @@ FixNPTGPU::FixNPTGPU(LAMMPS *lmp, int narg, char **arg) :
FixNHGPU(lmp, narg, arg)
{
if (!tstat_flag)
error->all(FLERR,"Temperature control must be used with fix npt/omp");
error->all(FLERR,"Temperature control must be used with fix npt/gpu");
if (!pstat_flag)
error->all(FLERR,"Pressure control must be used with fix npt/omp");
error->all(FLERR,"Pressure control must be used with fix npt/gpu");
// create a new compute temp style
// id = fix-ID + temp
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id)+"_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
id_press = utils::strdup(std::string(id)+"_press");
modify->add_compute(std::string(id_press)+" all pressure "+id_temp);
pcomputeflag = 1;
}

25
src/GPU/fix_nvt_gpu.cpp
View File

@ -11,11 +11,13 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <cstring>
#include "fix_nvt_gpu.h"
#include "error.h"
#include "group.h"
#include "modify.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -26,25 +28,14 @@ FixNVTGPU::FixNVTGPU(LAMMPS *lmp, int narg, char **arg) :
FixNHGPU(lmp, narg, arg)
{
if (!tstat_flag)
error->all(FLERR,"Temperature control must be used with fix nvt");
error->all(FLERR,"Temperature control must be used with fix nvt/gpu");
if (pstat_flag)
error->all(FLERR,"Pressure control can not be used with fix nvt");
error->all(FLERR,"Pressure control can not be used with fix nvt/gpu");
// create a new compute temp style
// id = fix-ID + temp
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp",id_temp,group->names[igroup]));
tcomputeflag = 1;
}

8
src/GRANULAR/fix_pour.cpp
View File

@ -948,18 +948,14 @@ void FixPour::options(int narg, char **arg)
} else if (strcmp(arg[iarg],"rigid") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
int n = strlen(arg[iarg+1]) + 1;
delete [] idrigid;
idrigid = new char[n];
strcpy(idrigid,arg[iarg+1]);
idrigid = utils::strdup(arg[iarg+1]);
rigidflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"shake") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix pour command");
int n = strlen(arg[iarg+1]) + 1;
delete [] idshake;
idshake = new char[n];
strcpy(idshake,arg[iarg+1]);
idshake = utils::strdup(arg[iarg+1]);
shakeflag = 1;
iarg += 2;

4
src/GRANULAR/fix_wall_gran.cpp
View File

@ -338,9 +338,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
} else if (strcmp(arg[iarg],"region") == 0) {
if (narg < iarg+2) error->all(FLERR,"Illegal fix wall/gran command");
wallstyle = REGION;
int n = strlen(arg[iarg+1]) + 1;
idregion = new char[n];
strcpy(idregion,arg[iarg+1]);
idregion = utils::strdup(arg[iarg+1]);
iarg += 2;
}

49
src/KIM/pair_kim.cpp
View File

@ -477,12 +477,10 @@ void PairKIM::coeff(int narg, char **arg)
if (paramname == str_name_str) break;
}
if (param_index >= numberOfParameters) {
auto msg = fmt::format("Wrong argument for pair coefficients.\n"
"This Model does not have the requested "
"'{}' parameter", paramname);
error->all(FLERR,msg);
}
if (param_index >= numberOfParameters)
error->all(FLERR,fmt::format("Wrong argument for pair coefficients.\n"
"This Model does not have the requested "
"'{}' parameter", paramname));
// Get the index_range for the requested parameter
int nlbound(0);
@ -492,12 +490,10 @@ void PairKIM::coeff(int narg, char **arg)
std::string argtostr(arg[i++]);
// Check to see if the indices range contains only integer numbers & :
if (argtostr.find_first_not_of("0123456789:") != std::string::npos) {
auto msg = fmt::format("Illegal index_range.\nExpected integer "
"parameter(s) instead of '{}' in "
"index_range", argtostr);
error->all(FLERR,msg);
}
if (argtostr.find_first_not_of("0123456789:") != std::string::npos)
error->all(FLERR,fmt::format("Illegal index_range.\nExpected integer"
" parameter(s) instead of '{}' in "
"index_range", argtostr));
std::string::size_type npos = argtostr.find(':');
if (npos != std::string::npos) {
@ -507,21 +503,17 @@ void PairKIM::coeff(int narg, char **arg)
nubound = atoi(words[1].c_str());
if (nubound < 1 || nubound > extent ||
nlbound < 1 || nlbound > nubound) {
auto msg = fmt::format("Illegal index_range '{}-{}' for '{}' "
"parameter with the extent of '{}'",
nlbound, nubound, paramname, extent);
error->all(FLERR,msg);
}
nlbound < 1 || nlbound > nubound)
error->all(FLERR,fmt::format("Illegal index_range '{}-{}' for '{}' "
"parameter with the extent of '{}'",
nlbound, nubound, paramname, extent));
} else {
nlbound = atoi(argtostr.c_str());
if (nlbound < 1 || nlbound > extent) {
auto msg = fmt::format("Illegal index '{}' for '{}' parameter "
"with the extent of '{}'", nlbound,
paramname, extent);
error->all(FLERR,msg);
}
if (nlbound < 1 || nlbound > extent)
error->all(FLERR,fmt::format("Illegal index '{}' for '{}' parameter "
"with the extent of '{}'", nlbound,
paramname, extent));
nubound = nlbound;
}
@ -551,11 +543,10 @@ void PairKIM::coeff(int narg, char **arg)
} else
error->all(FLERR,"Wrong parameter type to update");
} else {
auto msg = fmt::format("Wrong number of variable values for pair "
"coefficients.\n'{}' values are requested "
"for '{}' parameter", nubound - nlbound + 1,
paramname);
error->all(FLERR,msg);
error->all(FLERR,fmt::format("Wrong number of variable values for pair "
"coefficients.\n'{}' values are requested "
"for '{}' parameter", nubound - nlbound + 1,
paramname));
}
}

8
src/KOKKOS/atom_kokkos.cpp
View File

@ -257,9 +257,7 @@ int AtomKokkos::add_custom(const char *name, int flag)
nivector++;
iname = (char **) memory->srealloc(iname,nivector*sizeof(char *),
"atom:iname");
int n = strlen(name) + 1;
iname[index] = new char[n];
strcpy(iname[index],name);
iname[index] = utils::strdup(name);
ivector = (int **) memory->srealloc(ivector,nivector*sizeof(int *),
"atom:ivector");
memory->create(ivector[index],nmax,"atom:ivector");
@ -268,9 +266,7 @@ int AtomKokkos::add_custom(const char *name, int flag)
ndvector++;
dname = (char **) memory->srealloc(dname,ndvector*sizeof(char *),
"atom:dname");
int n = strlen(name) + 1;
dname[index] = new char[n];
strcpy(dname[index],name);
dname[index] = utils::strdup(name);
this->sync(Device,DVECTOR_MASK);
memoryKK->grow_kokkos(k_dvector,dvector,ndvector,nmax,
"atom:dvector");

9
src/KOKKOS/atom_vec_hybrid_kokkos.cpp
View File

@ -72,8 +72,7 @@ void AtomVecHybridKokkos::process_args(int narg, char **arg)
if (strcmp(arg[iarg],keywords[i]) == 0)
error->all(FLERR,"Atom style hybrid cannot use same atom style twice");
styles[nstyles] = atom->new_avec(arg[iarg],1,dummy);
keywords[nstyles] = new char[strlen(arg[iarg])+1];
strcpy(keywords[nstyles],arg[iarg]);
keywords[nstyles] = utils::strdup(arg[iarg]);
jarg = iarg + 1;
while (jarg < narg && !known_style(arg[jarg])) jarg++;
styles[nstyles]->process_args(jarg-iarg-1,&arg[iarg+1]);
@ -1175,10 +1174,8 @@ void AtomVecHybridKokkos::build_styles()
int n;
nallstyles = 0;
#define ATOM_CLASS
#define AtomStyle(key,Class) \
n = strlen(#key) + 1; \
allstyles[nallstyles] = new char[n]; \
strcpy(allstyles[nallstyles],#key); \
#define AtomStyle(key,Class) \
allstyles[nallstyles] = utils::strdup(#key); \
nallstyles++;
#include "style_atom.h" // IWYU pragma: keep
#undef AtomStyle

35
src/KOKKOS/fix_nph_kokkos.cpp
View File

@ -16,8 +16,6 @@
#include "modify.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -29,45 +27,26 @@ FixNPHKokkos<DeviceType>::FixNPHKokkos(LAMMPS *lmp, int narg, char **arg) :
{
this->kokkosable = 1;
if (this->tstat_flag)
this->error->all(FLERR,"Temperature control can not be used with fix nph");
this->error->all(FLERR,"Temperature control can not be used with fix nph/kk");
if (!this->pstat_flag)
this->error->all(FLERR,"Pressure control must be used with fix nph");
this->error->all(FLERR,"Pressure control must be used with fix nph/kk");
// create a new compute temp style
// id = fix-ID + temp
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(this->id) + 6;
this->id_temp = new char[n];
strcpy(this->id_temp,this->id);
strcat(this->id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = this->id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp/kk";
this->modify->add_compute(3,newarg);
delete [] newarg;
this->id_temp = utils::strdup(std::string(this->id) + "_temp");
this->modify->add_compute(fmt::format("{} all temp/kk",this->id_temp));
this->tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(this->id) + 7;
this->id_press = new char[n];
strcpy(this->id_press,this->id);
strcat(this->id_press,"_press");
newarg = new char*[4];
newarg[0] = this->id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = this->id_temp;
this->modify->add_compute(4,newarg);
delete [] newarg;
this->id_press = utils::strdup(std::string(this->id) + "_press");
this->modify->add_compute(fmt::format("{} all pressure {}",
this->id_press, this->id_temp));
this->pcomputeflag = 1;
}

35
src/KOKKOS/fix_npt_kokkos.cpp
View File

@ -16,8 +16,6 @@
#include "modify.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -29,45 +27,26 @@ FixNPTKokkos<DeviceType>::FixNPTKokkos(LAMMPS *lmp, int narg, char **arg) :
{
this->kokkosable = 1;
if (!this->tstat_flag)
this->error->all(FLERR,"Temperature control must be used with fix npt");
this->error->all(FLERR,"Temperature control must be used with fix npt/kk");
if (!this->pstat_flag)
this->error->all(FLERR,"Pressure control must be used with fix npt");
this->error->all(FLERR,"Pressure control must be used with fix npt/kk");
// create a new compute temp style
// id = fix-ID + temp
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(this->id) + 6;
this->id_temp = new char[n];
strcpy(this->id_temp,this->id);
strcat(this->id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = this->id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp/kk";
this->modify->add_compute(3,newarg);
delete [] newarg;
this->id_temp = utils::strdup(std::string(this->id) + "_temp");
this->modify->add_compute(fmt::format("{} all temp/kk",this->id_temp));
this->tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(this->id) + 7;
this->id_press = new char[n];
strcpy(this->id_press,this->id);
strcat(this->id_press,"_press");
newarg = new char*[4];
newarg[0] = this->id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = this->id_temp;
this->modify->add_compute(4,newarg);
delete [] newarg;
this->id_press = utils::strdup(std::string(this->id) + "_press");
this->modify->add_compute(fmt::format("{} all pressure {}",
this->id_press, this->id_temp));
this->pcomputeflag = 1;
}

18
src/KOKKOS/fix_nvt_kokkos.cpp
View File

@ -30,25 +30,15 @@ FixNVTKokkos<DeviceType>::FixNVTKokkos(LAMMPS *lmp, int narg, char **arg) :
{
this->kokkosable = 1;
if (!this->tstat_flag)
this->error->all(FLERR,"Temperature control must be used with fix nvt");
this->error->all(FLERR,"Temperature control must be used with fix nvt/kk");
if (this->pstat_flag)
this->error->all(FLERR,"Pressure control can not be used with fix nvt");
this->error->all(FLERR,"Pressure control can not be used with fix nvt/kk");
// create a new compute temp style
// id = fix-ID + temp
int n = strlen(this->id) + 6;
this->id_temp = new char[n];
strcpy(this->id_temp,this->id);
strcat(this->id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = this->id_temp;
newarg[1] = this->group->names[this->igroup];
newarg[2] = (char *) "temp/kk";
this->modify->add_compute(3,newarg);
delete [] newarg;
this->id_temp = utils::strdup(std::string(this->id)+"_temp");
this->modify->add_compute(fmt::format("{} all temp/kk",this->id_temp));
this->tcomputeflag = 1;
}

10
src/KOKKOS/pair_exp6_rx_kokkos.cpp
View File

@ -1800,14 +1800,8 @@ void PairExp6rxKokkos<DeviceType>::read_file(char *file)
}
params[nparams].ispecies = ispecies;
n = strlen(&atom->dname[ispecies][0]) + 1;
params[nparams].name = new char[n];
strcpy(params[nparams].name,&atom->dname[ispecies][0]);
n = strlen(words[1]) + 1;
params[nparams].potential = new char[n];
strcpy(params[nparams].potential,words[1]);
params[nparams].name = utils::strdup(&atom->dname[ispecies][0]);
params[nparams].potential = utils::strdup(words[1]);
if (strcmp(params[nparams].potential,"exp6") == 0) {
params[nparams].alpha = atof(words[2]);
params[nparams].epsilon = atof(words[3]);

18
src/KOKKOS/pair_table_rx_kokkos.cpp
View File

@ -1040,25 +1040,19 @@ void PairTableRXKokkos<DeviceType>::coeff(int narg, char **arg)
nspecies = atom->nspecies_dpd;
if (nspecies==0) error->all(FLERR,"There are no rx species specified.");
int n;
n = strlen(arg[4]) + 1;
site1 = new char[n];
strcpy(site1,arg[4]);
site1 = utils::strdup(arg[4]);
int ispecies;
for (ispecies = 0; ispecies < nspecies; ispecies++) {
for (ispecies = 0; ispecies < nspecies; ispecies++)
if (strcmp(site1,&atom->dname[ispecies][0]) == 0) break;
}
if (ispecies == nspecies && strcmp(site1,"1fluid") != 0)
error->all(FLERR,"Site1 name not recognized in pair coefficients");
n = strlen(arg[5]) + 1;
site2 = new char[n];
strcpy(site2,arg[5]);
for (ispecies = 0; ispecies < nspecies; ispecies++) {
site2 = utils::strdup(arg[5]);
for (ispecies = 0; ispecies < nspecies; ispecies++)
if (strcmp(site2,&atom->dname[ispecies][0]) == 0) break;
}
if (ispecies == nspecies && strcmp(site2,"1fluid") != 0)
error->all(FLERR,"Site2 name not recognized in pair coefficients");

5
src/LATTE/fix_latte.cpp
View File

@ -84,10 +84,7 @@ FixLatte::FixLatte(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR,"Fix latte does not yet support a LAMMPS calculation "
"of a Coulomb potential");
int n = strlen(arg[3]) + 1;
id_pe = new char[n];
strcpy(id_pe,arg[3]);
id_pe = utils::strdup(arg[3]);
int ipe = modify->find_compute(id_pe);
if (ipe < 0) error->all(FLERR,"Could not find fix latte compute ID");
if (modify->compute[ipe]->peatomflag == 0)

45
src/MC/fix_atom_swap.cpp
View File

@ -18,30 +18,31 @@
#include "fix_atom_swap.h"
#include "angle.h"
#include "atom.h"
#include "bond.h"
#include "comm.h"
#include "compute.h"
#include "dihedral.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "group.h"
#include "improper.h"
#include "kspace.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include "pair.h"
#include "random_park.h"
#include "region.h"
#include "update.h"
#include <cmath>
#include <cctype>
#include <cfloat>
#include <cstring>
#include "atom.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "comm.h"
#include "compute.h"
#include "group.h"
#include "domain.h"
#include "region.h"
#include "random_park.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "improper.h"
#include "kspace.h"
#include "memory.h"
#include "error.h"
#include "neighbor.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -138,9 +139,7 @@ void FixAtomSwap::options(int narg, char **arg)
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1)
error->all(FLERR,"Region ID for fix atom/swap does not exist");
int n = strlen(arg[iarg+1]) + 1;
idregion = new char[n];
strcpy(idregion,arg[iarg+1]);
idregion = utils::strdup(arg[iarg+1]);
regionflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"ke") == 0) {

41
src/MC/fix_bond_swap.cpp
View File

@ -13,27 +13,27 @@
#include "fix_bond_swap.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "neighbor.h"
#include "atom.h"
#include "bond.h"
#include "citeme.h"
#include "comm.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "comm.h"
#include "domain.h"
#include "modify.h"
#include "compute.h"
#include "neighbor.h"
#include "pair.h"
#include "random_mars.h"
#include "citeme.h"
#include "memory.h"
#include "error.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -86,11 +86,8 @@ FixBondSwap::FixBondSwap(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp style
// id = fix-ID + temp, compute group = fix group
std::string cmd = id + std::string("_temp");
id_temp = new char[cmd.size()+1];
strcpy(id_temp,cmd.c_str());
modify->add_compute(cmd + " all temp");
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tflag = 1;
// initialize atom list
@ -642,9 +639,7 @@ int FixBondSwap::modify_param(int narg, char **arg)
tflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
id_temp = utils::strdup(arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0)

53
src/MC/fix_gcmc.cpp
View File

@ -17,32 +17,33 @@
#include "fix_gcmc.h"
#include <cmath>
#include <cstring>
#include "angle.h"
#include "atom.h"
#include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "modify.h"
#include "fix.h"
#include "bond.h"
#include "comm.h"
#include "compute.h"
#include "group.h"
#include "domain.h"
#include "region.h"
#include "random_park.h"
#include "force.h"
#include "pair.h"
#include "bond.h"
#include "angle.h"
#include "dihedral.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "group.h"
#include "improper.h"
#include "kspace.h"
#include "math_extra.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "error.h"
#include "modify.h"
#include "molecule.h"
#include "neighbor.h"
#include "pair.h"
#include "random_park.h"
#include "region.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -301,9 +302,7 @@ void FixGCMC::options(int narg, char **arg)
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1)
error->all(FLERR,"Region ID for fix gcmc does not exist");
int n = strlen(arg[iarg+1]) + 1;
idregion = new char[n];
strcpy(idregion,arg[iarg+1]);
idregion = utils::strdup(arg[iarg+1]);
regionflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"maxangle") == 0) {
@ -327,18 +326,14 @@ void FixGCMC::options(int narg, char **arg)
iarg += 2;
} else if (strcmp(arg[iarg],"rigid") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
int n = strlen(arg[iarg+1]) + 1;
delete [] idrigid;
idrigid = new char[n];
strcpy(idrigid,arg[iarg+1]);
idrigid = utils::strdup(arg[iarg+1]);
rigidflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"shake") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix gcmc command");
int n = strlen(arg[iarg+1]) + 1;
delete [] idshake;
idshake = new char[n];
strcpy(idshake,arg[iarg+1]);
idshake = utils::strdup(arg[iarg+1]);
shakeflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"full_energy") == 0) {
@ -353,9 +348,7 @@ void FixGCMC::options(int narg, char **arg)
ngroupsmax*sizeof(char *),
"fix_gcmc:groupstrings");
}
int n = strlen(arg[iarg+1]) + 1;
groupstrings[ngroups] = new char[n];
strcpy(groupstrings[ngroups],arg[iarg+1]);
groupstrings[ngroups] = utils::strdup(arg[iarg+1]);
ngroups++;
iarg += 2;
} else if (strcmp(arg[iarg],"grouptype") == 0) {
@ -370,9 +363,7 @@ void FixGCMC::options(int narg, char **arg)
"fix_gcmc:grouptypestrings");
}
grouptypes[ngrouptypes] = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
int n = strlen(arg[iarg+2]) + 1;
grouptypestrings[ngrouptypes] = new char[n];
strcpy(grouptypestrings[ngrouptypes],arg[iarg+2]);
grouptypestrings[ngrouptypes] = utils::strdup(arg[iarg+2]);
ngrouptypes++;
iarg += 3;
} else if (strcmp(arg[iarg],"intra_energy") == 0) {

4
src/MC/fix_widom.cpp
View File

@ -224,9 +224,7 @@ void FixWidom::options(int narg, char **arg)
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1)
error->all(FLERR,"Region ID for fix widom does not exist");
int n = strlen(arg[iarg+1]) + 1;
idregion = new char[n];
strcpy(idregion,arg[iarg+1]);
idregion = utils::strdup(arg[iarg+1]);
regionflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"charge") == 0) {

26
src/MISC/compute_ti.cpp
View File

@ -17,17 +17,17 @@
#include "compute_ti.h"
#include <cstring>
#include "atom.h"
#include "update.h"
#include "domain.h"
#include "force.h"
#include "pair.h"
#include "kspace.h"
#include "input.h"
#include "variable.h"
#include "error.h"
#include "force.h"
#include "input.h"
#include "kspace.h"
#include "pair.h"
#include "update.h"
#include "variable.h"
#include <cstring>
using namespace LAMMPS_NS;
@ -77,21 +77,15 @@ ComputeTI::ComputeTI(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg],"tail") == 0) which[nterms] = TAIL;
else which[nterms] = PAIR;
int n = strlen(arg[iarg]) + 1;
pstyle[nterms] = new char[n];
strcpy(pstyle[nterms],arg[iarg]);
pstyle[nterms] = utils::strdup(arg[iarg]);
utils::bounds(FLERR,arg[iarg+1],1,atom->ntypes,ilo[nterms],ihi[nterms],error);
iarg += 1;
if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) {
int n = strlen(&arg[iarg+1][2]) + 1;
var1[nterms] = new char[n];
strcpy(var1[nterms],&arg[iarg+1][2]);
var1[nterms] = utils::strdup(&arg[iarg+1][2]);
} else error->all(FLERR,"Illegal compute ti command");
if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) {
int n = strlen(&arg[iarg+2][2]) + 1;
var2[nterms] = new char[n];
strcpy(var2[nterms],&arg[iarg+2][2]);
var2[nterms] = utils::strdup(&arg[iarg+2][2]);
} else error->all(FLERR,"Illegal compute ti command");
nterms++;

12
src/MISC/fix_deposit.cpp
View File

@ -697,9 +697,7 @@ void FixDeposit::options(int narg, char **arg)
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1)
error->all(FLERR,"Region ID for fix deposit does not exist");
int n = strlen(arg[iarg+1]) + 1;
idregion = new char[n];
strcpy(idregion,arg[iarg+1]);
idregion = utils::strdup(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"mol") == 0) {
@ -728,18 +726,14 @@ void FixDeposit::options(int narg, char **arg)
iarg += nmol+1;
} else if (strcmp(arg[iarg],"rigid") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
int n = strlen(arg[iarg+1]) + 1;
delete [] idrigid;
idrigid = new char[n];
strcpy(idrigid,arg[iarg+1]);
idrigid = utils::strdup(arg[iarg+1]);
rigidflag = 1;
iarg += 2;
} else if (strcmp(arg[iarg],"shake") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix deposit command");
int n = strlen(arg[iarg+1]) + 1;
delete [] idshake;
idshake = new char[n];
strcpy(idshake,arg[iarg+1]);
idshake = utils::strdup(arg[iarg+1]);
shakeflag = 1;
iarg += 2;

4
src/MISC/fix_evaporate.cpp
View File

@ -44,9 +44,7 @@ FixEvaporate::FixEvaporate(LAMMPS *lmp, int narg, char **arg) :
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
nflux = utils::inumeric(FLERR,arg[4],false,lmp);
iregion = domain->find_region(arg[5]);
int n = strlen(arg[5]) + 1;
idregion = new char[n];
strcpy(idregion,arg[5]);
idregion = utils::strdup(arg[5]);
int seed = utils::inumeric(FLERR,arg[6],false,lmp);
if (nevery <= 0 || nflux <= 0)

4
src/MISC/fix_oneway.cpp
View File

@ -42,9 +42,7 @@ FixOneWay::FixOneWay(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
nevery = utils::inumeric(FLERR,arg[3],false,lmp);
if (nevery < 1) error->all(FLERR,"Illegal fix oneway command");
int len = strlen(arg[4]);
regionstr = new char[len+1];
strcpy(regionstr,arg[4]);
regionstr = utils::strdup(arg[4]);
if (strcmp(arg[5], "x") == 0) direction = X|PLUS;
if (strcmp(arg[5], "X") == 0) direction = X|PLUS;

16
src/MISC/fix_orient_bcc.cpp
View File

@ -82,19 +82,11 @@ FixOrientBCC::FixOrientBCC(LAMMPS *lmp, int narg, char **arg) :
uxif_high = utils::numeric(FLERR,arg[8],false,lmp);
if (direction_of_motion == 0) {
int n = strlen(arg[9]) + 1;
chifilename = new char[n];
strcpy(chifilename,arg[9]);
n = strlen(arg[10]) + 1;
xifilename = new char[n];
strcpy(xifilename,arg[10]);
chifilename = utils::strdup(arg[9]);
xifilename = utils::strdup(arg[10]);
} else if (direction_of_motion == 1) {
int n = strlen(arg[9]) + 1;
xifilename = new char[n];
strcpy(xifilename,arg[9]);
n = strlen(arg[10]) + 1;
chifilename = new char[n];
strcpy(chifilename,arg[10]);
xifilename = utils::strdup(arg[9]);
chifilename = utils::strdup(arg[10]);
} else error->all(FLERR,"Illegal fix orient/bcc command");
// initializations

16
src/MISC/fix_orient_fcc.cpp
View File

@ -80,19 +80,11 @@ FixOrientFCC::FixOrientFCC(LAMMPS *lmp, int narg, char **arg) :
uxif_high = utils::numeric(FLERR,arg[8],false,lmp);
if (direction_of_motion == 0) {
int n = strlen(arg[9]) + 1;
chifilename = new char[n];
strcpy(chifilename,arg[9]);
n = strlen(arg[10]) + 1;
xifilename = new char[n];
strcpy(xifilename,arg[10]);
chifilename = utils::strdup(arg[9]);
xifilename = utils::strdup(arg[10]);
} else if (direction_of_motion == 1) {
int n = strlen(arg[9]) + 1;
xifilename = new char[n];
strcpy(xifilename,arg[9]);
n = strlen(arg[10]) + 1;
chifilename = new char[n];
strcpy(chifilename,arg[10]);
xifilename = utils::strdup(arg[9]);
chifilename = utils::strdup(arg[10]);
} else error->all(FLERR,"Illegal fix orient/fcc command");
// initializations

9
src/MLIAP/mliap_descriptor_snap.cpp
View File

@ -427,22 +427,19 @@ void MLIAPDescriptorSNAP::read_paramfile(char *paramfilename)
if (strcmp(keywd,"elems") == 0) {
for (int ielem = 0; ielem < nelements; ielem++) {
char* elemtmp = keyval;
int n = strlen(elemtmp) + 1;
elements[ielem] = new char[n];
strcpy(elements[ielem],elemtmp);
elements[ielem] = utils::strdup(keyval);
keyval = strtok(nullptr,"' \t\n\r\f");
}
elementsflag = 1;
} else if (strcmp(keywd,"radelems") == 0) {
for (int ielem = 0; ielem < nelements; ielem++) {
radelem[ielem] = atof(keyval);
radelem[ielem] = utils::numeric(FLERR,keyval,false,lmp);
keyval = strtok(nullptr,"' \t\n\r\f");
}
radelemflag = 1;
} else if (strcmp(keywd,"welems") == 0) {
for (int ielem = 0; ielem < nelements; ielem++) {
wjelem[ielem] = atof(keyval);
wjelem[ielem] = utils::numeric(FLERR,keyval,false,lmp);
keyval = strtok(nullptr,"' \t\n\r\f");
}
wjelemflag = 1;

2
src/MLIAP/mliap_model_nn.cpp
View File

@ -18,9 +18,9 @@
#include "mliap_model_nn.h"
#include "pair_mliap.h"
#include "mliap_data.h"
#include "error.h"
#include "comm.h"
#include "error.h"
#include "memory.h"
#include "tokenizer.h"

4
src/MOLECULE/angle_cosine_periodic.cpp
View File

@ -222,9 +222,9 @@ void AngleCosinePeriodic::coeff(int narg, char **arg)
/* ---------------------------------------------------------------------- */
double AngleCosinePeriodic::equilibrium_angle(int /*i*/)
double AngleCosinePeriodic::equilibrium_angle(int i)
{
return MY_PI;
return MY_PI*(1.0 - (b[i]>0)?0.0:1.0/static_cast<double>(multiplicity[i]));
}
/* ----------------------------------------------------------------------

4
src/REPLICA/compute_event_displace.cpp
View File

@ -206,7 +206,5 @@ void ComputeEventDisplace::reset_extra_compute_fix(const char *id_new)
id_event = nullptr;
if (id_new == nullptr) return;
int n = strlen(id_new) + 1;
id_event = new char[n];
strcpy(id_event,id_new);
id_event = utils::strdup(id_new);
}

22
src/REPLICA/prd.cpp
View File

@ -77,8 +77,7 @@ void PRD::command(int narg, char **arg)
t_dephase = utils::inumeric(FLERR,arg[3],false,lmp);
t_corr = utils::inumeric(FLERR,arg[4],false,lmp);
char *id_compute = new char[strlen(arg[5])+1];
strcpy(id_compute,arg[5]);
char *id_compute = utils::strdup(arg[5]);
int seed = utils::inumeric(FLERR,arg[6],false,lmp);
options(narg-7,&arg[7]);
@ -875,15 +874,8 @@ void PRD::options(int narg, char **arg)
temp_flag = 0;
stepmode = 0;
char *str = (char *) "geom";
int n = strlen(str) + 1;
loop_setting = new char[n];
strcpy(loop_setting,str);
str = (char *) "gaussian";
n = strlen(str) + 1;
dist_setting = new char[n];
strcpy(dist_setting,str);
loop_setting = utils::strdup("geom");
dist_setting = utils::strdup("gaussian");
int iarg = 0;
while (iarg < narg) {
@ -912,16 +904,12 @@ void PRD::options(int narg, char **arg)
else if (strcmp(arg[iarg+1],"local") == 0) loop_setting = nullptr;
else if (strcmp(arg[iarg+1],"geom") == 0) loop_setting = nullptr;
else error->all(FLERR,"Illegal prd command");
int n = strlen(arg[iarg+1]) + 1;
loop_setting = new char[n];
strcpy(loop_setting,arg[iarg+1]);
loop_setting = utils::strdup(arg[iarg+1]);
if (strcmp(arg[iarg+2],"uniform") == 0) dist_setting = nullptr;
else if (strcmp(arg[iarg+2],"gaussian") == 0) dist_setting = nullptr;
else error->all(FLERR,"Illegal prd command");
n = strlen(arg[iarg+2]) + 1;
dist_setting = new char[n];
strcpy(dist_setting,arg[iarg+2]);
dist_setting = utils::strdup(arg[iarg+2]);
iarg += 3;

19
src/REPLICA/tad.cpp
View File

@ -89,15 +89,12 @@ void TAD::command(int narg, char **arg)
delta_conf = utils::numeric(FLERR,arg[4],false,lmp);
tmax = utils::numeric(FLERR,arg[5],false,lmp);
char *id_compute = new char[strlen(arg[6])+1];
strcpy(id_compute,arg[6]);
char *id_compute = utils::strdup(arg[6]);
// quench minimizer is set by min_style command
// NEB minimizer is set by options, default = quickmin
int n = strlen(update->minimize_style) + 1;
min_style = new char[n];
strcpy(min_style,update->minimize_style);
min_style = utils::strdup(update->minimize_style);
options(narg-7,&arg[7]);
@ -585,9 +582,7 @@ void TAD::options(int narg, char **arg)
n2steps_neb = 100;
nevery_neb = 10;
int n = strlen("quickmin") + 1;
min_style_neb = new char[n];
strcpy(min_style_neb,"quickmin");
min_style_neb = utils::strdup("quickmin");
dt_neb = update->dt;
neb_logfilename = nullptr;
@ -619,9 +614,7 @@ void TAD::options(int narg, char **arg)
} else if (strcmp(arg[iarg],"neb_style") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal tad command");
delete [] min_style_neb;
int n = strlen(arg[iarg+1]) + 1;
min_style_neb = new char[n];
strcpy(min_style_neb,arg[iarg+1]);
min_style_neb = utils::strdup(arg[iarg+1]);
iarg += 2;
} else if (strcmp(arg[iarg],"neb_step") == 0) {
@ -635,9 +628,7 @@ void TAD::options(int narg, char **arg)
if (iarg+2 > narg) error->all(FLERR,"Illegal tad command");
if (strcmp(arg[iarg+1],"none") == 0) neb_logfilename = nullptr;
else {
int n = strlen(arg[iarg+1]) + 1;
neb_logfilename = new char[n];
strcpy(neb_logfilename,arg[iarg+1]);
neb_logfilename = utils::strdup(arg[iarg+1]);
}
iarg += 2;
} else error->all(FLERR,"Illegal tad command");

16
src/RIGID/compute_erotate_rigid.cpp
View File

@ -12,14 +12,16 @@
------------------------------------------------------------------------- */
#include "compute_erotate_rigid.h"
#include <cstring>
#include "update.h"
#include "force.h"
#include "modify.h"
#include "error.h"
#include "fix.h"
#include "fix_rigid.h"
#include "fix_rigid_small.h"
#include "error.h"
#include "force.h"
#include "modify.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
@ -33,9 +35,7 @@ ComputeERotateRigid::ComputeERotateRigid(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
extscalar = 1;
int n = strlen(arg[3]) + 1;
rfix = new char[n];
strcpy(rfix,arg[3]);
rfix = utils::strdup(arg[3]);
}
/* ---------------------------------------------------------------------- */

16
src/RIGID/compute_ke_rigid.cpp
View File

@ -12,14 +12,16 @@
------------------------------------------------------------------------- */
#include "compute_ke_rigid.h"
#include <cstring>
#include "update.h"
#include "force.h"
#include "modify.h"
#include "error.h"
#include "fix.h"
#include "fix_rigid.h"
#include "fix_rigid_small.h"
#include "error.h"
#include "force.h"
#include "modify.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
@ -33,9 +35,7 @@ ComputeKERigid::ComputeKERigid(LAMMPS *lmp, int narg, char **arg) :
scalar_flag = 1;
extscalar = 1;
int n = strlen(arg[3]) + 1;
rfix = new char[n];
strcpy(rfix,arg[3]);
rfix = utils::strdup(arg[3]);
}
/* ---------------------------------------------------------------------- */

4
src/RIGID/compute_rigid_local.cpp
View File

@ -40,9 +40,7 @@ ComputeRigidLocal::ComputeRigidLocal(LAMMPS *lmp, int narg, char **arg) :
local_flag = 1;
nvalues = narg - 4;
int n = strlen(arg[3]) + 1;
idrigid = new char[n];
strcpy(idrigid,arg[3]);
idrigid = utils::strdup(arg[3]);
rstyle = new int[nvalues];

22
src/RIGID/fix_ehex.cpp
View File

@ -24,17 +24,19 @@
#include "fix_ehex.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "domain.h"
#include "region.h"
#include "group.h"
#include "force.h"
#include "update.h"
#include "modify.h"
#include "memory.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "memory.h"
#include "modify.h"
#include "region.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include "fix_shake.h"
using namespace LAMMPS_NS;
@ -91,9 +93,7 @@ FixEHEX::FixEHEX(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
iregion = domain->find_region(arg[iarg+1]);
if (iregion == -1)
error->all(FLERR,"Region ID for fix ehex does not exist");
int n = strlen(arg[iarg+1]) + 1;
idregion = new char[n];
strcpy(idregion,arg[iarg+1]);
idregion = utils::strdup(arg[iarg+1]);
iarg += 2;
}

8
src/RIGID/fix_rigid_nh.cpp
View File

@ -1202,9 +1202,7 @@ int FixRigidNH::modify_param(int narg, char **arg)
tcomputeflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
id_temp = utils::strdup(arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0)
@ -1236,9 +1234,7 @@ int FixRigidNH::modify_param(int narg, char **arg)
pcomputeflag = 0;
}
delete [] id_press;
int n = strlen(arg[1]) + 1;
id_press = new char[n];
strcpy(id_press,arg[1]);
id_press = utils::strdup(arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");

8
src/RIGID/fix_rigid_nh_small.cpp
View File

@ -1319,9 +1319,7 @@ int FixRigidNHSmall::modify_param(int narg, char **arg)
tcomputeflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
id_temp = utils::strdup(arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0)
@ -1353,9 +1351,7 @@ int FixRigidNHSmall::modify_param(int narg, char **arg)
pcomputeflag = 0;
}
delete [] id_press;
int n = strlen(arg[1]) + 1;
id_press = new char[n];
strcpy(id_press,arg[1]);
id_press = utils::strdup(arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");

31
src/RIGID/fix_rigid_nph.cpp
View File

@ -18,9 +18,9 @@
------------------------------------------------------------------------- */
#include "fix_rigid_nph.h"
#include <cstring>
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
@ -55,34 +55,15 @@ FixRigidNPH::FixRigidNPH(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id)+"_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
id_press = utils::strdup(std::string(id)+"_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press,id_temp));
pcomputeflag = 1;
}

31
src/RIGID/fix_rigid_nph_small.cpp
View File

@ -18,9 +18,9 @@
------------------------------------------------------------------------- */
#include "fix_rigid_nph_small.h"
#include <cstring>
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
@ -58,34 +58,15 @@ FixRigidNPHSmall::FixRigidNPHSmall(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id)+"_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
id_press = utils::strdup(std::string(id)+"_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press,id_temp));
pcomputeflag = 1;
}

31
src/RIGID/fix_rigid_npt.cpp
View File

@ -18,9 +18,9 @@
------------------------------------------------------------------------- */
#include "fix_rigid_npt.h"
#include <cstring>
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
@ -65,34 +65,15 @@ FixRigidNPT::FixRigidNPT(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id)+"_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
id_press = utils::strdup(std::string(id)+"_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press,id_temp));
pcomputeflag = 1;
}

31
src/RIGID/fix_rigid_npt_small.cpp
View File

@ -18,9 +18,9 @@
------------------------------------------------------------------------- */
#include "fix_rigid_npt_small.h"
#include <cstring>
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
@ -69,34 +69,15 @@ FixRigidNPTSmall::FixRigidNPTSmall(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id)+"_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
id_press = utils::strdup(std::string(id)+"_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press,id_temp));
pcomputeflag = 1;
}

25
src/SHOCK/fix_msst.cpp
View File

@ -19,19 +19,20 @@
#include "fix_msst.h"
#include <cstring>
#include <cmath>
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "modify.h"
#include "fix_external.h"
#include "compute.h"
#include "kspace.h"
#include "update.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#include "fix_external.h"
#include "force.h"
#include "kspace.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -840,9 +841,7 @@ int FixMSST::modify_param(int narg, char **arg)
tflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
id_temp = utils::strdup(arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
@ -864,9 +863,7 @@ int FixMSST::modify_param(int narg, char **arg)
pflag = 0;
}
delete [] id_press;
int n = strlen(arg[1]) + 1;
id_press = new char[n];
strcpy(id_press,arg[1]);
id_press = utils::strdup(arg[1]);
int icompute = modify->find_compute(id_press);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");

55
src/SHOCK/fix_nphug.cpp
View File

@ -12,15 +12,17 @@
------------------------------------------------------------------------- */
#include "fix_nphug.h"
#include <cstring>
#include <cmath>
#include "modify.h"
#include "error.h"
#include "update.h"
#include "compute.h"
#include "force.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "modify.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -118,51 +120,22 @@ FixNPHug::FixNPHug(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id)+"_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
id_press = utils::strdup(std::string(id)+"_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press,id_temp));
pcomputeflag = 1;
// create a new compute potential energy compute
n = strlen(id) + 4;
id_pe = new char[n];
strcpy(id_pe,id);
strcat(id_pe,"_pe");
newarg = new char*[3];
newarg[0] = id_pe;
newarg[1] = (char*)"all";
newarg[2] = (char*)"pe";
modify->add_compute(3,newarg);
delete [] newarg;
id_pe = utils::strdup(std::string(id)+"_pe");
modify->add_compute(fmt::format("{} all pe",id_pe));
peflag = 1;
}

28
src/SPIN/fix_neb_spin.cpp
View File

@ -23,18 +23,19 @@
#include "fix_neb_spin.h"
#include <cmath>
#include <cstring>
#include "universe.h"
#include "update.h"
#include "atom.h"
#include "comm.h"
#include "modify.h"
#include "compute.h"
#include "group.h"
#include "memory.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "memory.h"
#include "modify.h"
#include "universe.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -120,17 +121,8 @@ FixNEBSpin::FixNEBSpin(LAMMPS *lmp, int narg, char **arg) :
// create a new compute pe style
// id = fix-ID + pe, compute group = all
int n = strlen(id) + 4;
id_pe = new char[n];
strcpy(id_pe,id);
strcat(id_pe,"_pe");
char **newarg = new char*[3];
newarg[0] = id_pe;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pe";
modify->add_compute(3,newarg);
delete [] newarg;
id_pe = utils::strdup(std::string(id)+"_pe");
modify->add_compute(std::string(id_pe)+" all pe");
// initialize local storage

8
src/USER-BOCS/compute_pressure_bocs.cpp
View File

@ -59,9 +59,7 @@ ComputePressureBocs::ComputePressureBocs(LAMMPS *lmp, int narg, char **arg) :
if (strcmp(arg[3],"NULL") == 0) id_temp = nullptr;
else {
int n = strlen(arg[3]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[3]);
id_temp = utils::strdup(arg[3]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0)
@ -446,7 +444,5 @@ void ComputePressureBocs::virial_compute(int n, int ndiag)
void ComputePressureBocs::reset_extra_compute_fix(const char *id_new)
{
delete [] id_temp;
int n = strlen(id_new) + 1;
id_temp = new char[n];
strcpy(id_temp,id_new);
id_temp = utils::strdup(id_new);
}

51
src/USER-BOCS/fix_bocs.cpp
View File

@ -16,11 +16,6 @@
#include "fix_bocs.h"
#include <cstring>
#include <cmath>
#include <vector>
#include "atom.h"
#include "citeme.h"
#include "comm.h"
@ -29,7 +24,6 @@
#include "domain.h"
#include "error.h"
#include "fix_deform.h"
#include "force.h"
#include "group.h"
#include "irregular.h"
@ -40,6 +34,10 @@
#include "respa.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include <vector>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -70,7 +68,8 @@ const int NUM_INPUT_DATA_COLUMNS = 2; // columns in the pressure correction
FixBocs::FixBocs(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
rfix(nullptr), id_dilate(nullptr), irregular(nullptr), id_temp(nullptr), id_press(nullptr),
rfix(nullptr), id_dilate(nullptr), irregular(nullptr),
id_temp(nullptr), id_press(nullptr),
eta(nullptr), eta_dot(nullptr), eta_dotdot(nullptr),
eta_mass(nullptr), etap(nullptr), etap_dot(nullptr), etap_dotdot(nullptr),
etap_mass(nullptr)
@ -403,38 +402,16 @@ FixBocs::FixBocs(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id)+"_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "PRESSURE/BOCS";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
id_press = utils::strdup(std::string(id)+"_press");
modify->add_compute(fmt::format("{} all PRESSURE/BOCS {}",id_press,id_temp));
pcomputeflag = 1;
/*~ MRD End of stuff copied from fix_npt.cpp~*/
@ -1548,9 +1525,7 @@ int FixBocs::modify_param(int narg, char **arg)
tcomputeflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
id_temp = utils::strdup(arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0)
@ -1582,9 +1557,7 @@ int FixBocs::modify_param(int narg, char **arg)
pcomputeflag = 0;
}
delete [] id_press;
int n = strlen(arg[1]) + 1;
id_press = new char[n];
strcpy(id_press,arg[1]);
id_press = utils::strdup(arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");

189
src/USER-COLVARS/group_ndx.cpp
View File

@ -23,6 +23,8 @@
#include "error.h"
#include "group.h"
#include <cmath>
using namespace LAMMPS_NS;
/* ----------------------------------------------------------------------
@ -40,96 +42,6 @@ static int cmptagint(const void *p1, const void *p2)
}
}
/* ----------------------------------------------------------------------
helper function. writes out one group to a gromacs style index file
---------------------------------------------------------------------- */
static void write_group(FILE *fp, int gid, Atom *atom, Group *group, int me,
int np, MPI_Comm world, FILE *screen, FILE *logfile)
{
char fmt[16];
tagint *sendlist, *recvlist;
bigint num = group->count(gid);
int lnum, cols;
if (me == 0) {
if (screen) fprintf(screen, " writing group %s... ", group->names[gid]);
if (logfile) fprintf(logfile, " writing group %s... ", group->names[gid]);
// the "all" group in LAMMPS is called "System" in gromacs
if (gid == 0) {
fputs("[ System ]\n", fp);
} else {
fprintf(fp,"[ %s ]\n", group->names[gid]);
}
// derive format string for index lists
bigint j = atom->natoms;
int i=0;
while (j > 0) {
++i;
j /= 10;
}
snprintf(fmt,16,"%%%dd ", i);
cols = 80 / (i+1);
}
if (num > 0) {
const int * const mask = atom->mask;
const tagint * const tag = atom->tag;
const int groupbit = group->bitmask[gid];
const int nlocal = atom->nlocal;
int i;
sendlist = new tagint[nlocal];
recvlist = new tagint[num];
lnum = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) sendlist[lnum++] = tag[i];
int nrecv,allrecv;
if (me == 0) {
MPI_Status status;
MPI_Request request;
for (i=0; i < lnum; i++)
recvlist[i] = sendlist[i];
allrecv = lnum;
for (i=1; i < np; ++i) {
MPI_Irecv(recvlist+allrecv,num-allrecv,MPI_LMP_TAGINT,i,0, world,&request);
MPI_Send(&nrecv,0,MPI_INT,i,0,world);
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_LMP_TAGINT,&nrecv);
allrecv += nrecv;
}
// sort received list
qsort((void *)recvlist, num, sizeof(tagint), cmptagint);
} else {
MPI_Recv(&nrecv,0,MPI_INT,0,0,world,MPI_STATUS_IGNORE);
MPI_Rsend(sendlist,lnum,MPI_LMP_TAGINT,0,0,world);
}
delete [] sendlist;
}
if (me == 0) {
int i, j;
for (i=0, j=0; i < num; ++i) {
fprintf(fp,fmt,recvlist[i]);
++j;
if (j == cols) {
fputs("\n",fp);
j = 0;
}
}
if (j > 0) fputs("\n",fp);
if (screen) fputs("done\n",screen);
if (logfile) fputs("done\n",logfile);
}
if (num > 0) delete[] recvlist;
}
/* ---------------------------------------------------------------------- */
void Group2Ndx::command(int narg, char **arg)
@ -144,31 +56,100 @@ void Group2Ndx::command(int narg, char **arg)
if (comm->me == 0) {
fp = fopen(arg[0], "w");
if (fp == nullptr)
error->one(FLERR,"Cannot open index file for writing");
if (screen)
fprintf(screen, "Writing groups to index file %s:\n",arg[0]);
if (logfile)
fprintf(logfile,"Writing groups to index file %s:\n",arg[0]);
error->one(FLERR,fmt::format("Cannot open index file for writing: {}",
utils::getsyserror()));
utils::logmesg(lmp,fmt::format("Writing groups to index file {}:\n",arg[0]));
}
if (narg == 1) { // write out all groups
for (int i=0; i < group->ngroup; ++i) {
write_group(fp,i,atom,group,comm->me,comm->nprocs,world,screen,logfile);
write_group(fp,i);
}
} else { // write only selected groups
for (int i=1; i < narg; ++i) {
int gid = group->find(arg[i]);
if (gid < 0) error->all(FLERR, "Non-existing group requested");
write_group(fp,gid,atom,group,comm->me,comm->nprocs,world,screen,logfile);
write_group(fp,gid);
}
}
if (comm->me == 0) {
if (screen) fputs("\n",screen);
if (logfile) fputs("\n",logfile);
fclose(fp);
}
if (comm->me == 0) fclose(fp);
}
/* ----------------------------------------------------------------------
write out one group to a Gromacs style index file
---------------------------------------------------------------------- */
void Group2Ndx::write_group(FILE *fp, int gid)
{
tagint *sendlist, *recvlist;
bigint gcount = group->count(gid);
int lnum, width, cols;
if (comm->me == 0) {
utils::logmesg(lmp,fmt::format(" writing group {}...",group->names[gid]));
// the "all" group in LAMMPS is called "System" in Gromacs
if (gid == 0) {
fputs("[ System ]\n", fp);
} else {
fmt::print(fp,"[ {} ]\n", group->names[gid]);
}
width = log10((double) atom->natoms)+2;
cols = 80 / width;
}
if (gcount > 0) {
const int * const mask = atom->mask;
const tagint * const tag = atom->tag;
const int groupbit = group->bitmask[gid];
const int nlocal = atom->nlocal;
int i;
sendlist = new tagint[nlocal];
recvlist = new tagint[gcount];
lnum = 0;
for (i = 0; i < nlocal; i++)
if (mask[i] & groupbit) sendlist[lnum++] = tag[i];
int nrecv=0;
bigint allrecv;
if (comm->me == 0) {
MPI_Status status;
MPI_Request request;
for (i=0; i < lnum; i++)
recvlist[i] = sendlist[i];
allrecv = lnum;
for (i=1; i < comm->nprocs; ++i) {
MPI_Irecv(recvlist+allrecv,gcount-allrecv,MPI_LMP_TAGINT,i,0, world,&request);
MPI_Send(&nrecv,0,MPI_INT,i,0,world); // block rank "i" until we are ready to receive
MPI_Wait(&request,&status);
MPI_Get_count(&status,MPI_LMP_TAGINT,&nrecv);
allrecv += nrecv;
}
// sort received list
qsort((void *)recvlist, allrecv, sizeof(tagint), cmptagint);
} else {
MPI_Recv(&nrecv,0,MPI_INT,0,0,world,MPI_STATUS_IGNORE);
MPI_Rsend(sendlist,lnum,MPI_LMP_TAGINT,0,0,world);
}
delete [] sendlist;
}
if (comm->me == 0) {
int i, j;
for (i=0, j=0; i < gcount; ++i) {
fmt::print(fp,"{:>{}}",recvlist[i],width);
++j;
if (j == cols) {
fputs("\n",fp);
j = 0;
}
}
if (j > 0) fputs("\n",fp);
utils::logmesg(lmp,"done\n");
}
if (gcount > 0) delete[] recvlist;
}

2
src/USER-COLVARS/group_ndx.h
View File

@ -30,6 +30,8 @@ class Group2Ndx : protected Pointers {
public:
Group2Ndx(class LAMMPS *lmp) : Pointers(lmp) {};
void command(int, char **);
private:
void write_group(FILE *, int);
};
}

233
src/USER-COLVARS/ndx_group.cpp
View File

@ -22,64 +22,52 @@
#include "comm.h"
#include "error.h"
#include "group.h"
#include "memory.h"
#include "tokenizer.h"
#include <cstring>
#include <cstdlib>
using namespace LAMMPS_NS;
#define BUFLEN 4096
#define DELTA 16384
static char *find_section(FILE *fp, const char *name)
// read file until next section "name" or any next section if name == ""
static std::string find_section(FILE *fp, const std::string &name)
{
char linebuf[BUFLEN];
char *n,*p,*t,*r;
while ((p = fgets(linebuf,BUFLEN,fp))) {
t = strtok(p," \t\n\r\f");
if ((t != nullptr) && *t == '[') {
t = strtok(nullptr," \t\n\r\f");
if (t != nullptr) {
n = t;
t = strtok(nullptr," \t\n\r\f");
if ((t != nullptr) && *t == ']') {
if ((name == nullptr) || strcmp(name,n) == 0) {
int l = strlen(n);
r = new char[l+1];
strncpy(r,n,l+1);
return r;
}
}
}
}
std::string pattern = "^\\s*\\[\\s+\\S+\\s+\\]\\s*$";
if (!name.empty())
pattern = fmt::format("^\\s*\\[\\s+{}\\s+\\]\\s*$",name);
fgets(linebuf,BUFLEN,fp);
while (!feof(fp)) {
if (utils::strmatch(linebuf,pattern))
return Tokenizer(linebuf).as_vector()[1];
fgets(linebuf,BUFLEN,fp);
}
return nullptr;
return "";
}
static tagint *read_section(FILE *fp, bigint &num)
static std::vector<tagint> read_section(FILE *fp, std::string &name)
{
char linebuf[BUFLEN];
char *p,*t;
tagint *tagbuf;
bigint nmax;
std::vector<tagint> tagbuf;
std::string pattern = "^\\s*\\[\\s+\\S+\\s+\\]\\s*$";
num = 0;
nmax = DELTA;
tagbuf = (tagint *)malloc(sizeof(tagint)*nmax);
while ((p = fgets(linebuf,BUFLEN,fp))) {
t = strtok(p," \t\n\r\f");
while (t != nullptr) {
// start of a new section. we are done here.
if (*t == '[') return tagbuf;
tagbuf[num++] = ATOTAGINT(t);
if (num == nmax) {
nmax += DELTA;
tagbuf = (tagint *)realloc(tagbuf,sizeof(tagint)*nmax);
}
t = strtok(nullptr," \t\n\r\f");
while (fgets(linebuf,BUFLEN,fp)) {
// start of new section. we are done, update "name"
if (utils::strmatch(linebuf,pattern)) {
name = Tokenizer(linebuf).as_vector()[1];
return tagbuf;
}
ValueTokenizer values(linebuf);
while (values.has_next())
tagbuf.push_back(values.next_tagint());
}
// set empty name to indicate end of file
name = "";
return tagbuf;
}
@ -90,151 +78,122 @@ void Ndx2Group::command(int narg, char **arg)
int len;
bigint num;
FILE *fp;
char *name = nullptr;
tagint *tags;
tagint *tagbuf;
std::string name = "", next;
if (narg < 1) error->all(FLERR,"Illegal ndx2group command");
if (atom->tag_enable == 0)
error->all(FLERR,"Must have atom IDs for ndx2group command");
error->all(FLERR,"Must have atom IDs for ndx2group command");
if (comm->me == 0) {
fp = fopen(arg[0], "r");
if (fp == nullptr)
error->one(FLERR,"Cannot open index file for reading");
if (screen)
fprintf(screen, "Reading groups from index file %s:\n",arg[0]);
if (logfile)
fprintf(logfile,"Reading groups from index file %s:\n",arg[0]);
error->one(FLERR,fmt::format("Cannot open index file for reading: {}",
utils::getsyserror()));
utils::logmesg(lmp,fmt::format("Reading groups from index file {}:\n",arg[0]));
}
if (narg == 1) { // restore all groups
if (narg == 1) { // restore all groups
do {
if (comm->me == 0) {
len = 0;
if (comm->me == 0) {
name = find_section(fp,"");
while (!name.empty()) {
// find the next section.
// if we had processed a section, before we need to step back
if (name != nullptr) {
rewind(fp);
char *tmp = find_section(fp,name);
delete[] tmp;
delete[] name;
name = nullptr;
// skip over group "all", which is called "System" in gromacs
if (name == "System") {
name = find_section(fp,"");
continue;
}
name = find_section(fp,nullptr);
if (name != nullptr) {
len=strlen(name)+1;
// skip over group "all", which is called "System" in gromacs
if (strcmp(name,"System") == 0) continue;
if (screen)
fprintf(screen," Processing group '%s'\n",name);
if (logfile)
fprintf(logfile," Processing group '%s'\n",name);
}
utils::logmesg(lmp,fmt::format(" Processing group '{}'\n",name));
len = name.size()+1;
MPI_Bcast(&len,1,MPI_INT,0,world);
if (len > 0) {
MPI_Bcast(name,len,MPI_CHAR,0,world);
if (len > 1) {
MPI_Bcast((void *)name.c_str(),len,MPI_CHAR,0,world);
// read tags for atoms in group and broadcast
num = 0;
tags = read_section(fp,num);
std::vector<tagint> tags = read_section(fp,next);
num = tags.size();
MPI_Bcast(&num,1,MPI_LMP_BIGINT,0,world);
MPI_Bcast(tags,num,MPI_LMP_TAGINT,0,world);
create(name,num,tags);
free(tags);
}
} else {
MPI_Bcast(&len,1,MPI_INT,0,world);
if (len > 0) {
delete[] name;
name = new char[len];
MPI_Bcast(name,len,MPI_CHAR,0,world);
MPI_Bcast(&num,1,MPI_LMP_BIGINT,0,world);
tags = (tagint *)malloc(sizeof(tagint)*(num ? num : 1));
MPI_Bcast(tags,num,MPI_LMP_TAGINT,0,world);
create(name,num,tags);
free(tags);
MPI_Bcast((void *)tags.data(),num,MPI_LMP_TAGINT,0,world);
create(name,tags);
name = next;
}
}
} while (len);
len = -1;
MPI_Bcast(&len,1,MPI_INT,0,world);
} else {
while (1) {
MPI_Bcast(&len,1,MPI_INT,0,world);
if (len < 0) break;
if (len > 1) {
char *buf = new char[len];
MPI_Bcast(buf,len,MPI_CHAR,0,world);
MPI_Bcast(&num,1,MPI_LMP_BIGINT,0,world);
tagint *tbuf = new tagint[num];
MPI_Bcast(tbuf,num,MPI_LMP_TAGINT,0,world);
create(buf,std::vector<tagint>(tbuf,tbuf+num));
delete[] buf;
delete[] tbuf;
}
}
}
} else { // restore selected groups
for (int idx=1; idx < narg; ++idx) {
for (int idx=1; idx < narg; ++idx) {
if (comm->me == 0) {
len = 0;
// find named section, search from beginning of file
if (name != nullptr) delete[] name;
rewind(fp);
name = find_section(fp,arg[idx]);
if (name != nullptr) len=strlen(name)+1;
if (screen)
fprintf(screen," %s group '%s'\n",
len ? "Processing" : "Skipping",arg[idx]);
if (logfile)
fprintf(logfile,"%s group '%s'\n",
len ? "Processing" : "Skipping",arg[idx]);
utils::logmesg(lmp,fmt::format(" {} group '{}'\n", name.size()
? "Processing" : "Skipping",arg[idx]));
len = name.size()+1;
MPI_Bcast(&len,1,MPI_INT,0,world);
if (len > 0) {
MPI_Bcast(name,len,MPI_CHAR,0,world);
if (len > 1) {
MPI_Bcast((void *)name.c_str(),len,MPI_CHAR,0,world);
// read tags for atoms in group and broadcast
num = 0;
tags = read_section(fp,num);
std::vector<tagint> tags = read_section(fp,name);
num = tags.size();
MPI_Bcast(&num,1,MPI_LMP_BIGINT,0,world);
MPI_Bcast(tags,num,MPI_LMP_TAGINT,0,world);
create(name,num,tags);
free(tags);
MPI_Bcast((void *)tags.data(),num,MPI_LMP_TAGINT,0,world);
create(name,tags);
}
} else {
MPI_Bcast(&len,1,MPI_INT,0,world);
if (len > 0) {
delete[] name;
name = new char[len];
MPI_Bcast(name,len,MPI_CHAR,0,world);
if (len > 1) {
char *buf = new char[len];
MPI_Bcast(buf,len,MPI_CHAR,0,world);
MPI_Bcast(&num,1,MPI_LMP_BIGINT,0,world);
tags = (tagint *)malloc(sizeof(tagint)*(num ? num : 1));
MPI_Bcast(tags,num,MPI_LMP_TAGINT,0,world);
create(name,num,tags);
free(tags);
tagint *tbuf = new tagint[num];
MPI_Bcast(tbuf,num,MPI_LMP_TAGINT,0,world);
create(buf,std::vector<tagint>(tbuf,tbuf+num));
delete[] buf;
delete[] tbuf;
}
}
}
}
delete[] name;
if (comm->me == 0) {
if (screen) fputs("\n",screen);
if (logfile) fputs("\n",logfile);
fclose(fp);
}
if (comm->me == 0) fclose(fp);
}
/* ---------------------------------------------------------------------- */
void Ndx2Group::create(char *name, bigint num, tagint *tags)
void Ndx2Group::create(const std::string &name, const std::vector<tagint> &tags)
{
// wipe out all members if the group exists. gid==0 is group "all"
int gid = group->find(name);
if (gid > 0) group->assign(std::string(name) + " clear");
if (gid > 0) group->assign(name + " clear");
// map from global to local
const int nlocal = atom->nlocal;
int *flags = (int *)calloc(nlocal,sizeof(int));
for (bigint i=0; i < num; ++i) {
for (bigint i=0; i < tags.size(); ++i) {
const int id = atom->map(tags[i]);
if (id < nlocal && id >= 0)
flags[id] = 1;
if (id < nlocal && id >= 0) flags[id] = 1;
}
group->create(name,flags);
free(flags);

3
src/USER-COLVARS/ndx_group.h
View File

@ -23,6 +23,7 @@ CommandStyle(ndx2group,Ndx2Group)
#define LMP_NDX_GROUP_H
#include "pointers.h"
#include <vector>
namespace LAMMPS_NS {
@ -30,7 +31,7 @@ class Ndx2Group : protected Pointers {
public:
Ndx2Group(class LAMMPS *lmp) : Pointers(lmp) {};
void command(int, char **);
void create(char *, bigint, tagint *);
void create(const std::string &, const std::vector<tagint> &);
};
}

1
src/USER-DIFFRACTION/compute_saed.cpp
View File

@ -30,6 +30,7 @@
#include <cmath>
#include <cstring>
#include <strings.h> // for strcasecmp()
#include "omp_compat.h"
using namespace LAMMPS_NS;

1
src/USER-DIFFRACTION/compute_xrd.cpp
View File

@ -31,6 +31,7 @@
#include <cmath>
#include <cstring>
#include <strings.h> // for strcasecmp()
#include "omp_compat.h"
using namespace LAMMPS_NS;

7
src/USER-DPD/fix_eos_table_rx.cpp
View File

@ -17,16 +17,15 @@
#include "fix_eos_table_rx.h"
#include <cstring>
#include <cmath>
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "comm.h"
#include "modify.h"
#include <cmath>
#include <cstring>
#define MAXLINE 1024

83
src/USER-DPD/fix_rx.cpp
View File

@ -13,28 +13,26 @@
#include "fix_rx.h"
#include <cstring>
#include <cmath>
#include <cfloat> // DBL_EPSILON
#include "atom.h"
#include "error.h"
#include "group.h"
#include "modify.h"
#include "force.h"
#include "memory.h"
#include "comm.h"
#include "update.h"
#include "domain.h"
#include "neighbor.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "math_special.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "math_special.h"
#include "neighbor.h"
#include "pair_dpd_fdt_energy.h"
#include "update.h"
#include <vector> // std::vector<>
#include <algorithm> // std::max
#include <cfloat> // DBL_EPSILON
#include <cmath>
#include <cstring>
#include <vector> // std::vector<>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -257,7 +255,6 @@ void FixRX::post_constructor()
bool match;
char **tmpspecies = new char*[maxspecies];
int tmpmaxstrlen = 0;
for (int jj=0; jj < maxspecies; jj++)
tmpspecies[jj] = nullptr;
@ -316,9 +313,7 @@ void FixRX::post_constructor()
if (!match) {
if (nUniqueSpecies+1>=maxspecies)
error->all(FLERR,"Exceeded the maximum number of species permitted in fix rx.");
tmpspecies[nUniqueSpecies] = new char[strlen(word)+1];
strcpy(tmpspecies[nUniqueSpecies],word);
tmpmaxstrlen = MAX(tmpmaxstrlen,(int)strlen(word));
tmpspecies[nUniqueSpecies] = utils::strdup(word);
nUniqueSpecies++;
}
word = strtok(nullptr, " \t\n\r\f");
@ -335,61 +330,31 @@ void FixRX::post_constructor()
id_fix_species = nullptr;
id_fix_species_old = nullptr;
n = strlen(id) + strlen("_SPECIES") + 1;
id_fix_species = new char[n];
n = strlen(id) + strlen("_SPECIES_OLD") + 1;
id_fix_species_old = new char[n];
id_fix_species = utils::strdup(std::string(id)+"_SPECIES");
id_fix_species_old = utils::strdup(std::string(id)+"_SPECIES_OLD");
strcpy(id_fix_species,id);
strcat(id_fix_species,"_SPECIES");
strcpy(id_fix_species_old,id);
strcat(id_fix_species_old,"_SPECIES_OLD");
char **newarg = new char*[nspecies+5];
char **newarg2 = new char*[nspecies+5];
newarg[0] = id_fix_species;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "property/atom";
newarg2[0] = id_fix_species_old;
newarg2[1] = group->names[igroup];
newarg2[2] = (char *) "property/atom";
char *str1 = new char[tmpmaxstrlen+3];
char *str2 = new char[tmpmaxstrlen+6];
const std::string fmtstr = "{} {} property/atom ";
auto newcmd1 = fmt::format(fmtstr,id_fix_species,group->names[igroup]);
auto newcmd2 = fmt::format(fmtstr,id_fix_species_old,group->names[igroup]);
for (int ii=0; ii<nspecies; ii++) {
strcpy(str1,"d_");
strcpy(str2,"d_");
strcat(str1,tmpspecies[ii]);
strcat(str2,tmpspecies[ii]);
strcat(str2,"Old");
newarg[ii+3] = new char[strlen(str1)+1];
newarg2[ii+3] = new char[strlen(str2)+1];
strcpy(newarg[ii+3],str1);
strcpy(newarg2[ii+3],str2);
newcmd1 += fmt::format(" d_{}",tmpspecies[ii]);
newcmd2 += fmt::format(" d_{}Old",tmpspecies[ii]);
}
delete[] str1;
delete[] str2;
newarg[nspecies+3] = (char *) "ghost";
newarg[nspecies+4] = (char *) "yes";
newarg2[nspecies+3] = (char *) "ghost";
newarg2[nspecies+4] = (char *) "yes";
newcmd1 += " ghost yes";
newcmd2 += " ghost yes";
modify->add_fix(nspecies+5,newarg,1);
modify->add_fix(newcmd1);
fix_species = (FixPropertyAtom *) modify->fix[modify->nfix-1];
restartFlag = modify->fix[modify->nfix-1]->restart_reset;
modify->add_fix(nspecies+5,newarg2,1);
modify->add_fix(newcmd2);
fix_species_old = (FixPropertyAtom *) modify->fix[modify->nfix-1];
if (nspecies==0) error->all(FLERR,"There are no rx species specified.");
for (int jj=0;jj<nspecies;jj++) {
delete[] tmpspecies[jj];
delete[] newarg[jj+3];
delete[] newarg2[jj+3];
}
delete[] newarg;
delete[] newarg2;
delete[] tmpspecies;
read_file( kineticsFile );

46
src/USER-DPD/pair_exp6_rx.cpp
View File

@ -13,20 +13,19 @@
#include "pair_exp6_rx.h"
#include <cmath>
#include <cstring>
#include <cfloat>
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "neigh_list.h"
#include "math_special.h"
#include "memory.h"
#include "error.h"
#include "modify.h"
#include "fix.h"
#include "neigh_list.h"
#include <cfloat>
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathSpecial;
@ -596,9 +595,7 @@ void PairExp6rx::coeff(int narg, char **arg)
if (nspecies==0) error->all(FLERR,"There are no rx species specified.");
read_file(arg[2]);
n = strlen(arg[3]) + 1;
site1 = new char[n];
strcpy(site1,arg[3]);
site1 = utils::strdup(arg[3]);
int ispecies;
for (ispecies = 0; ispecies < nspecies; ispecies++) {
@ -607,9 +604,7 @@ void PairExp6rx::coeff(int narg, char **arg)
if (ispecies == nspecies && strcmp(site1,"1fluid") != 0)
error->all(FLERR,"Site1 name not recognized in pair coefficients");
n = strlen(arg[4]) + 1;
site2 = new char[n];
strcpy(site2,arg[4]);
site2 = utils::strdup(arg[4]);
for (ispecies = 0; ispecies < nspecies; ispecies++) {
if (strcmp(site2,&atom->dname[ispecies][0]) == 0) break;
@ -806,17 +801,12 @@ void PairExp6rx::read_file(char *file)
params[nparams].ispecies = ispecies;
n = strlen(&atom->dname[ispecies][0]) + 1;
params[nparams].name = new char[n];
strcpy(params[nparams].name,&atom->dname[ispecies][0]);
n = strlen(words[1]) + 1;
params[nparams].potential = new char[n];
strcpy(params[nparams].potential,words[1]);
params[nparams].name = utils::strdup(&atom->dname[ispecies][0]);
params[nparams].potential = utils::strdup(words[1]);
if (strcmp(params[nparams].potential,"exp6") == 0) {
params[nparams].alpha = atof(words[2]);
params[nparams].epsilon = atof(words[3]);
params[nparams].rm = atof(words[4]);
params[nparams].alpha = utils::numeric(FLERR,words[2],false,lmp);
params[nparams].epsilon = utils::numeric(FLERR,words[3],false,lmp);
params[nparams].rm = utils::numeric(FLERR,words[4],false,lmp);
if (params[nparams].epsilon <= 0.0 || params[nparams].rm <= 0.0 ||
params[nparams].alpha < 0.0)
error->all(FLERR,"Illegal exp6/rx parameters. Rm and Epsilon must be greater than zero. Alpha cannot be negative.");
@ -842,11 +832,9 @@ void PairExp6rx::read_file2(char *file)
fp = nullptr;
if (comm->me == 0) {
fp = fopen(file,"r");
if (fp == nullptr) {
char str[128];
snprintf(str,128,"Cannot open polynomial file %s",file);
error->one(FLERR,str);
}
if (fp == nullptr)
error->one(FLERR,fmt::format("Cannot open polynomial file {}: {}",
file,utils::getsyserror()));
}
// one set of params can span multiple lines

22
src/USER-DPD/pair_multi_lucy.cpp
View File

@ -21,21 +21,19 @@
The Journal of Chemical Physics, 2016, 144, 104501.
------------------------------------------------------------------------------------------- */
#include "pair_multi_lucy.h"
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "neigh_list.h"
#include <cmath>
#include "math_const.h"
#include <cstring>
#include "pair_multi_lucy.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "citeme.h"
using namespace LAMMPS_NS;

35
src/USER-DPD/pair_multi_lucy_rx.cpp
View File

@ -21,23 +21,21 @@
The Journal of Chemical Physics, 2016, 144, 104501.
------------------------------------------------------------------------------------------- */
#include "pair_multi_lucy_rx.h"
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include <cmath>
#include "math_const.h"
#include <cstring>
#include "pair_multi_lucy_rx.h"
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "citeme.h"
#include "modify.h"
#include "fix.h"
using namespace LAMMPS_NS;
@ -387,14 +385,9 @@ void PairMultiLucyRX::coeff(int narg, char **arg)
bcast_table(tb);
nspecies = atom->nspecies_dpd;
int n;
n = strlen(arg[4]) + 1;
site1 = new char[n];
strcpy(site1,arg[4]);
n = strlen(arg[5]) + 1;
site2 = new char[n];
strcpy(site2,arg[5]);
site1 = utils::strdup(arg[4]);
site2 = utils::strdup(arg[5]);
// set table cutoff

10
src/USER-DPD/pair_table_rx.cpp
View File

@ -322,10 +322,8 @@ void PairTableRX::coeff(int narg, char **arg)
nspecies = atom->nspecies_dpd;
if (nspecies==0) error->all(FLERR,"There are no rx species specified.");
int n;
n = strlen(arg[4]) + 1;
site1 = new char[n];
strcpy(site1,arg[4]);
site1 = utils::strdup(arg[4]);
int ispecies;
for (ispecies = 0; ispecies < nspecies; ispecies++) {
@ -334,9 +332,7 @@ void PairTableRX::coeff(int narg, char **arg)
if (ispecies == nspecies && strcmp(site1,"1fluid") != 0)
error->all(FLERR,"Site1 name not recognized in pair coefficients");
n = strlen(arg[5]) + 1;
site2 = new char[n];
strcpy(site2,arg[5]);
site2 = utils::strdup(arg[5]);
for (ispecies = 0; ispecies < nspecies; ispecies++) {
if (strcmp(site2,&atom->dname[ispecies][0]) == 0) break;

17
src/USER-DRUDE/compute_temp_drude.cpp
View File

@ -13,15 +13,16 @@
#include "compute_temp_drude.h"
#include <cstring>
#include "atom.h"
#include "update.h"
#include "force.h"
#include "modify.h"
#include "fix_drude.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "comm.h"
#include "fix_drude.h"
#include "force.h"
#include "modify.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
@ -116,9 +117,7 @@ int ComputeTempDrude::modify_param(int narg, char **arg)
if (strcmp(arg[0],"temp") == 0) {
if (narg < 2) error->all(FLERR,"Illegal fix_modify command");
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
id_temp = utils::strdup(arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0)

32
src/USER-DRUDE/fix_tgnh_drude.cpp
View File

@ -16,23 +16,25 @@
---------------------------------------------------------------------------------------- */
#include "fix_tgnh_drude.h"
#include <cstring>
#include <cmath>
#include "atom.h"
#include "comm.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "fix_deform.h"
#include "force.h"
#include "group.h"
#include "comm.h"
#include "neighbor.h"
#include "irregular.h"
#include "modify.h"
#include "fix_deform.h"
#include "compute.h"
#include "kspace.h"
#include "update.h"
#include "respa.h"
#include "domain.h"
#include "memory.h"
#include "error.h"
#include "modify.h"
#include "neighbor.h"
#include "respa.h"
#include "update.h"
#include <cstring>
#include <cmath>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -1428,9 +1430,7 @@ int FixTGNHDrude::modify_param(int narg, char **arg)
tcomputeflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
id_temp = utils::strdup(arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0)
@ -1462,9 +1462,7 @@ int FixTGNHDrude::modify_param(int narg, char **arg)
pcomputeflag = 0;
}
delete [] id_press;
int n = strlen(arg[1]) + 1;
id_press = new char[n];
strcpy(id_press,arg[1]);
id_press = utils::strdup(arg[1]);
int icompute = modify->find_compute(arg[1]);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify pressure ID");

19
src/USER-DRUDE/fix_tgnpt_drude.cpp
View File

@ -12,10 +12,9 @@
------------------------------------------------------------------------- */
#include "fix_tgnpt_drude.h"
#include <cstring>
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -35,23 +34,15 @@ FixTGNPTDrude::FixTGNPTDrude(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp, tcmd.c_str());
tcmd += " all temp";
modify->add_compute(tcmd);
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
std::string pcmd = id + std::string("_press");
id_press = new char[pcmd.size()+1];
strcpy(id_press, pcmd.c_str());
pcmd += " all pressure " + std::string(id_temp);
modify->add_compute(pcmd);
id_press = utils::strdup(std::string(id) + "_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
pcomputeflag = 1;
}

13
src/USER-DRUDE/fix_tgnvt_drude.cpp
View File

@ -12,11 +12,12 @@
------------------------------------------------------------------------- */
#include "fix_tgnvt_drude.h"
#include <cstring>
#include "error.h"
#include "group.h"
#include "modify.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -34,11 +35,7 @@ FixTGNVTDrude::FixTGNVTDrude(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp style
// id = fix-ID + temp
std::string tcmd = id + std::string("_temp");
id_temp = new char[tcmd.size()+1];
strcpy(id_temp, tcmd.c_str());
tcmd += fmt::format(" {} temp", group->names[igroup]);
modify->add_compute(tcmd);
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp",id_temp,group->names[igroup]));
tcomputeflag = 1;
}

17
src/USER-EFF/compute_temp_region_eff.cpp
View File

@ -16,17 +16,18 @@
------------------------------------------------------------------------- */
#include <cstring>
#include "compute_temp_region_eff.h"
#include "atom.h"
#include "update.h"
#include "force.h"
#include "domain.h"
#include "region.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "memory.h"
#include "error.h"
#include "region.h"
#include "update.h"
#include <cstring>
using namespace LAMMPS_NS;
@ -43,9 +44,7 @@ ComputeTempRegionEff::ComputeTempRegionEff(LAMMPS *lmp, int narg, char **arg) :
iregion = domain->find_region(arg[3]);
if (iregion == -1)
error->all(FLERR,"Region ID for compute temp/region/eff does not exist");
int n = strlen(arg[3]) + 1;
idregion = new char[n];
strcpy(idregion,arg[3]);
idregion = utils::strdup(arg[3]);
scalar_flag = vector_flag = 1;
size_vector = 6;

32
src/USER-EFF/fix_nph_eff.cpp
View File

@ -11,10 +11,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <cstring>
#include "fix_nph_eff.h"
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -34,35 +34,15 @@ FixNPHEff::FixNPHEff(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp/eff";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} all temp/eff",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
id_press = utils::strdup(std::string(id) + "_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
pcomputeflag = 1;
}

32
src/USER-EFF/fix_npt_eff.cpp
View File

@ -11,10 +11,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <cstring>
#include "fix_npt_eff.h"
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -34,35 +34,15 @@ FixNPTEff::FixNPTEff(LAMMPS *lmp, int narg, char **arg) :
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp/eff";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} all temp/eff",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
id_press = utils::strdup(std::string(id) + "_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
pcomputeflag = 1;
}

21
src/USER-EFF/fix_nvt_eff.cpp
View File

@ -11,11 +11,13 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <cstring>
#include "fix_nvt_eff.h"
#include "error.h"
#include "group.h"
#include "modify.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -33,17 +35,8 @@ FixNVTEff::FixNVTEff(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp style
// id = fix-ID + temp
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "temp/eff";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp/eff",
id_temp,group->names[igroup]));
tcomputeflag = 1;
}

32
src/USER-EFF/fix_nvt_sllod_eff.cpp
View File

@ -11,19 +11,20 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <cmath>
#include <cstring>
#include "fix_nvt_sllod_eff.h"
#include "math_extra.h"
#include "atom.h"
#include "compute.h"
#include "domain.h"
#include "group.h"
#include "modify.h"
#include "error.h"
#include "fix.h"
#include "fix_deform.h"
#include "compute.h"
#include "error.h"
#include "group.h"
#include "math_extra.h"
#include "modify.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -45,18 +46,9 @@ FixNVTSllodEff::FixNVTSllodEff(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp style
// id = fix-ID + temp
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "temp/deform/eff";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} tmp/deform/eff",
id_temp,group->names[igroup]));
tcomputeflag = 1;
}

33
src/USER-EFF/fix_temp_rescale_eff.cpp
View File

@ -15,18 +15,19 @@
Contributing author: Andres Jaramillo-Botero
------------------------------------------------------------------------- */
#include <cstring>
#include <cmath>
#include "fix_temp_rescale_eff.h"
#include "atom.h"
#include "force.h"
#include "group.h"
#include "update.h"
#include "comm.h"
#include "modify.h"
#include "compute.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "modify.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -56,17 +57,9 @@ FixTempRescaleEff::FixTempRescaleEff(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp/eff
// id = fix-ID + temp, compute group = fix group
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[6];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "temp/eff";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp/eff",
id_temp,group->names[igroup]));
tflag = 1;
energy = 0.0;
@ -172,9 +165,7 @@ int FixTempRescaleEff::modify_param(int narg, char **arg)
tflag = 0;
}
delete [] id_temp;
int n = strlen(arg[1]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[1]);
id_temp = utils::strdup(arg[1]);
int icompute = modify->find_compute(id_temp);
if (icompute < 0) error->all(FLERR,"Could not find fix_modify temperature ID");

36
src/USER-INTEL/fix_npt_intel.cpp
View File

@ -11,10 +11,10 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <cstring>
#include "fix_npt_intel.h"
#include "modify.h"
#include "error.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace FixConst;
@ -25,44 +25,24 @@ FixNPTIntel::FixNPTIntel(LAMMPS *lmp, int narg, char **arg) :
FixNHIntel(lmp, narg, arg)
{
if (!tstat_flag)
error->all(FLERR,"Temperature control must be used with fix npt/omp");
error->all(FLERR,"Temperature control must be used with fix npt/intel");
if (!pstat_flag)
error->all(FLERR,"Pressure control must be used with fix npt/omp");
error->all(FLERR,"Pressure control must be used with fix npt/intl");
// create a new compute temp style
// id = fix-ID + temp
// compute group = all since pressure is always global (group all)
// and thus its KE/temperature contribution should use group all
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = (char *) "all";
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} all temp",id_temp));
tcomputeflag = 1;
// create a new compute pressure style
// id = fix-ID + press, compute group = all
// pass id_temp as 4th arg to pressure constructor
n = strlen(id) + 7;
id_press = new char[n];
strcpy(id_press,id);
strcat(id_press,"_press");
newarg = new char*[4];
newarg[0] = id_press;
newarg[1] = (char *) "all";
newarg[2] = (char *) "pressure";
newarg[3] = id_temp;
modify->add_compute(4,newarg);
delete [] newarg;
id_press = utils::strdup(std::string(id) + "_press");
modify->add_compute(fmt::format("{} all pressure {}",id_press, id_temp));
pcomputeflag = 1;
}

24
src/USER-INTEL/fix_nvt_intel.cpp
View File

@ -11,11 +11,13 @@
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include <cstring>
#include "fix_nvt_intel.h"
#include "error.h"
#include "group.h"
#include "modify.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -26,25 +28,15 @@ FixNVTIntel::FixNVTIntel(LAMMPS *lmp, int narg, char **arg) :
FixNHIntel(lmp, narg, arg)
{
if (!tstat_flag)
error->all(FLERR,"Temperature control must be used with fix nvt");
error->all(FLERR,"Temperature control must be used with fix nvt/intel");
if (pstat_flag)
error->all(FLERR,"Pressure control can not be used with fix nvt");
error->all(FLERR,"Pressure control can not be used with fix nvt/intel");
// create a new compute temp style
// id = fix-ID + temp
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp",id_temp,group->names[igroup]));
tcomputeflag = 1;
}

15
src/USER-INTEL/fix_nvt_sllod_intel.cpp
View File

@ -46,18 +46,9 @@ FixNVTSllodIntel::FixNVTSllodIntel(LAMMPS *lmp, int narg, char **arg) :
// create a new compute temp style
// id = fix-ID + temp
int n = strlen(id) + 6;
id_temp = new char[n];
strcpy(id_temp,id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "temp/deform";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp/deform",
id_temp,group->names[igroup]));
tcomputeflag = 1;
}

24
src/USER-LB/fix_lb_fluid.cpp
View File

@ -17,22 +17,22 @@
------------------------------------------------------------------------- */
#include "fix_lb_fluid.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "group.h"
#include "memory.h"
#include "modify.h"
#include "random_mars.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include <algorithm>
#include <utility>
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "domain.h"
#include "atom.h"
#include "group.h"
#include "random_mars.h"
#include "update.h"
#include "force.h"
#include "modify.h"
using namespace LAMMPS_NS;
using namespace FixConst;

18
src/USER-MANIFOLD/fix_manifoldforce.cpp
View File

@ -58,22 +58,16 @@ FixManifoldForce::FixManifoldForce(LAMMPS *lmp, int narg, char **arg) :
ptr_m = create_manifold(m_name,lmp,narg,arg);
// Construct manifold from factory:
if (!ptr_m) {
char msg[2048];
snprintf(msg,2048,"Manifold pointer for manifold '%s' was NULL for some reason", arg[3]);
error->all(FLERR,msg);
}
if (!ptr_m)
error->all(FLERR,fmt::format("Manifold pointer for manifold '{}' "
"was NULL for some reason", arg[3]));
// After constructing the manifold, you can safely make
// room for the parameters
nvars = ptr_m->nparams();
if (narg < nvars+4) {
char msg[2048];
sprintf(msg,"Manifold %s needs at least %d argument(s)!",
m_name, nvars);
error->all(FLERR,msg);
}
if (narg < nvars+4)
error->all(FLERR,fmt::format("Manifold {} needs at least {} "
"argument(s)!", m_name, nvars));
ptr_m->params = new double[nvars];
if (ptr_m->params == nullptr) {

27
src/USER-MANIFOLD/fix_nve_manifold_rattle.cpp
View File

@ -31,10 +31,8 @@
------------------------------------------------------------------------- */
#include "fix_nve_manifold_rattle.h"
#include <cstring>
#include "atom.h"
#include "force.h"
#include "update.h"
@ -44,18 +42,17 @@
#include "citeme.h"
#include "comm.h"
#include <cstring>
#include "manifold_factory.h"
#include "manifold.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace user_manifold;
enum { CONST, EQUAL }; // For treating the variables.
static const char* cite_fix_nve_manifold_rattle =
"fix nve/manifold/rattle command:\n\n"
"@article{paquay-2016,\n"
@ -95,9 +92,7 @@ FixNVEManifoldRattle::FixNVEManifoldRattle( LAMMPS *lmp, int &narg, char **arg,
max_iter = utils::numeric( FLERR, arg[4] ,false,lmp);
ptr_m = create_manifold(arg[5], lmp, narg, arg);
if (!ptr_m) {
error->all(FLERR,"Error creating manifold pointer");
}
if (!ptr_m) error->all(FLERR,"Error creating manifold pointer");
nvars = ptr_m->nparams();
tstrs = new char*[nvars];
@ -105,17 +100,13 @@ FixNVEManifoldRattle::FixNVEManifoldRattle( LAMMPS *lmp, int &narg, char **arg,
tstyle = new int[nvars];
is_var = new int[nvars];
if (!tstrs || !tvars || !tstyle || !is_var) {
if (!tstrs || !tvars || !tstyle || !is_var)
error->all(FLERR, "Error creating manifold arg arrays");
}
// Check if you have enough args:
if (6 + nvars > narg) {
char msg[2048];
sprintf(msg, "Not enough args for manifold %s, %d expected but got %d\n",
ptr_m->id(), nvars, narg - 6);
error->all(FLERR, msg);
}
if (6 + nvars > narg)
error->all(FLERR,fmt::format("Not enough args for manifold {}, {} expected "
"but got {}\n",ptr_m->id(),nvars, narg - 6));
// Loop over manifold args:
for (int i = 0; i < nvars; ++i) {
int len = 0, offset = 0;
@ -154,9 +145,7 @@ FixNVEManifoldRattle::FixNVEManifoldRattle( LAMMPS *lmp, int &narg, char **arg,
}
argi += 2;
} else if (error_on_unknown_keyword) {
char msg[2048];
sprintf(msg,"Error parsing arg \"%s\".\n", arg[argi]);
error->all(FLERR, msg);
error->all(FLERR,fmt::format("Error parsing arg \"{}\".\n",arg[argi]));
} else {
argi += 1;
}

35
src/USER-MANIFOLD/fix_nvt_manifold_rattle.cpp
View File

@ -131,24 +131,10 @@ FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
}
// Create temperature compute:
const char *fix_id = arg[1];
int n = strlen(fix_id)+6;
id_temp = new char[n];
strcpy(id_temp,fix_id);
strcat(id_temp,"_temp");
char **newarg = new char*[3];
newarg[0] = id_temp;
newarg[1] = group->names[igroup];
newarg[2] = (char*) "temp";
modify->add_compute(3,newarg);
delete [] newarg;
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp",id_temp,group->names[igroup]));
int icompute = modify->find_compute(id_temp);
if (icompute < 0) {
error->all(FLERR,"Temperature ID for fix nvt/manifold/rattle "
"does not exist");
}
temperature = modify->compute[icompute];
if (temperature->tempbias) which = BIAS;
else which = NOBIAS;
@ -167,8 +153,6 @@ FixNVTManifoldRattle::FixNVTManifoldRattle(LAMMPS *lmp, int narg, char **arg,
for (int ich = 0; ich < mtchain; ++ich) {
eta[ich] = eta_dot[ich] = eta_dotdot[ich] = 0.0;
}
}
/* ---------------------------------------------------------------------- */
@ -185,9 +169,6 @@ FixNVTManifoldRattle::~FixNVTManifoldRattle()
if (id_temp) delete[] id_temp;
}
int FixNVTManifoldRattle::setmask()
{
int mask = 0;
@ -198,7 +179,6 @@ int FixNVTManifoldRattle::setmask()
return mask;
}
/* --------------------------------------------------------------------------
Check that force modification happens before position and velocity update.
Make sure respa is not used.
@ -219,8 +199,6 @@ void FixNVTManifoldRattle::init()
}
void FixNVTManifoldRattle::setup(int /*vflag*/)
{
compute_temp_target();
@ -359,9 +337,6 @@ void FixNVTManifoldRattle::nh_v_temp()
}
}
// Most of this logic is based on fix_nh:
void FixNVTManifoldRattle::initial_integrate(int /*vflag*/)
{
@ -381,8 +356,6 @@ void FixNVTManifoldRattle::final_integrate()
nhc_temp_integrate();
}
/* ---------------------------------------------------------------------- */
void FixNVTManifoldRattle::reset_dt()
{
@ -395,10 +368,6 @@ void FixNVTManifoldRattle::reset_dt()
}
double FixNVTManifoldRattle::memory_usage()
{
double bytes = FixNVEManifoldRattle::memory_usage();

16
src/USER-MGPT/pair_mgpt.cpp
View File

@ -23,20 +23,18 @@
#include "pair_mgpt.h"
#include <cmath>
#include <cstring>
#include <cassert>
#include "atom.h"
#include "force.h"
#include "comm.h"
#include "neighbor.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "memory.h"
#include "error.h"
#include "neighbor.h"
#include <cassert>
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;

13
src/USER-MISC/compute_gyration_shape.cpp
View File

@ -15,17 +15,18 @@
* Contributing author: Evangelos Voyiatzis (Royal DSM)
* ------------------------------------------------------------------------- */
#include "compute_gyration_shape.h"
#include <cmath>
#include <cstring>
#include "error.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "math_extra.h"
#include "math_special.h"
#include "modify.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
@ -41,9 +42,7 @@ ComputeGyrationShape::ComputeGyrationShape(LAMMPS *lmp, int narg, char **arg) :
extvector = 0;
// ID of compute gyration
int n = strlen(arg[3]) + 1;
id_gyration = new char[n];
strcpy(id_gyration,arg[3]);
id_gyration = utils::strdup(arg[3]);
init();

21
src/USER-MISC/compute_gyration_shape_chunk.cpp
View File

@ -15,17 +15,18 @@
* Contributing author: Evangelos Voyiatzis (Royal DSM)
* ------------------------------------------------------------------------- */
#include "compute_gyration_shape_chunk.h"
#include "error.h"
#include "math_eigen.h"
#include "math_extra.h"
#include "math_special.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include "error.h"
#include "math_extra.h"
#include "math_eigen.h"
#include "math_special.h"
#include "modify.h"
#include "memory.h"
#include "update.h"
using namespace LAMMPS_NS;
@ -37,9 +38,7 @@ ComputeGyrationShapeChunk::ComputeGyrationShapeChunk(LAMMPS *lmp, int narg, char
if (narg != 4) error->all(FLERR,"Illegal compute gyration/shape/chunk command");
// ID of compute gyration
int n = strlen(arg[3]) + 1;
id_gyration_chunk = new char[n];
strcpy(id_gyration_chunk,arg[3]);
id_gyration_chunk = utils::strdup(arg[3]);
init();

24
src/USER-MISC/compute_hma.cpp
View File

@ -78,11 +78,7 @@ ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) :
if (narg < 4) error->all(FLERR,"Illegal compute hma command");
if (igroup) error->all(FLERR,"Compute hma must use group all");
if (strcmp(arg[3],"NULL") == 0) {error->all(FLERR,"fix ID specifying the set temperature of canonical simulation is required");}
else {
int n = strlen(arg[3]) + 1;
id_temp = new char[n];
strcpy(id_temp,arg[3]);
}
else id_temp = utils::strdup(arg[3]);
create_attribute = 1;
extscalar = 1;
@ -92,23 +88,11 @@ ComputeHMA::ComputeHMA(LAMMPS *lmp, int narg, char **arg) :
// our new fix's id (id_fix)= compute-ID + COMPUTE_STORE
// our new fix's group = same as compute group
int n = strlen(id) + strlen("_COMPUTE_STORE") + 1;
id_fix = new char[n];
strcpy(id_fix,id);
strcat(id_fix,"_COMPUTE_STORE");
char **newarg = new char*[6];
newarg[0] = id_fix;
newarg[1] = group->names[igroup];
newarg[2] = (char *) "STORE";
newarg[3] = (char *) "peratom";
newarg[4] = (char *) "1";
newarg[5] = (char *) "3";
modify->add_fix(6,newarg);
id_fix = utils::strdup(std::string(id)+"_COMPUTE_STORE");
modify->add_fix(fmt::format("{} {} STORE peratom 1 3",
id_fix, group->names[igroup]));
fix = (FixStore *) modify->fix[modify->nfix-1];
delete [] newarg;
// calculate xu,yu,zu for fix store array
// skip if reset from restart file

8
src/USER-MISC/compute_pressure_grem.cpp
View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
#include "compute_pressure_grem.h"
#include <cstring>
#include "update.h"
#include "domain.h"
#include "modify.h"
@ -21,6 +21,8 @@
#include "kspace.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ----------------------------------------------------------------------
@ -30,9 +32,7 @@ using namespace LAMMPS_NS;
ComputePressureGrem::ComputePressureGrem(LAMMPS *lmp, int narg, char **arg) :
ComputePressure(lmp, narg-1, arg)
{
int len = strlen(arg[narg-1])+1;
fix_grem = new char[len];
strcpy(fix_grem,arg[narg-1]);
fix_grem = utils::strdup(arg[narg-1]);
}
/* ---------------------------------------------------------------------- */

73
src/USER-MISC/dihedral_table.cpp
View File

@ -464,7 +464,6 @@ DihedralTable::~DihedralTable()
if (allocated) {
memory->destroy(setflag);
//memory->destroy(phi0); <- equilibrium angles not supported
memory->destroy(tabindex);
}
}
@ -718,49 +717,18 @@ void DihedralTable::compute(int eflag, int vflag)
}
} // void DihedralTable::compute()
// single() calculates the dihedral-angle energy of atoms i1, i2, i3, i4.
double DihedralTable::single(int type, int i1, int i2, int i3, int i4)
{
double vb12[g_dim];
double vb23[g_dim];
double vb34[g_dim];
double n123[g_dim];
double n234[g_dim];
double **x = atom->x;
double phi = Phi(x[i1], x[i2], x[i3], x[i4], domain,
vb12, vb23, vb34, n123, n234);
double u=0.0;
u_lookup(type, phi, u); //Calculate the energy, and store it in "u"
return u;
}
/* ---------------------------------------------------------------------- */
void DihedralTable::allocate()
{
allocated = 1;
int n = atom->ndihedraltypes;
memory->create(tabindex,n+1,"dihedral:tabindex");
//memory->create(phi0,n+1,"dihedral:phi0"); <-equilibrium angles not supported
memory->create(setflag,n+1,"dihedral:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
@ -791,16 +759,13 @@ void DihedralTable::settings(int narg, char **arg)
tables = nullptr;
}
/* ----------------------------------------------------------------------
set coeffs for one type
------------------------------------------------------------------------- */
void DihedralTable::coeff(int narg, char **arg)
{
if (narg != 3) error->all(FLERR,"Illegal dihedral_coeff command");
if (narg != 3) error->all(FLERR,"Incorrect args for dihedral coefficients");
if (!allocated) allocate();
int ilo,ihi;
@ -815,25 +780,22 @@ void DihedralTable::coeff(int narg, char **arg)
if (me == 0) read_table(tb,arg[1],arg[2]);
bcast_table(tb);
// --- check the angle data for range errors ---
// --- and resolve issues with periodicity ---
if (tb->ninput < 2) {
error->one(FLERR,fmt::format("Invalid dihedral table length ({}).",
arg[2]));
} else if ((tb->ninput == 2) && (tabstyle == SPLINE)) {
error->one(FLERR,fmt::format("Invalid dihedral spline table length. "
"(Try linear)\n ({}).",arg[2]));
}
if (tb->ninput < 2)
error->all(FLERR,fmt::format("Invalid dihedral table length: {}",arg[2]));
else if ((tb->ninput == 2) && (tabstyle == SPLINE))
error->all(FLERR,fmt::format("Invalid dihedral spline table length: {} "
"(Try linear)",arg[2]));
// check for monotonicity
for (int i=0; i < tb->ninput-1; i++) {
if (tb->phifile[i] >= tb->phifile[i+1]) {
auto err_msg = fmt::format("Dihedral table values are not increasing "
"({}, {}th entry)",arg[2],i+1);
"({}, entry {})",arg[2],i+1);
if (i==0)
err_msg += std::string("\n(This is probably a mistake with your table format.)\n");
err_msg += "\n(This is probably a mistake with your table format.)\n";
error->all(FLERR,err_msg);
}
}
@ -971,15 +933,12 @@ void DihedralTable::coeff(int narg, char **arg)
}
ntables++;
if (count == 0)
error->all(FLERR,"Illegal dihedral_coeff command");
} //DihedralTable::coeff()
if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients");
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void DihedralTable::write_restart(FILE *fp)
{
@ -1065,7 +1024,6 @@ void DihedralTable::read_table(Table *tb, char *file, char *keyword)
error->one(FLERR,"Did not find keyword in table file");
}
// read args on 2nd line of section
// allocate table arrays for file values
@ -1183,10 +1141,8 @@ void DihedralTable::spline_table(Table *tb)
}
} // for (int i=0; i<tb->ninput; i++)
if ((num_disagreements > tb->ninput/2) && (num_disagreements > 2)) {
std::string msg("Dihedral table has inconsistent forces and energies. (Try \"NOF\".)\n");
error->all(FLERR, msg);
}
if ((num_disagreements > tb->ninput/2) && (num_disagreements > 2))
error->all(FLERR,"Dihedral table has inconsistent forces and energies. (Try \"NOF\".)\n");
} // check for consistency if (! tb->f_unspecified)
@ -1281,7 +1237,6 @@ void DihedralTable::compute_table(Table *tb)
void DihedralTable::param_extract(Table *tb, char *line)
{
//tb->theta0 = 180.0; <- equilibrium angles not supported
tb->ninput = 0;
tb->f_unspecified = false; //default
tb->use_degrees = true; //default

8
src/USER-MISC/dihedral_table.h
View File

@ -25,7 +25,6 @@ DihedralStyle(table,DihedralTable)
#define LMP_DIHEDRAL_TABLE_H
#include "dihedral.h"
namespace LAMMPS_NS {
class DihedralTable : public Dihedral {
@ -34,7 +33,7 @@ class DihedralTable : public Dihedral {
virtual ~DihedralTable();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
virtual void coeff(int, char **);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
@ -43,13 +42,11 @@ class DihedralTable : public Dihedral {
protected:
int tabstyle,tablength;
// double *phi0; <- equilibrium angles not supported
std::string checkU_fname;
std::string checkF_fname;
struct Table {
int ninput;
//double phi0; <-equilibrium angles not supported
int f_unspecified; // boolean (but MPI does not like type "bool")
int use_degrees; // boolean (but MPI does not like type "bool")
double *phifile,*efile,*ffile;
@ -62,7 +59,7 @@ class DihedralTable : public Dihedral {
Table *tables;
int *tabindex;
void allocate();
virtual void allocate();
void null_table(Table *);
void free_table(Table *);
void read_table(Table *, char *, char *);
@ -152,6 +149,5 @@ class DihedralTable : public Dihedral {
} // namespace LAMMPS_NS
#endif //#ifndef LMP_DIHEDRAL_TABLE_H
#endif //#ifdef DIHEDRAL_CLASS ... #else

842
src/USER-MISC/dihedral_table_cut.cpp
View File

@ -18,28 +18,23 @@
#include "dihedral_table_cut.h"
#include <cctype>
#include "atom.h"
#include "citeme.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "neighbor.h"
#include "update.h"
#include <cmath>
#include <cstring>
#include <fstream> // IWYU pragma: keep
#include <sstream> // IWYU pragma: keep
#include "atom.h"
#include "neighbor.h"
#include "update.h"
#include "comm.h"
#include "force.h"
#include "citeme.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
using namespace MathConst;
using namespace std;
static const char cite_dihedral_tablecut[] =
"dihedral_style table/cut command:\n\n"
@ -89,272 +84,6 @@ enum { //GSL status return codes.
};
static int solve_cyc_tridiag( const double diag[], size_t d_stride,
const double offdiag[], size_t o_stride,
const double b[], size_t b_stride,
double x[], size_t x_stride,
size_t N, bool warn)
{
int status = GSL_SUCCESS;
double * delta = (double *) malloc (N * sizeof (double));
double * gamma = (double *) malloc (N * sizeof (double));
double * alpha = (double *) malloc (N * sizeof (double));
double * c = (double *) malloc (N * sizeof (double));
double * z = (double *) malloc (N * sizeof (double));
if (delta == 0 || gamma == 0 || alpha == 0 || c == 0 || z == 0) {
if (warn)
fprintf(stderr,"Internal Cyclic Spline Error: failed to allocate working space\n");
if (delta) free(delta);
if (gamma) free(gamma);
if (alpha) free(alpha);
if (c) free(c);
if (z) free(z);
return GSL_ENOMEM;
}
else
{
size_t i, j;
double sum = 0.0;
/* factor */
if (N == 1)
{
x[0] = b[0] / diag[0];
free(delta);
free(gamma);
free(alpha);
free(c);
free(z);
return GSL_SUCCESS;
}
alpha[0] = diag[0];
gamma[0] = offdiag[0] / alpha[0];
delta[0] = offdiag[o_stride * (N-1)] / alpha[0];
if (alpha[0] == 0) {
status = GSL_EZERODIV;
}
for (i = 1; i < N - 2; i++)
{
alpha[i] = diag[d_stride * i] - offdiag[o_stride * (i-1)] * gamma[i - 1];
gamma[i] = offdiag[o_stride * i] / alpha[i];
delta[i] = -delta[i - 1] * offdiag[o_stride * (i-1)] / alpha[i];
if (alpha[i] == 0) {
status = GSL_EZERODIV;
}
}
for (i = 0; i < N - 2; i++)
{
sum += alpha[i] * delta[i] * delta[i];
}
alpha[N - 2] = diag[d_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * gamma[N - 3];
gamma[N - 2] = (offdiag[o_stride * (N - 2)] - offdiag[o_stride * (N - 3)] * delta[N - 3]) / alpha[N - 2];
alpha[N - 1] = diag[d_stride * (N - 1)] - sum - alpha[(N - 2)] * gamma[N - 2] * gamma[N - 2];
/* update */
z[0] = b[0];
for (i = 1; i < N - 1; i++)
{
z[i] = b[b_stride * i] - z[i - 1] * gamma[i - 1];
}
sum = 0.0;
for (i = 0; i < N - 2; i++)
{
sum += delta[i] * z[i];
}
z[N - 1] = b[b_stride * (N - 1)] - sum - gamma[N - 2] * z[N - 2];
for (i = 0; i < N; i++)
{
c[i] = z[i] / alpha[i];
}
/* backsubstitution */
x[x_stride * (N - 1)] = c[N - 1];
x[x_stride * (N - 2)] = c[N - 2] - gamma[N - 2] * x[x_stride * (N - 1)];
if (N >= 3)
{
for (i = N - 3, j = 0; j <= N - 3; j++, i--)
{
x[x_stride * i] = c[i] - gamma[i] * x[x_stride * (i + 1)] - delta[i] * x[x_stride * (N - 1)];
}
}
}
free (z);
free (c);
free (alpha);
free (gamma);
free (delta);
if ((status == GSL_EZERODIV) && warn)
fprintf(stderr, "Internal Cyclic Spline Error: Matrix must be positive definite.\n");
return status;
} //solve_cyc_tridiag()
/* ----------------------------------------------------------------------
spline and splint routines modified from Numerical Recipes
------------------------------------------------------------------------- */
static int cyc_spline(double const *xa,
double const *ya,
int n,
double period,
double *y2a, bool warn)
{
double *diag = new double[n];
double *offdiag = new double[n];
double *rhs = new double[n];
double xa_im1, xa_ip1;
// In the cyclic case, there are n equations with n unknows.
// The for loop sets up the equations we need to solve.
// Later we invoke the GSL tridiagonal matrix solver to solve them.
for (int i=0; i < n; i++) {
// I have to lookup xa[i+1] and xa[i-1]. This gets tricky because of
// periodic boundary conditions. We handle that now.
int im1 = i-1;
if (im1<0) {
im1 += n;
xa_im1 = xa[im1] - period;
}
else
xa_im1 = xa[im1];
int ip1 = i+1;
if (ip1>=n) {
ip1 -= n;
xa_ip1 = xa[ip1] + period;
}
else
xa_ip1 = xa[ip1];
// Recall that we want to find the y2a[] parameters (there are n of them).
// To solve for them, we have a linear equation with n unknowns
// (in the cyclic case that is). For details, the non-cyclic case is
// explained in equation 3.3.7 in Numerical Recipes in C, p. 115.
diag[i] = (xa_ip1 - xa_im1) / 3.0;
offdiag[i] = (xa_ip1 - xa[i]) / 6.0;
rhs[i] = ((ya[ip1] - ya[i]) / (xa_ip1 - xa[i])) -
((ya[i] - ya[im1]) / (xa[i] - xa_im1));
}
// Because this matrix is tridiagonal (and cyclic), we can use the following
// cheap method to invert it.
if (solve_cyc_tridiag(diag, 1,
offdiag, 1,
rhs, 1,
y2a, 1,
n, warn) != GSL_SUCCESS) {
if (warn)
fprintf(stderr,"Error in inverting matrix for splines.\n");
delete [] diag;
delete [] offdiag;
delete [] rhs;
return 1;
}
delete [] diag;
delete [] offdiag;
delete [] rhs;
return 0;
} // cyc_spline()
/* ---------------------------------------------------------------------- */
// cyc_splint(): Evaluates a spline at position x, with n control
// points located at xa[], ya[], and with parameters y2a[]
// The xa[] must be monotonically increasing and their
// range should not exceed period (ie xa[n-1] < xa[0] + period).
// x must lie in the range: [(xa[n-1]-period), (xa[0]+period)]
// "period" is typically 2*PI.
static double cyc_splint(double const *xa,
double const *ya,
double const *y2a,
int n,
double period,
double x)
{
int klo = -1;
int khi = n;
int k;
double xlo = xa[n-1] - period;
double xhi = xa[0] + period;
while (khi-klo > 1) {
k = (khi+klo) >> 1; //(k=(khi+klo)/2)
if (xa[k] > x) {
khi = k;
xhi = xa[k];
}
else {
klo = k;
xlo = xa[k];
}
}
if (khi == n) khi = 0;
if (klo ==-1) klo = n-1;
double h = xhi-xlo;
double a = (xhi-x) / h;
double b = (x-xlo) / h;
double y = a*ya[klo] + b*ya[khi] +
((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
return y;
} // cyc_splint()
static double cyc_lin(double const *xa,
double const *ya,
int n,
double period,
double x)
{
int klo = -1;
int khi = n;
int k;
double xlo = xa[n-1] - period;
double xhi = xa[0] + period;
while (khi-klo > 1) {
k = (khi+klo) >> 1; //(k=(khi+klo)/2)
if (xa[k] > x) {
khi = k;
xhi = xa[k];
}
else {
klo = k;
xlo = xa[k];
}
}
if (khi == n) khi = 0;
if (klo ==-1) klo = n-1;
double h = xhi-xlo;
double a = (xhi-x) / h;
double b = (x-xlo) / h;
double y = a*ya[klo] + b*ya[khi];
return y;
} // cyc_lin()
// cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x,
// with n control points at xa[], ya[], with parameters y2a[].
@ -405,35 +134,26 @@ static double cyc_splintD(double const *xa,
} // cyc_splintD()
// --------------------------------------------
// ------- Calculate the dihedral angle -------
// --------------------------------------------
static const int g_dim=3;
/* ---------------------------------------------------------------------- */
DihedralTableCut::DihedralTableCut(LAMMPS *lmp) : Dihedral(lmp)
DihedralTableCut::DihedralTableCut(LAMMPS *lmp) : DihedralTable(lmp)
{
if (lmp->citeme) lmp->citeme->add(cite_dihedral_tablecut);
ntables = 0;
tables = nullptr;
checkU_fname = checkF_fname = nullptr;
aat_k = aat_theta0_1 = aat_theta0_2 = nullptr;
}
/* ---------------------------------------------------------------------- */
DihedralTableCut::~DihedralTableCut()
{
if (allocated) {
memory->destroy(aat_k);
memory->destroy(aat_theta0_1);
memory->destroy(aat_theta0_2);
for (int m = 0; m < ntables; m++) free_table(&tables[m]);
memory->sfree(tables);
memory->sfree(checkU_fname);
memory->sfree(checkF_fname);
memory->destroy(setflag);
memory->destroy(tabindex);
}
memory->destroy(aat_k);
memory->destroy(aat_theta0_1);
memory->destroy(aat_theta0_2);
}
/* ---------------------------------------------------------------------- */
@ -702,7 +422,7 @@ void DihedralTableCut::compute(int eflag, int vflag)
for (i = 0; i < 4; i++)
for (j = 0; j < 3; j++)
fabcd[i][j] -= - gt*gtt*fpphi*dphidr[i][j]
fabcd[i][j] -= gt*gtt*fpphi*dphidr[i][j]
- gt*gptt*fphi*dthetadr[1][i][j] + gpt*gtt*fphi*dthetadr[0][i][j];
// apply force to each of 4 atoms
@ -751,35 +471,7 @@ void DihedralTableCut::allocate()
memory->create(tabindex,n+1,"dihedral:tabindex");
memory->create(setflag,n+1,"dihedral:setflag");
for (int i = 1; i <= n; i++)
setflag[i] = 0;
}
void DihedralTableCut::settings(int narg, char **arg)
{
if (narg != 2) error->all(FLERR,"Illegal dihedral_style command");
if (strcmp(arg[0],"linear") == 0) tabstyle = LINEAR;
else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
else error->all(FLERR,"Unknown table style in dihedral style table_cut");
tablength = utils::inumeric(FLERR,arg[1],false,lmp);
if (tablength < 3)
error->all(FLERR,"Illegal number of dihedral table entries");
// delete old tables, since cannot just change settings
for (int m = 0; m < ntables; m++) free_table(&tables[m]);
memory->sfree(tables);
if (allocated) {
memory->destroy(setflag);
memory->destroy(tabindex);
}
allocated = 0;
ntables = 0;
tables = nullptr;
for (int i = 1; i <= n; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
@ -790,9 +482,9 @@ void DihedralTableCut::settings(int narg, char **arg)
void DihedralTableCut::coeff(int narg, char **arg)
{
if (narg != 7) error->all(FLERR,"Incorrect args for dihedral coefficients");
if (!allocated) allocate();
int ilo,ihi;
utils::bounds(FLERR,arg[0],1,atom->ndihedraltypes,ilo,ihi,error);
@ -820,29 +512,19 @@ void DihedralTableCut::coeff(int narg, char **arg)
// --- check the angle data for range errors ---
// --- and resolve issues with periodicity ---
if (tb->ninput < 2) {
string err_msg;
err_msg = string("Invalid dihedral table length (")
+ string(arg[5]) + string(").");
error->one(FLERR,err_msg);
}
else if ((tb->ninput == 2) && (tabstyle == SPLINE)) {
string err_msg;
err_msg = string("Invalid dihedral spline table length. (Try linear)\n (")
+ string(arg[5]) + string(").");
error->one(FLERR,err_msg);
}
if (tb->ninput < 2)
error->all(FLERR,fmt::format("Invalid dihedral table length: {}",arg[5]));
else if ((tb->ninput == 2) && (tabstyle == SPLINE))
error->all(FLERR,fmt::format("Invalid dihedral spline table length: {} "
"(Try linear)",arg[5]));
// check for monotonicity
for (int i=0; i < tb->ninput-1; i++) {
if (tb->phifile[i] >= tb->phifile[i+1]) {
stringstream i_str;
i_str << i+1;
string err_msg =
string("Dihedral table values are not increasing (") +
string(arg[5]) + string(", ")+i_str.str()+string("th entry)");
auto err_msg =fmt::format("Dihedral table values are not increasing "
"({}, entry {})",arg[5],i+1);
if (i==0)
err_msg += string("\n(This is probably a mistake with your table format.)\n");
err_msg += "\n(This is probably a mistake with your table format.)\n";
error->all(FLERR,err_msg);
}
}
@ -851,20 +533,13 @@ void DihedralTableCut::coeff(int narg, char **arg)
double philo = tb->phifile[0];
double phihi = tb->phifile[tb->ninput-1];
if (tb->use_degrees) {
if ((phihi - philo) >= 360) {
string err_msg;
err_msg = string("Dihedral table angle range must be < 360 degrees (")
+string(arg[5]) + string(").");
error->all(FLERR,err_msg);
}
}
else {
if ((phihi - philo) >= MY_2PI) {
string err_msg;
err_msg = string("Dihedral table angle range must be < 2*PI radians (")
+ string(arg[5]) + string(").");
error->all(FLERR,err_msg);
}
if ((phihi - philo) >= 360)
error->all(FLERR,fmt::format("Dihedral table angle range must be < 360 "
"degrees ({})",arg[5]));
} else {
if ((phihi - philo) >= MY_2PI)
error->all(FLERR,fmt::format("Dihedral table angle range must be < 2*PI "
"radians ({})",arg[5]));
}
// convert phi from degrees to radians
@ -932,10 +607,9 @@ void DihedralTableCut::coeff(int narg, char **arg)
// Optional: allow the user to print out the interpolated spline tables
if (me == 0) {
if (checkU_fname && (strlen(checkU_fname) != 0))
{
ofstream checkU_file;
checkU_file.open(checkU_fname, ios::out);
if (!checkU_fname.empty()) {
std::ofstream checkU_file;
checkU_file.open(checkU_fname, std::ios::out);
for (int i=0; i < tablength; i++) {
double phi = i*MY_2PI/tablength;
double u = tb->e[i];
@ -945,12 +619,10 @@ void DihedralTableCut::coeff(int narg, char **arg)
}
checkU_file.close();
}
if (checkF_fname && (strlen(checkF_fname) != 0))
{
ofstream checkF_file;
checkF_file.open(checkF_fname, ios::out);
for (int i=0; i < tablength; i++)
{
if (!checkF_fname.empty()) {
std::ofstream checkF_file;
checkF_file.open(checkF_fname, std::ios::out);
for (int i=0; i < tablength; i++) {
double phi = i*MY_2PI/tablength;
double f;
if ((tabstyle == SPLINE) && (tb->f_unspecified)) {
@ -992,429 +664,3 @@ void DihedralTableCut::coeff(int narg, char **arg)
if (count == 0) error->all(FLERR,"Incorrect args for dihedral coefficients");
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void DihedralTableCut::write_restart(FILE *fp)
{
write_restart_settings(fp);
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
void DihedralTableCut::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void DihedralTableCut::write_restart_settings(FILE *fp)
{
fwrite(&tabstyle,sizeof(int),1,fp);
fwrite(&tablength,sizeof(int),1,fp);
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
void DihedralTableCut::read_restart_settings(FILE *fp)
{
if (comm->me == 0) {
utils::sfread(FLERR,&tabstyle,sizeof(int),1,fp,nullptr,error);
utils::sfread(FLERR,&tablength,sizeof(int),1,fp,nullptr,error);
}
MPI_Bcast(&tabstyle,1,MPI_INT,0,world);
MPI_Bcast(&tablength,1,MPI_INT,0,world);
}
/* ---------------------------------------------------------------------- */
void DihedralTableCut::null_table(Table *tb)
{
tb->phifile = tb->efile = tb->ffile = nullptr;
tb->e2file = tb->f2file = nullptr;
tb->phi = tb->e = tb->de = nullptr;
tb->f = tb->df = tb->e2 = tb->f2 = nullptr;
}
/* ---------------------------------------------------------------------- */
void DihedralTableCut::free_table(Table *tb)
{
memory->destroy(tb->phifile);
memory->destroy(tb->efile);
memory->destroy(tb->ffile);
memory->destroy(tb->e2file);
memory->destroy(tb->f2file);
memory->destroy(tb->phi);
memory->destroy(tb->e);
memory->destroy(tb->de);
memory->destroy(tb->f);
memory->destroy(tb->df);
memory->destroy(tb->e2);
memory->destroy(tb->f2);
}
/* ----------------------------------------------------------------------
read table file, only called by proc 0
------------------------------------------------------------------------- */
static const int MAXLINE=2048;
void DihedralTableCut::read_table(Table *tb, char *file, char *keyword)
{
char line[MAXLINE];
// open file
FILE *fp = utils::open_potential(file,lmp,nullptr);
if (fp == nullptr) {
string err_msg = string("Cannot open file ") + string(file);
error->one(FLERR,err_msg);
}
// loop until section found with matching keyword
while (1) {
if (fgets(line,MAXLINE,fp) == nullptr) {
string err_msg=string("Did not find keyword \"")
+string(keyword)+string("\" in dihedral table file.");
error->one(FLERR, err_msg);
}
if (strspn(line," \t\n\r") == strlen(line)) continue; // blank line
if (line[0] == '#') continue; // comment
char *word = strtok(line," \t\n\r");
if (strcmp(word,keyword) == 0) break; // matching keyword
utils::sfgets(FLERR,line,MAXLINE,fp,file,error); // no match, skip section
param_extract(tb,line);
utils::sfgets(FLERR,line,MAXLINE,fp,file,error);
for (int i = 0; i < tb->ninput; i++)
utils::sfgets(FLERR,line,MAXLINE,fp,file,error);
}
// read args on 2nd line of section
// allocate table arrays for file values
utils::sfgets(FLERR,line,MAXLINE,fp,file,error);
param_extract(tb,line);
memory->create(tb->phifile,tb->ninput,"dihedral:phifile");
memory->create(tb->efile,tb->ninput,"dihedral:efile");
memory->create(tb->ffile,tb->ninput,"dihedral:ffile");
// read a,e,f table values from file
int itmp;
for (int i = 0; i < tb->ninput; i++) {
// Read the next line. Make sure the file is long enough.
utils::sfgets(FLERR,line,MAXLINE,fp,file,error);
// Skip blank lines and delete text following a '#' character
char *pe = strchr(line, '#');
if (pe != nullptr) *pe = '\0'; //terminate string at '#' character
char *pc = line;
while ((*pc != '\0') && isspace(*pc))
pc++;
if (*pc != '\0') { //If line is not a blank line
stringstream line_ss(line);
if (tb->f_unspecified) {
line_ss >> itmp;
line_ss >> tb->phifile[i];
line_ss >> tb->efile[i];
} else {
line_ss >> itmp;
line_ss >> tb->phifile[i];
line_ss >> tb->efile[i];
line_ss >> tb->ffile[i];
}
if (! line_ss) {
stringstream err_msg;
err_msg << "Read error in table "<< keyword<<", near line "<<i+1<<"\n"
<< " (Check to make sure the number of columns is correct.)";
if ((! tb->f_unspecified) && (i==0))
err_msg << "\n (This sometimes occurs if users forget to specify the \"NOF\" option.)\n";
error->one(FLERR, err_msg.str().c_str());
}
} else //if it is a blank line, then skip it.
i--;
} //for (int i = 0; (i < tb->ninput) && fp; i++) {
fclose(fp);
}
/* ----------------------------------------------------------------------
build spline representation of e,f over entire range of read-in table
this function sets these values in e2file,f2file.
I also perform a crude check for force & energy consistency.
------------------------------------------------------------------------- */
void DihedralTableCut::spline_table(Table *tb)
{
memory->create(tb->e2file,tb->ninput,"dihedral:e2file");
memory->create(tb->f2file,tb->ninput,"dihedral:f2file");
if (cyc_spline(tb->phifile, tb->efile, tb->ninput,
MY_2PI,tb->e2file,comm->me == 0))
error->one(FLERR,"Error computing dihedral spline tables");
if (! tb->f_unspecified) {
if (cyc_spline(tb->phifile, tb->ffile, tb->ninput,
MY_2PI, tb->f2file, comm->me == 0))
error->one(FLERR,"Error computing dihedral spline tables");
}
// CHECK to help make sure the user calculated forces in a way
// which is grossly numerically consistent with the energy table.
if (! tb->f_unspecified) {
int num_disagreements = 0;
for (int i=0; i<tb->ninput; i++) {
// Calculate what the force should be at the control points
// by using linear interpolation of the derivatives of the energy:
double phi_i = tb->phifile[i];
// First deal with periodicity
double phi_im1, phi_ip1;
int im1 = i-1;
if (im1 < 0) {
im1 += tb->ninput;
phi_im1 = tb->phifile[im1] - MY_2PI;
}
else
phi_im1 = tb->phifile[im1];
int ip1 = i+1;
if (ip1 >= tb->ninput) {
ip1 -= tb->ninput;
phi_ip1 = tb->phifile[ip1] + MY_2PI;
}
else
phi_ip1 = tb->phifile[ip1];
// Now calculate the midpoints above and below phi_i = tb->phifile[i]
double phi_lo= 0.5*(phi_im1 + phi_i); //midpoint between phi_im1 and phi_i
double phi_hi= 0.5*(phi_i + phi_ip1); //midpoint between phi_i and phi_ip1
// Use a linear approximation to the derivative at these two midpoints
double dU_dphi_lo =
(tb->efile[i] - tb->efile[im1])
/
(phi_i - phi_im1);
double dU_dphi_hi =
(tb->efile[ip1] - tb->efile[i])
/
(phi_ip1 - phi_i);
// Now calculate the derivative at position
// phi_i (=tb->phifile[i]) using linear interpolation
double a = (phi_i - phi_lo) / (phi_hi - phi_lo);
double b = (phi_hi - phi_i) / (phi_hi - phi_lo);
double dU_dphi = a*dU_dphi_lo + b*dU_dphi_hi;
double f = -dU_dphi;
// Alternately, we could use spline interpolation instead:
// double f = - splintD(tb->phifile, tb->efile, tb->e2file,
// tb->ninput, MY_2PI, tb->phifile[i]);
// This is the way I originally did it, but I trust
// the ugly simple linear way above better.
// Recall this entire block of code doess not calculate
// anything important. It does not have to be perfect.
// We are only checking for stupid user errors here.
if ((f != 0.0) &&
(tb->ffile[i] != 0.0) &&
((f/tb->ffile[i] < 0.5) || (f/tb->ffile[i] > 2.0))) {
num_disagreements++;
}
} // for (int i=0; i<tb->ninput; i++)
if ((num_disagreements > tb->ninput/2) && (num_disagreements > 2)) {
string msg("Dihedral table has inconsistent forces and energies. (Try \"NOF\".)\n");
error->all(FLERR,msg.c_str());
}
} // check for consistency if (! tb->f_unspecified)
} // DihedralTable::spline_table()
/* ----------------------------------------------------------------------
compute a,e,f vectors from splined values
------------------------------------------------------------------------- */
void DihedralTableCut::compute_table(Table *tb)
{
//delta = table spacing in dihedral angle for tablength (cyclic) bins
tb->delta = MY_2PI / tablength;
tb->invdelta = 1.0/tb->delta;
tb->deltasq6 = tb->delta*tb->delta / 6.0;
// N evenly spaced bins in dihedral angle from 0 to 2*PI
// phi,e,f = value at lower edge of bin
// de,df values = delta values of e,f (cyclic, in this case)
// phi,e,f,de,df are arrays containing "tablength" number of entries
memory->create(tb->phi,tablength,"dihedral:phi");
memory->create(tb->e,tablength,"dihedral:e");
memory->create(tb->de,tablength,"dihedral:de");
memory->create(tb->f,tablength,"dihedral:f");
memory->create(tb->df,tablength,"dihedral:df");
memory->create(tb->e2,tablength,"dihedral:e2");
memory->create(tb->f2,tablength,"dihedral:f2");
if (tabstyle == SPLINE) {
// Use cubic spline interpolation to calculate the entries in the
// internal table. (This is true regardless...even if tabstyle!=SPLINE.)
for (int i = 0; i < tablength; i++) {
double phi = i*tb->delta;
tb->phi[i] = phi;
tb->e[i]= cyc_splint(tb->phifile,tb->efile,tb->e2file,tb->ninput,MY_2PI,phi);
if (! tb->f_unspecified)
tb->f[i] = cyc_splint(tb->phifile,tb->ffile,tb->f2file,tb->ninput,MY_2PI,phi);
}
} // if (tabstyle == SPLINE)
else if (tabstyle == LINEAR) {
if (! tb->f_unspecified) {
for (int i = 0; i < tablength; i++) {
double phi = i*tb->delta;
tb->phi[i] = phi;
tb->e[i]= cyc_lin(tb->phifile,tb->efile,tb->ninput,MY_2PI,phi);
tb->f[i]= cyc_lin(tb->phifile,tb->ffile,tb->ninput,MY_2PI,phi);
}
}
else {
for (int i = 0; i < tablength; i++) {
double phi = i*tb->delta;
tb->phi[i] = phi;
tb->e[i]= cyc_lin(tb->phifile,tb->efile,tb->ninput,MY_2PI,phi);
}
// In the linear case, if the user did not specify the forces, then we
// must generate the "f" array. Do this using linear interpolation
// of the e array (which itself was generated above)
for (int i = 0; i < tablength; i++) {
int im1 = i-1; if (im1 < 0) im1 += tablength;
int ip1 = i+1; if (ip1 >= tablength) ip1 -= tablength;
double dedx = (tb->e[ip1] - tb->e[im1]) / (2.0 * tb->delta);
// (This is the average of the linear slopes on either side of the node.
// Note that the nodes in the internal table are evenly spaced.)
tb->f[i] = -dedx;
}
}
// Fill the linear interpolation tables (de, df)
for (int i = 0; i < tablength; i++) {
int ip1 = i+1; if (ip1 >= tablength) ip1 -= tablength;
tb->de[i] = tb->e[ip1] - tb->e[i];
tb->df[i] = tb->f[ip1] - tb->f[i];
}
} // else if (tabstyle == LINEAR)
cyc_spline(tb->phi, tb->e, tablength, MY_2PI, tb->e2, comm->me == 0);
if (! tb->f_unspecified)
cyc_spline(tb->phi, tb->f, tablength, MY_2PI, tb->f2, comm->me == 0);
}
/* ----------------------------------------------------------------------
extract attributes from parameter line in table section
format of line: N value NOF DEGREES RADIANS
N is required, other params are optional
------------------------------------------------------------------------- */
void DihedralTableCut::param_extract(Table *tb, char *line)
{
//tb->theta0 = 180.0; <- equilibrium angles not supported
tb->ninput = 0;
tb->f_unspecified = false; //default
tb->use_degrees = true; //default
char *word = strtok(line," \t\n\r\f");
while (word) {
if (strcmp(word,"N") == 0) {
word = strtok(nullptr," \t\n\r\f");
tb->ninput = atoi(word);
}
else if (strcmp(word,"NOF") == 0) {
tb->f_unspecified = true;
}
else if ((strcmp(word,"DEGREES") == 0) || (strcmp(word,"degrees") == 0)) {
tb->use_degrees = true;
}
else if ((strcmp(word,"RADIANS") == 0) || (strcmp(word,"radians") == 0)) {
tb->use_degrees = false;
}
else if (strcmp(word,"CHECKU") == 0) {
word = strtok(nullptr," \t\n\r\f");
memory->sfree(checkU_fname);
memory->create(checkU_fname,strlen(word)+1,"dihedral_table:checkU");
strcpy(checkU_fname, word);
}
else if (strcmp(word,"CHECKF") == 0) {
word = strtok(nullptr," \t\n\r\f");
memory->sfree(checkF_fname);
memory->create(checkF_fname,strlen(word)+1,"dihedral_table:checkF");
strcpy(checkF_fname, word);
}
// COMMENTING OUT: equilibrium angles are not supported
//else if (strcmp(word,"EQ") == 0) {
// word = strtok(nullptr," \t\n\r\f");
// tb->theta0 = atof(word);
//}
else {
string err_msg("Invalid keyword in dihedral angle table parameters");
err_msg += string(" (") + string(word) + string(")");
error->one(FLERR, err_msg);
}
word = strtok(nullptr," \t\n\r\f");
}
if (tb->ninput == 0)
error->one(FLERR,"Dihedral table parameters did not set N");
} // DihedralTable::param_extract()
/* ----------------------------------------------------------------------
broadcast read-in table info from proc 0 to other procs
this function communicates these values in Table:
ninput,phifile,efile,ffile,
f_unspecified,use_degrees
------------------------------------------------------------------------- */
void DihedralTableCut::bcast_table(Table *tb)
{
MPI_Bcast(&tb->ninput,1,MPI_INT,0,world);
int me;
MPI_Comm_rank(world,&me);
if (me > 0) {
memory->create(tb->phifile,tb->ninput,"dihedral:phifile");
memory->create(tb->efile,tb->ninput,"dihedral:efile");
memory->create(tb->ffile,tb->ninput,"dihedral:ffile");
}
MPI_Bcast(tb->phifile,tb->ninput,MPI_DOUBLE,0,world);
MPI_Bcast(tb->efile,tb->ninput,MPI_DOUBLE,0,world);
MPI_Bcast(tb->ffile,tb->ninput,MPI_DOUBLE,0,world);
MPI_Bcast(&tb->f_unspecified,1,MPI_INT,0,world);
MPI_Bcast(&tb->use_degrees,1,MPI_INT,0,world);
// COMMENTING OUT: equilibrium angles are not supported
//MPI_Bcast(&tb->theta0,1,MPI_DOUBLE,0,world);
}

118
src/USER-MISC/dihedral_table_cut.h
View File

@ -20,131 +20,21 @@ DihedralStyle(table/cut,DihedralTableCut)
#ifndef LMP_DIHEDRAL_TABLE_CUT_H
#define LMP_DIHEDRAL_TABLE_CUT_H
#include "dihedral.h"
#include "dihedral_table.h"
namespace LAMMPS_NS {
class DihedralTableCut : public Dihedral {
class DihedralTableCut : public DihedralTable {
public:
DihedralTableCut(class LAMMPS *);
virtual ~DihedralTableCut();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
virtual void coeff(int, char **);
protected:
double *aat_k,*aat_theta0_1,*aat_theta0_2;
void allocate();
int tabstyle,tablength;
char *checkU_fname;
char *checkF_fname;
struct Table {
int ninput;
int f_unspecified; // boolean (but MPI does not like type "bool")
int use_degrees; // boolean (but MPI does not like type "bool")
double *phifile,*efile,*ffile;
double *e2file,*f2file;
double delta,invdelta,deltasq6;
double *phi,*e,*de,*f,*df,*e2,*f2;
};
int ntables;
Table *tables;
int *tabindex;
void null_table(Table *);
void free_table(Table *);
void read_table(Table *, char *, char *);
void bcast_table(Table *);
void spline_table(Table *);
void compute_table(Table *);
void param_extract(Table *, char *);
// --------------------------------------------
// ------------ inline functions --------------
// --------------------------------------------
// -----------------------------------------------------------
// uf_lookup()
// quickly calculate the potential u and force f at angle x,
// using the internal tables tb->e and tb->f (evenly spaced)
// -----------------------------------------------------------
enum{LINEAR,SPLINE};
inline void uf_lookup(int type, double x, double &u, double &f)
{
Table *tb = &tables[tabindex[type]];
double x_over_delta = x*tb->invdelta;
int i = static_cast<int> (x_over_delta);
double a;
double b = x_over_delta - i;
// Apply periodic boundary conditions to indices i and i+1
if (i >= tablength) i -= tablength;
int ip1 = i+1; if (ip1 >= tablength) ip1 -= tablength;
switch(tabstyle) {
case LINEAR:
u = tb->e[i] + b * tb->de[i];
f = -(tb->f[i] + b * tb->df[i]); //<--works even if tb->f_unspecified==true
break;
case SPLINE:
a = 1.0 - b;
u = a * tb->e[i] + b * tb->e[ip1] +
((a*a*a-a)*tb->e2[i] + (b*b*b-b)*tb->e2[ip1]) *
tb->deltasq6;
if (tb->f_unspecified)
//Formula below taken from equation3.3.5 of "numerical recipes in c"
//"f"=-derivative of e with respect to x (or "phi" in this case)
f = -((tb->e[i]-tb->e[ip1])*tb->invdelta +
((3.0*a*a-1.0)*tb->e2[i]+(1.0-3.0*b*b)*tb->e2[ip1])*tb->delta/6.0);
else
f = -(a * tb->f[i] + b * tb->f[ip1] +
((a*a*a-a)*tb->f2[i] + (b*b*b-b)*tb->f2[ip1]) *
tb->deltasq6);
break;
} // switch(tabstyle)
} // uf_lookup()
// ----------------------------------------------------------
// u_lookup()
// quickly calculate the potential u at angle x using tb->e
//-----------------------------------------------------------
inline void u_lookup(int type, double x, double &u)
{
Table *tb = &tables[tabindex[type]];
int N = tablength;
// i = static_cast<int> ((x - tb->lo) * tb->invdelta); <-general version
double x_over_delta = x*tb->invdelta;
int i = static_cast<int> (x_over_delta);
double b = x_over_delta - i;
// Apply periodic boundary conditions to indices i and i+1
if (i >= N) i -= N;
int ip1 = i+1; if (ip1 >= N) ip1 -= N;
if (tabstyle == LINEAR) {
u = tb->e[i] + b * tb->de[i];
}
else if (tabstyle == SPLINE) {
double a = 1.0 - b;
u = a * tb->e[i] + b * tb->e[ip1] +
((a*a*a-a)*tb->e2[i] + (b*b*b-b)*tb->e2[ip1]) *
tb->deltasq6;
}
} // u_lookup()
virtual void allocate();
};
}

Some files were not shown because too many files have changed in this diff Show More