ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source 
Collected small changes and fixes
This commit is contained in:
Axel Kohlmeyer
2025-01-09 19:33:52 -05:00
committed by GitHub
parent cb930d6470 a769b3d6b2
commit 4330801a2f
9 changed files with 375 additions and 209 deletions
Download Patch File Download Diff File Expand all files Collapse all files

72
doc/src/Developer_unittest.rst
View File

@ -227,12 +227,12 @@ Tests for the C-style library interface
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Tests for validating the LAMMPS C-style library interface are in the
``unittest/c-library`` folder. They are implemented either to be used
for utility functions or for LAMMPS commands, but use the functions
implemented in the ``src/library.cpp`` file as much as possible. There
may be some overlap with other tests, but only in as much as is required
to test the C-style library API. The tests are distributed over
multiple test programs which try to match the grouping of the
``unittest/c-library`` folder. They text either utility functions or
LAMMPS commands, but use the functions implemented in
``src/library.cpp`` as much as possible. There may be some overlap with
other tests as far as the LAMMPS functionality is concerned, but the
focus is on testing the C-style library API. The tests are distributed
over multiple test programs which try to match the grouping of the
functions in the source code and :ref:`in the manual <lammps_c_api>`.
This group of tests also includes tests invoking LAMMPS in parallel
@ -258,7 +258,7 @@ Tests for the Python module and package
The ``unittest/python`` folder contains primarily tests for classes and
functions in the LAMMPS python module but also for commands in the
PYTHON package. These tests are only enabled if the necessary
PYTHON package. These tests are only enabled, if the necessary
prerequisites are detected or enabled during configuration and
compilation of LAMMPS (shared library build enabled, Python interpreter
found, Python development files found).
@ -272,29 +272,30 @@ Tests for the Fortran interface
Tests for using the Fortran module are in the ``unittest/fortran``
folder. Since they are also using the GoogleTest library, they require
implementing test wrappers in C++ that will call fortran functions
which provide a C function interface through ISO_C_BINDINGS that will in
turn call the functions in the LAMMPS Fortran module.
test wrappers written in C++ that will call fortran functions with a C
function interface through ISO_C_BINDINGS which will in turn call the
functions in the LAMMPS Fortran module.
Tests for the C++-style library interface
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The tests in the ``unittest/cplusplus`` folder are somewhat similar to
the tests for the C-style library interface, but do not need to test the
several convenience and utility functions that are only available through
the C-style interface. Instead it can focus on the more generic features
that are used internally. This part of the unit tests is currently still
mostly in the planning stage.
convenience and utility functions that are only available through the
C-style library interface. Instead they focus on the more generic
features that are used in LAMMPS internally. This part of the unit
tests is currently still mostly in the planning stage.
Tests for reading and writing file formats
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The ``unittest/formats`` folder contains test programs for reading and
writing files like data files, restart files, potential files or dump files.
This covers simple things like the file i/o convenience functions in the
``utils::`` namespace to complex tests of atom styles where creating and
deleting atoms with different properties is tested in different ways
and through script commands or reading and writing of data or restart files.
writing files like data files, restart files, potential files or dump
files. This covers simple things like the file i/o convenience
functions in the ``utils::`` namespace to complex tests of atom styles
where creating and deleting of atoms with different properties is tested
in different ways and through script commands or reading and writing of
data or restart files.
Tests for styles computing or modifying forces
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@ -443,7 +444,7 @@ file for a style that is similar to one to be tested. The file name should
follow the naming conventions described above and after copying the file,
the first step is to replace the style names where needed. The coefficient
values do not have to be meaningful, just in a reasonable range for the
given system. It does not matter if some forces are large, as long as
given system. It does not matter if some forces are large, for as long as
they do not diverge.
The template input files define a large number of index variables at the top
@ -531,19 +532,20 @@ Python module.
Troubleshooting failed unit tests
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
The are by default no unit tests for newly added features (e.g. pair, fix,
or compute styles) unless your pull request also includes tests for the
added features. If you are modifying some features, you may see failures
for existing tests, if your modifications have some unexpected side effects
or your changes render the existing test invalid. If you are adding an
accelerated version of an existing style, then only tests for INTEL,
KOKKOS (with OpenMP only), OPENMP, and OPT will be run automatically.
Tests for the GPU package are time consuming and thus are only run
*after* a merge, or when a special label, ``gpu_unit_tests`` is added
to the pull request. After the test has started, it is often best to
remove the label since every PR activity will re-trigger the test (that
is a limitation of triggering a test with a label). Support for unit
tests when using KOKKOS with GPU acceleration is currently not supported.
There are by default no unit tests for newly added features (e.g. pair,
fix, or compute styles) unless your pull request also includes tests for
these added features. If you are modifying some existing LAMMPS
features, you may see failures for existing tests, if your modifications
have some unexpected side effects or your changes render the existing
test invalid. If you are adding an accelerated version of an existing
style, then only tests for INTEL, KOKKOS (with OpenMP only), OPENMP, and
OPT will be run automatically. Tests for the GPU package are time
consuming and thus are only run *after* a merge, or when a special
label, ``gpu_unit_tests`` is added to the pull request. After the test
has started, it is often best to remove the label since every PR
activity will re-trigger the test (that is a limitation of triggering a
test with a label). Support for unit tests using KOKKOS with GPU
acceleration is currently not supported.
When you see a failed build on GitHub, click on ``Details`` to be taken
to the corresponding LAMMPS Jenkins CI web page. Click on the "Exit"
@ -588,7 +590,7 @@ While the epsilon (relative precision) for a single, `IEEE 754 compliant
<https://en.wikipedia.org/wiki/IEEE_754>`_, double precision floating
point operation is at about 2.2e-16, the achievable precision for the
tests is lower due to most numbers being sums over intermediate results
and the non-associativity of floating point math leading to larger
for which the non-associativity of floating point math leads to larger
errors. As a rule of thumb, the test epsilon can often be in the range
5.0e-14 to 1.0e-13. But for "noisy" force kernels, e.g. those a larger
amount of arithmetic operations involving `exp()`, `log()` or `sin()`
@ -602,7 +604,7 @@ of floating point operations or that some or most intermediate operations
may be done using approximations or with single precision floating point
math.
To rerun the failed unit test individually, change to the ``build`` directory
To rerun a failed unit test individually, change to the ``build`` directory
and run the test with verbose output. For example,
.. code-block:: bash

50
doc/src/Errors_debug.rst
View File

@ -235,3 +235,53 @@ from GDB. In addition you get a more specific hint about what cause the
segmentation fault, i.e. that it is a NULL pointer dereference. To find
out which pointer exactly was NULL, you need to use the debugger, though.
Debugging when LAMMPS appears to be stuck
=========================================
Sometimes the LAMMPS calculation appears to be stuck, that is the LAMMPS
process or processes are active, but there is no visible progress. This
can have multiple reasons:
- The selected styles are slow and require a lot of CPU time and the
system is large. When extrapolating the expected speed from smaller
systems, one has to factor in that not all models scale linearly with
system size, e.g. :doc:`kspace styles like ewald or pppm
<kspace_style>`. There is very little that can be done in this case.
- The output interval is not set or set to a large value with the
:doc:`thermo <thermo>` command. I the first case, there will be output
only at the first and last step.
- The output is block-buffered and instead of line-buffered. The output
will only be written to the screen after 4096 or 8192 characters of
output have accumulated. This most often happens for files but also
with MPI parallel executables for output to the screen, since the
output to the screen is handled by the MPI library so that output from
all processes can be shown. This can be suppressed by using the
``-nonblock`` or ``-nb`` command-line flag, which turns off buffering
for screen and logfile output.
- An MPI parallel calculation has a bug where a collective MPI function
is called (e.g. ``MPI_Barrier()``, ``MPI_Bcast()``,
``MPI_Allreduce()`` and so on) before pending point-to-point
communications are completed or when the collective function is only
called from a subset of the MPI processes. This also applies to some
internal LAMMPS functions like ``Error::all()`` which uses
``MPI_Barrier()`` and thus ``Error::one()`` must be called, if the
error condition does not happen on all MPI processes simultaneously.
- Some function in LAMMPS has a bug where a ``for`` or ``while`` loop
does not trigger the exit condition and thus will loop forever. This
can happen when the wrong variable is incremented or when one value in
a comparison becomes ``NaN`` due to an overflow.
In the latter two cases, further information and stack traces (see above)
can be obtain by attaching a debugger to a running process. For that the
process ID (PID) is needed; this can be found on Linux machines with the
``top``, ``htop``, ``ps``, or ``pstree`` commands.
Then running the (GNU) debugger ``gdb`` with the ``-p`` flag followed by
the process id will attach the process to the debugger and stop
execution of that specific process. From there on it is possible to
issue all debugger commands in the same way as when LAMMPS was started
from the debugger (see above). Most importantly it is possible to
obtain a stack trace with the ``where`` command and thus determine where
in the execution of a timestep this process is. Also internal data can
be printed and execution single stepped or continued. When the debugger
is exited, the calculation will resume normally.

6
doc/src/variable.rst
View File

@ -71,9 +71,9 @@ Syntax
feature functions = is_available(category,feature), is_active(category,feature), is_defined(category,id)
atom value = id[i], mass[i], type[i], mol[i], x[i], y[i], z[i], vx[i], vy[i], vz[i], fx[i], fy[i], fz[i], q[i]
atom vector = id, mass, type, mol, radius, q, x, y, z, vx, vy, vz, fx, fy, fz
custom atom property = i_name, d_name, i_name[i], d_name[i], i2_name[i], d2_name[i], i2_name[i][j], d_name[i][j]
compute references = c_ID, c_ID[i], c_ID[i][j], C_ID, C_ID[i]
fix references = f_ID, f_ID[i], f_ID[i][j], F_ID, F_ID[i]
custom atom property = i_name, d_name, i_name[i], d_name[i], i2_name[i], d2_name[i], i2_name[i][j], d2_name[i][j]
compute references = c_ID, c_ID[i], c_ID[i][j], C_ID, C_ID[i], C_ID[i][j]
fix references = f_ID, f_ID[i], f_ID[i][j], F_ID, F_ID[i], F_ID[i][j]
variable references = v_name, v_name[i]
vector initialization = [1,3,7,10] (for *vector* variables only)

7
src/MISC/fix_ipi.cpp
View File

@ -438,7 +438,7 @@ void FixIPI::final_integrate()
char header[MSGLEN+1];
double vir[9], pot=0.0;
double forceconv, potconv, posconv, pressconv, posconv3;
char retstr[1024];
char retstr[1024] = { '\0' };
// conversions from LAMMPS units to atomic units, which are used by i-PI
potconv=3.1668152e-06/force->boltz;
@ -488,7 +488,7 @@ void FixIPI::final_integrate()
vir[1] = comp_p->vector[3]*pressconv*myvol;
vir[2] = comp_p->vector[4]*pressconv*myvol;
vir[5] = comp_p->vector[5]*pressconv*myvol;
retstr[0]=0;
retstr[0] = '\0';
}
if (master) {
@ -511,7 +511,8 @@ void FixIPI::final_integrate()
writebuffer(ipisock,(char*) &nat,4, error);
writebuffer(ipisock,(char*) buffer, bsize*8, error);
writebuffer(ipisock,(char*) vir,9*8, error);
nat=strlen(retstr); writebuffer(ipisock,(char*) &nat,4, error);
nat=strlen(retstr);
writebuffer(ipisock,(char*) &nat,4, error);
writebuffer(ipisock,(char*) retstr, nat, error);
}
else

4
src/NETCDF/dump_netcdf_mpiio.cpp
View File

@ -639,8 +639,8 @@ void DumpNetCDFMPIIO::write()
nme = count();
int *block_sizes = new int[comm->nprocs];
MPI_Allgather(&nme, 1, MPI_INT, block_sizes, 1, MPI_INT, world);
blocki = 0;
for (int i = 0; i < comm->me; i++) blocki += block_sizes[i];
blocki = (MPI_Offset) 0;
for (int i = 0; i < comm->me; i++) blocki += (MPI_Offset) (block_sizes[i]);
delete[] block_sizes;
// ensure buf is sized for packing and communicating

6
src/NETCDF/dump_netcdf_mpiio.h
View File

@ -49,11 +49,11 @@ class DumpNetCDFMPIIO : public DumpCustom {
int quantity; // type of the quantity
};
int framei; // current frame index
int blocki; // current block index
MPI_Offset framei; // current frame index
MPI_Offset blocki; // current block index
int ndata; // number of data blocks to expect
bigint ntotalgr; // # of atoms
MPI_Offset ntotalgr; // # of atoms
int n_perat; // # of netcdf per-atom properties
nc_perat_t *perat; // per-atom properties

3
src/REAXFF/fix_reaxff_species.cpp
View File

@ -261,7 +261,8 @@ FixReaxFFSpecies::FixReaxFFSpecies(LAMMPS *lmp, int narg, char **arg) :
posflag = 1;
posfreq = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
if (posfreq < nfreq || (posfreq % nfreq != 0))
error->all(FLERR, "Incompatible fix reaxff/species postion frequency {}", posfreq);
error->all(FLERR, "Incompatible fix reaxff/species position frequency {}",
posfreq);
filepos = new char[255];
strcpy(filepos, arg[iarg + 2]);

288
src/VTK/dump_vtk.cpp
View File

@ -107,7 +107,7 @@ static constexpr int DELTA = 1048576;
DumpVTK::DumpVTK(LAMMPS *lmp, int narg, char **arg) :
DumpCustom(lmp, narg, arg)
{
if (narg == 5) error->all(FLERR,"No dump vtk arguments specified");
if (narg < 6) utils::missing_cmd_args(FLERR,"dump vtk", error);
pack_choice.clear();
vtype.clear();
@ -214,17 +214,18 @@ void DumpVTK::init_style()
// check that fix frequency is acceptable
for (int i = 0; i < ncompute; i++) {
int icompute = modify->find_compute(id_compute[i]);
if (icompute < 0) error->all(FLERR,"Could not find dump vtk compute ID");
compute[i] = modify->compute[icompute];
compute[i] = modify->get_compute_by_id(id_compute[i]);
if (!compute[i]) error->all(FLERR,"Could not find dump vtk compute ID {}", id_compute[i]);
}
for (int i = 0; i < nfix; i++) {
int ifix = modify->find_fix(id_fix[i]);
if (ifix < 0) error->all(FLERR,"Could not find dump vtk fix ID");
fix[i] = modify->fix[ifix];
if (nevery % modify->fix[ifix]->peratom_freq)
error->all(FLERR,"Dump vtk and fix not computed at compatible times");
fix[i] = modify->get_fix_by_id(id_fix[i]);
if (!fix[i]) {
error->all(FLERR,"Could not find dump vtk fix ID {}", id_fix[i]);
} else {
if (nevery % fix[i]->peratom_freq)
error->all(FLERR,"Dump vtk and fix ID {} not called at compatible times", id_fix[i]);
}
}
for (int i = 0; i < nvariable; i++) {
@ -359,8 +360,7 @@ int DumpVTK::count()
nstride = 1;
} else if (thresh_array[ithresh] == MOL) {
if (!atom->molecule_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR, "Threshold for an atom property that isn't allocated");
tagint *molecule = atom->molecule;
for (i = 0; i < nlocal; i++) dchoose[i] = molecule[i];
ptr = dchoose;
@ -654,64 +654,54 @@ int DumpVTK::count()
} else if (thresh_array[ithresh] == RADIUS) {
if (!atom->radius_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR,"Threshold for an atom property that isn't allocated");
ptr = atom->radius;
nstride = 1;
} else if (thresh_array[ithresh] == DIAMETER) {
if (!atom->radius_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR,"Threshold for an atom property that isn't allocated");
double *radius = atom->radius;
for (i = 0; i < nlocal; i++) dchoose[i] = 2.0*radius[i];
ptr = dchoose;
nstride = 1;
} else if (thresh_array[ithresh] == OMEGAX) {
if (!atom->omega_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR,"Threshold for an atom property that isn't allocated");
ptr = &atom->omega[0][0];
nstride = 3;
} else if (thresh_array[ithresh] == OMEGAY) {
if (!atom->omega_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR,"Threshold for an atom property that isn't allocated");
ptr = &atom->omega[0][1];
nstride = 3;
} else if (thresh_array[ithresh] == OMEGAZ) {
if (!atom->omega_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR,"Threshold for an atom property that isn't allocated");
ptr = &atom->omega[0][2];
nstride = 3;
} else if (thresh_array[ithresh] == ANGMOMX) {
if (!atom->angmom_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR,"Threshold for an atom property that isn't allocated");
ptr = &atom->angmom[0][0];
nstride = 3;
} else if (thresh_array[ithresh] == ANGMOMY) {
if (!atom->angmom_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR,"Threshold for an atom property that isn't allocated");
ptr = &atom->angmom[0][1];
nstride = 3;
} else if (thresh_array[ithresh] == ANGMOMZ) {
if (!atom->angmom_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR,"Threshold for an atom property that isn't allocated");
ptr = &atom->angmom[0][2];
nstride = 3;
} else if (thresh_array[ithresh] == TQX) {
if (!atom->torque_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR,"Threshold for an atom property that isn't allocated");
ptr = &atom->torque[0][0];
nstride = 3;
} else if (thresh_array[ithresh] == TQY) {
if (!atom->torque_flag)
error->all(FLERR,
"Threshold for an atom property that isn't allocated");
error->all(FLERR,"Threshold for an atom property that isn't allocated");
ptr = &atom->torque[0][1];
nstride = 3;
} else if (thresh_array[ithresh] == TQZ) {
@ -928,7 +918,7 @@ void DumpVTK::write()
if (nmax > maxbuf) {
if ((bigint) nmax * size_one > MAXSMALLINT)
error->all(FLERR,"Too much per-proc info for dump");
error->all(FLERR,"Too much per processor data info for dump");
maxbuf = nmax;
memory->destroy(buf);
memory->create(buf,maxbuf*size_one,"dump:buf");
@ -983,9 +973,10 @@ void DumpVTK::write()
void DumpVTK::pack(tagint *ids)
{
int n = 0;
for (std::map<int,FnPtrPack>::iterator it=pack_choice.begin(); it!=pack_choice.end(); ++it, ++n) {
current_pack_choice_key = it->first; // work-around for pack_compute, pack_fix, pack_variable
(this->*(it->second))(n);
for (auto &choice : pack_choice) {
current_pack_choice_key = choice.first; // work-around for pack_compute, pack_fix, pack_variable
(this->*(choice.second))(n);
++n;
}
if (ids) {
@ -1108,10 +1099,10 @@ void DumpVTK::buf2arrays(int n, double *mybuf)
pid[0] = points->InsertNextPoint(mybuf[iatom*size_one],mybuf[iatom*size_one+1],mybuf[iatom*size_one+2]);
int j=3; // 0,1,2 = x,y,z handled just above
for (std::map<int, vtkSmartPointer<vtkAbstractArray> >::iterator it=myarrays.begin(); it!=myarrays.end(); ++it) {
vtkAbstractArray *paa = it->second;
if (it->second->GetNumberOfComponents() == 3) {
switch (vtype[it->first]) {
for (auto &it : myarrays) {
vtkAbstractArray *paa = it.second;
if (it.second->GetNumberOfComponents() == 3) {
switch (vtype[it.first]) {
case Dump::INT:
{
int iv3[3] = { static_cast<int>(mybuf[iatom*size_one+j ]),
@ -1130,7 +1121,7 @@ void DumpVTK::buf2arrays(int n, double *mybuf)
}
j+=3;
} else {
switch (vtype[it->first]) {
switch (vtype[it.first]) {
case Dump::INT:
{
vtkIntArray *pia = static_cast<vtkIntArray*>(paa);
@ -1317,8 +1308,8 @@ void DumpVTK::write_vtk(int n, double *mybuf)
unstructuredGrid->SetPoints(points);
unstructuredGrid->SetCells(VTK_VERTEX, pointsCells);
for (std::map<int, vtkSmartPointer<vtkAbstractArray> >::iterator it=myarrays.begin(); it!=myarrays.end(); ++it) {
unstructuredGrid->GetPointData()->AddArray(it->second);
for (const auto & it : myarrays) {
unstructuredGrid->GetPointData()->AddArray(it.second);
}
vtkSmartPointer<vtkUnstructuredGridWriter> writer = vtkSmartPointer<vtkUnstructuredGridWriter>::New();
@ -1327,8 +1318,8 @@ void DumpVTK::write_vtk(int n, double *mybuf)
polyData->SetPoints(points);
polyData->SetVerts(pointsCells);
for (std::map<int, vtkSmartPointer<vtkAbstractArray> >::iterator it=myarrays.begin(); it!=myarrays.end(); ++it) {
polyData->GetPointData()->AddArray(it->second);
for (auto &it : myarrays) {
polyData->GetPointData()->AddArray(it.second);
}
vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
@ -1384,8 +1375,8 @@ void DumpVTK::write_vtp(int n, double *mybuf)
polyData->SetPoints(points);
polyData->SetVerts(pointsCells);
for (std::map<int, vtkSmartPointer<vtkAbstractArray> >::iterator it=myarrays.begin(); it!=myarrays.end(); ++it) {
polyData->GetPointData()->AddArray(it->second);
for (auto &it : myarrays) {
polyData->GetPointData()->AddArray(it.second);
}
vtkSmartPointer<vtkXMLPolyDataWriter> writer = vtkSmartPointer<vtkXMLPolyDataWriter>::New();
@ -1448,8 +1439,8 @@ void DumpVTK::write_vtu(int n, double *mybuf)
unstructuredGrid->SetPoints(points);
unstructuredGrid->SetCells(VTK_VERTEX, pointsCells);
for (std::map<int, vtkSmartPointer<vtkAbstractArray> >::iterator it=myarrays.begin(); it!=myarrays.end(); ++it) {
unstructuredGrid->GetPointData()->AddArray(it->second);
for (auto &it : myarrays) {
unstructuredGrid->GetPointData()->AddArray(it.second);
}
vtkSmartPointer<vtkXMLUnstructuredGridWriter> writer = vtkSmartPointer<vtkXMLUnstructuredGridWriter>::New();
@ -1540,9 +1531,8 @@ int DumpVTK::parse_fields(int narg, char **arg)
name[Z] = "z";
// customize by adding to if statement
int i;
for (int iarg = 5; iarg < narg; iarg++) {
i = iarg-5;
for (int iarg = 0; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"id") == 0) {
pack_choice[ID] = &DumpVTK::pack_id;
@ -1550,7 +1540,7 @@ int DumpVTK::parse_fields(int narg, char **arg)
name[ID] = arg[iarg];
} else if (strcmp(arg[iarg],"mol") == 0) {
if (!atom->molecule_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property 'mol' that isn't allocated");
pack_choice[MOL] = &DumpVTK::pack_molecule;
vtype[MOL] = Dump::INT;
name[MOL] = arg[iarg];
@ -1665,98 +1655,97 @@ int DumpVTK::parse_fields(int narg, char **arg)
name[FZ] = arg[iarg];
} else if (strcmp(arg[iarg],"q") == 0) {
if (!atom->q_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[Q] = &DumpVTK::pack_q;
vtype[Q] = Dump::DOUBLE;
name[Q] = arg[iarg];
} else if (strcmp(arg[iarg],"mux") == 0) {
if (!atom->mu_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[MUX] = &DumpVTK::pack_mux;
vtype[MUX] = Dump::DOUBLE;
name[MUX] = arg[iarg];
} else if (strcmp(arg[iarg],"muy") == 0) {
if (!atom->mu_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[MUY] = &DumpVTK::pack_muy;
vtype[MUY] = Dump::DOUBLE;
name[MUY] = arg[iarg];
} else if (strcmp(arg[iarg],"muz") == 0) {
if (!atom->mu_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[MUZ] = &DumpVTK::pack_muz;
vtype[MUZ] = Dump::DOUBLE;
name[MUZ] = arg[iarg];
} else if (strcmp(arg[iarg],"mu") == 0) {
if (!atom->mu_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[MU] = &DumpVTK::pack_mu;
vtype[MU] = Dump::DOUBLE;
name[MU] = arg[iarg];
} else if (strcmp(arg[iarg],"radius") == 0) {
if (!atom->radius_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[RADIUS] = &DumpVTK::pack_radius;
vtype[RADIUS] = Dump::DOUBLE;
name[RADIUS] = arg[iarg];
} else if (strcmp(arg[iarg],"diameter") == 0) {
if (!atom->radius_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[DIAMETER] = &DumpVTK::pack_diameter;
vtype[DIAMETER] = Dump::DOUBLE;
name[DIAMETER] = arg[iarg];
} else if (strcmp(arg[iarg],"omegax") == 0) {
if (!atom->omega_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[OMEGAX] = &DumpVTK::pack_omegax;
vtype[OMEGAX] = Dump::DOUBLE;
name[OMEGAX] = arg[iarg];
} else if (strcmp(arg[iarg],"omegay") == 0) {
if (!atom->omega_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[OMEGAY] = &DumpVTK::pack_omegay;
vtype[OMEGAY] = Dump::DOUBLE;
name[OMEGAY] = arg[iarg];
} else if (strcmp(arg[iarg],"omegaz") == 0) {
if (!atom->omega_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[OMEGAZ] = &DumpVTK::pack_omegaz;
vtype[OMEGAZ] = Dump::DOUBLE;
name[OMEGAZ] = arg[iarg];
} else if (strcmp(arg[iarg],"angmomx") == 0) {
if (!atom->angmom_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[ANGMOMX] = &DumpVTK::pack_angmomx;
vtype[ANGMOMX] = Dump::DOUBLE;
name[ANGMOMX] = arg[iarg];
} else if (strcmp(arg[iarg],"angmomy") == 0) {
if (!atom->angmom_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[ANGMOMY] = &DumpVTK::pack_angmomy;
vtype[ANGMOMY] = Dump::DOUBLE;
name[ANGMOMY] = arg[iarg];
} else if (strcmp(arg[iarg],"angmomz") == 0) {
if (!atom->angmom_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[ANGMOMZ] = &DumpVTK::pack_angmomz;
vtype[ANGMOMZ] = Dump::DOUBLE;
name[ANGMOMZ] = arg[iarg];
} else if (strcmp(arg[iarg],"tqx") == 0) {
if (!atom->torque_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[TQX] = &DumpVTK::pack_tqx;
vtype[TQX] = Dump::DOUBLE;
name[TQX] = arg[iarg];
} else if (strcmp(arg[iarg],"tqy") == 0) {
if (!atom->torque_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[TQY] = &DumpVTK::pack_tqy;
vtype[TQY] = Dump::DOUBLE;
name[TQY] = arg[iarg];
} else if (strcmp(arg[iarg],"tqz") == 0) {
if (!atom->torque_flag)
error->all(FLERR,"Dumping an atom property that isn't allocated");
error->all(FLERR,"Dumping atom property '{}' that isn't allocated", arg[iarg]);
pack_choice[TQZ] = &DumpVTK::pack_tqz;
vtype[TQZ] = Dump::DOUBLE;
name[TQZ] = arg[iarg];
@ -1766,7 +1755,7 @@ int DumpVTK::parse_fields(int narg, char **arg)
int n,flag,cols;
ArgInfo argi(arg[iarg],ArgInfo::COMPUTE|ArgInfo::FIX|ArgInfo::VARIABLE
|ArgInfo::DNAME|ArgInfo::INAME);
argindex[ATTRIBUTES+i] = argi.get_index1();
argindex[ATTRIBUTES+iarg] = argi.get_index1();
auto aname = argi.get_name();
switch (argi.get_type()) {
@ -1779,107 +1768,114 @@ int DumpVTK::parse_fields(int narg, char **arg)
// if no trailing [], then arg is set to 0, else arg is int between []
case ArgInfo::COMPUTE:
pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_compute;
vtype[ATTRIBUTES+i] = Dump::DOUBLE;
{
pack_choice[ATTRIBUTES+iarg] = &DumpVTK::pack_compute;
vtype[ATTRIBUTES+iarg] = Dump::DOUBLE;
n = modify->find_compute(aname);
if (n < 0) error->all(FLERR,"Could not find dump vtk compute ID: {}",aname);
if (modify->compute[n]->peratom_flag == 0)
error->all(FLERR,"Dump vtk compute {} does not compute per-atom info",aname);
if (argi.get_dim() == 0 && modify->compute[n]->size_peratom_cols > 0)
error->all(FLERR,"Dump vtk compute {} does not calculate per-atom vector",aname);
if (argi.get_dim() > 0 && modify->compute[n]->size_peratom_cols == 0)
error->all(FLERR,"Dump vtk compute {} does not calculate per-atom array",aname);
if (argi.get_dim() > 0 &&
argi.get_index1() > modify->compute[n]->size_peratom_cols)
error->all(FLERR,"Dump vtk compute {} vector is accessed out-of-range",aname);
field2index[ATTRIBUTES+i] = add_compute(aname);
name[ATTRIBUTES+i] = arg[iarg];
auto icompute = modify->get_compute_by_id(aname);
if (!icompute) {
error->all(FLERR,"Could not find dump vtk compute ID: {}",aname);
} else {
if (icompute->peratom_flag == 0)
error->all(FLERR,"Dump vtk compute {} does not compute per-atom info",aname);
if (argi.get_dim() == 0 && icompute->size_peratom_cols > 0)
error->all(FLERR,"Dump vtk compute {} does not calculate per-atom vector",aname);
if (argi.get_dim() > 0 && icompute->size_peratom_cols == 0)
error->all(FLERR,"Dump vtk compute {} does not calculate per-atom array",aname);
if (argi.get_dim() > 0 && argi.get_index1() > icompute->size_peratom_cols)
error->all(FLERR,"Dump vtk compute {} vector is accessed out-of-range",aname);
field2index[ATTRIBUTES+iarg] = add_compute(aname);
name[ATTRIBUTES+iarg] = arg[iarg];
}
break;
}
// fix value = f_ID
// if no trailing [], then arg is set to 0, else arg is between []
case ArgInfo::FIX:
pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_fix;
vtype[ATTRIBUTES+i] = Dump::DOUBLE;
{
pack_choice[ATTRIBUTES+iarg] = &DumpVTK::pack_fix;
vtype[ATTRIBUTES+iarg] = Dump::DOUBLE;
n = modify->find_fix(aname);
if (n < 0) error->all(FLERR,"Could not find dump vtk fix ID: {}",aname);
if (modify->fix[n]->peratom_flag == 0)
error->all(FLERR,"Dump vtk fix {} does not compute per-atom info",aname);
if (argi.get_dim() == 0 && modify->fix[n]->size_peratom_cols > 0)
error->all(FLERR,"Dump vtk fix {} does not compute per-atom vector",aname);
if (argi.get_dim() > 0 && modify->fix[n]->size_peratom_cols == 0)
error->all(FLERR,"Dump vtk fix {} does not compute per-atom array",aname);
if (argi.get_dim() > 0 &&
argi.get_index1() > modify->fix[n]->size_peratom_cols)
error->all(FLERR,"Dump vtk fix {} vector is accessed out-of-range",aname);
auto ifix = modify->get_fix_by_id(aname);
if (!ifix) {
error->all(FLERR,"Could not find dump vtk fix ID: {}",aname);
} else {
if (ifix->peratom_flag == 0)
error->all(FLERR,"Dump vtk fix {} does not compute per-atom info",aname);
if (argi.get_dim() == 0 && ifix->size_peratom_cols > 0)
error->all(FLERR,"Dump vtk fix {} does not compute per-atom vector",aname);
if (argi.get_dim() > 0 && ifix->size_peratom_cols == 0)
error->all(FLERR,"Dump vtk fix {} does not compute per-atom array",aname);
if (argi.get_dim() > 0 && argi.get_index1() > ifix->size_peratom_cols)
error->all(FLERR,"Dump vtk fix {} vector is accessed out-of-range",aname);
field2index[ATTRIBUTES+i] = add_fix(aname);
name[ATTRIBUTES+i] = arg[iarg];
field2index[ATTRIBUTES+iarg] = add_fix(aname);
name[ATTRIBUTES+iarg] = arg[iarg];
}
break;
}
// variable value = v_name
case ArgInfo::VARIABLE:
pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_variable;
vtype[ATTRIBUTES+i] = Dump::DOUBLE;
pack_choice[ATTRIBUTES+iarg] = &DumpVTK::pack_variable;
vtype[ATTRIBUTES+iarg] = Dump::DOUBLE;
n = input->variable->find(aname);
if (n < 0) error->all(FLERR,"Could not find dump vtk variable name {}",aname);
if (input->variable->atomstyle(n) == 0)
error->all(FLERR,"Dump vtk variable {} is not atom-style variable",aname);
field2index[ATTRIBUTES+i] = add_variable(aname);
name[ATTRIBUTES+i] = arg[iarg];
field2index[ATTRIBUTES+iarg] = add_variable(aname);
name[ATTRIBUTES+iarg] = arg[iarg];
break;
// custom per-atom floating point vector or array = d_ID d2_ID
case ArgInfo::DNAME:
pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_custom;
vtype[ATTRIBUTES+i] = Dump::DOUBLE;
pack_choice[ATTRIBUTES+iarg] = &DumpVTK::pack_custom;
vtype[ATTRIBUTES+iarg] = Dump::DOUBLE;
n = atom->find_custom(aname,flag,cols);
if (n < 0)
error->all(FLERR,"Could not find custom per-atom property ID: {}", aname);
if (argindex[ATTRIBUTES+i] == 0) {
if (argindex[ATTRIBUTES+iarg] == 0) {
if (!flag || cols)
error->all(FLERR,"Property double vector {} for dump vtk does not exist",aname);
} else {
if (!flag || !cols)
error->all(FLERR,"Property double array {} for dump vtk does not exist",aname);
if (argindex[ATTRIBUTES+i] > atom->dcols[n])
if (argindex[ATTRIBUTES+iarg] > atom->dcols[n])
error->all(FLERR,"Dump vtk property array {} is accessed out-of-range",aname);
}
field2index[ATTRIBUTES+i] = add_custom(aname,1);
name[ATTRIBUTES+i] = arg[iarg];
field2index[ATTRIBUTES+iarg] = add_custom(aname,1);
name[ATTRIBUTES+iarg] = arg[iarg];
break;
// custom per-atom integer vector or array = i_ID or i2_ID
case ArgInfo::INAME:
pack_choice[ATTRIBUTES+i] = &DumpVTK::pack_custom;
vtype[ATTRIBUTES+i] = Dump::INT;
pack_choice[ATTRIBUTES+iarg] = &DumpVTK::pack_custom;
vtype[ATTRIBUTES+iarg] = Dump::INT;
n = atom->find_custom(aname,flag,cols);
if (n < 0)
error->all(FLERR,"Could not find custom per-atom property ID: {}", aname);
if (argindex[ATTRIBUTES+i] == 0) {
if (argindex[ATTRIBUTES+iarg] == 0) {
if (flag || cols)
error->all(FLERR,"Property integer vector {} for dump vtk does not exist",aname);
} else {
if (flag || !cols)
error->all(FLERR,"Property integer array {} for dump vtk does not exist",aname);
if (argindex[ATTRIBUTES+i] > atom->icols[n])
if (argindex[ATTRIBUTES+iarg] > atom->icols[n])
error->all(FLERR,"Dump vtk property array {} is accessed out-of-range",aname);
}
field2index[ATTRIBUTES+i] = add_custom(aname,0);
name[ATTRIBUTES+i] = arg[iarg];
field2index[ATTRIBUTES+iarg] = add_custom(aname,0);
name[ATTRIBUTES+iarg] = arg[iarg];
break;
// no match
@ -2232,42 +2228,50 @@ int DumpVTK::modify_param(int narg, char **arg)
// if no trailing [], then arg is set to 0, else arg is between []
case ArgInfo::COMPUTE:
{
thresh_array[nthresh] = COMPUTE;
n = modify->find_compute(aname);
if (n < 0) error->all(FLERR,"Could not find dump modify compute ID: {}",aname);
auto icompute = modify->get_compute_by_id(aname);
if (!icompute) {
error->all(FLERR,"Could not find dump modify compute ID: {}",aname);
} else {
if (modify->compute[n]->peratom_flag == 0)
error->all(FLERR,"Dump modify compute ID {} does not compute per-atom info",aname);
if (argi.get_dim() == 0 && icompute->size_peratom_cols > 0)
error->all(FLERR,"Dump modify compute ID {} does not compute per-atom vector",aname);
if (argi.get_index1() > 0 && icompute->size_peratom_cols == 0)
error->all(FLERR,"Dump modify compute ID {} does not compute per-atom array",aname);
if (argi.get_index1() > 0 &&
argi.get_index1() > icompute->size_peratom_cols)
error->all(FLERR,"Dump modify compute ID {} vector is not large enough",aname);
if (modify->compute[n]->peratom_flag == 0)
error->all(FLERR,"Dump modify compute ID {} does not compute per-atom info",aname);
if (argi.get_dim() == 0 && modify->compute[n]->size_peratom_cols > 0)
error->all(FLERR,"Dump modify compute ID {} does not compute per-atom vector",aname);
if (argi.get_index1() > 0 && modify->compute[n]->size_peratom_cols == 0)
error->all(FLERR,"Dump modify compute ID {} does not compute per-atom array",aname);
if (argi.get_index1() > 0 &&
argi.get_index1() > modify->compute[n]->size_peratom_cols)
error->all(FLERR,"Dump modify compute ID {} vector is not large enough",aname);
field2index[ATTRIBUTES+nfield+nthresh] = add_compute(aname);
field2index[ATTRIBUTES+nfield+nthresh] = add_compute(aname);
}
break;
}
// fix value = f_ID
// if no trailing [], then arg is set to 0, else arg is between []
case ArgInfo::FIX:
{
thresh_array[nthresh] = FIX;
n = modify->find_fix(aname);
if (n < 0) error->all(FLERR,"Could not find dump modify fix ID: {}",aname);
auto ifix = modify->get_fix_by_id(aname);
if (!ifix) {
error->all(FLERR,"Could not find dump modify fix ID: {}",aname);
} else {
if (ifix->peratom_flag == 0)
error->all(FLERR,"Dump modify fix ID {} does not compute per-atom info",aname);
if (argi.get_dim() == 0 && ifix->size_peratom_cols > 0)
error->all(FLERR,"Dump modify fix ID {} does not compute per-atom vector",aname);
if (argi.get_index1() > 0 && ifix->size_peratom_cols == 0)
error->all(FLERR,"Dump modify fix ID {} does not compute per-atom array",aname);
if (argi.get_index1() > 0 && argi.get_index1() > ifix->size_peratom_cols)
error->all(FLERR,"Dump modify fix ID {} vector is not large enough",aname);
if (modify->fix[n]->peratom_flag == 0)
error->all(FLERR,"Dump modify fix ID {} does not compute per-atom info",aname);
if (argi.get_dim() == 0 && modify->fix[n]->size_peratom_cols > 0)
error->all(FLERR,"Dump modify fix ID {} does not compute per-atom vector",aname);
if (argi.get_index1() > 0 && modify->fix[n]->size_peratom_cols == 0)
error->all(FLERR,"Dump modify fix ID {} does not compute per-atom array",aname);
if (argi.get_index1() > 0 && argi.get_index1() > modify->fix[n]->size_peratom_cols)
error->all(FLERR,"Dump modify fix ID {} vector is not large enough",aname);
field2index[ATTRIBUTES+nfield+nthresh] = add_fix(aname);
field2index[ATTRIBUTES+nfield+nthresh] = add_fix(aname);
}
break;
}
// variable value = v_ID

148
tools/lammps-gui/update-wham-2.0.11.patch
View File

@ -55,30 +55,31 @@ index 0000000..b4f0fe6
+ TYPE DOC
+ PERMISSIONS OWNER_READ GROUP_READ WORLD_READ
+)
diff --git a/wham/wham.c b/wham/wham.c
index 487871b..526908e 100644
--- a/wham/wham.c
+++ b/wham/wham.c
@@ -21,7 +21,7 @@
#include "wham.h"
diff --git a/wham-2d/wham-2d.c b/wham-2d/wham-2d.c
index fb6e059..2c5594f 100644
--- a/wham-2d/wham-2d.c
+++ b/wham-2d/wham-2d.c
@@ -25,7 +25,7 @@
#include <time.h>
#include "wham-2d.h"
-#define COMMAND_LINE "Command line: wham [P|Ppi|Pval] hist_min hist_max num_bins tol temperature numpad metadatafile freefile [num_MC_trials randSeed]\n"
+#define COMMAND_LINE "Command line: wham [P|Ppi|Pval] [units <real|metal|lj|...>] hist_min hist_max num_bins tol temperature numpad metadatafile freefile [num_MC_trials randSeed]\n"
double HIST_MAX,HIST_MIN,BIN_WIDTH,TOL;
double *HISTOGRAM;
@@ -29,6 +29,7 @@ double kT;
int NUM_BINS;
int PERIODIC;
double PERIOD;
-#define COMMAND_LINE "Command line: wham-2d Px[=0|pi|val] hist_min_x hist_max_x num_bins_x Py[=0|pi|val] hist_min_y hist_max_y num_bins_y tol temperature numpad metadatafile freefile use_mask\n"
+#define COMMAND_LINE "Command line: wham-2d [units <real|metal|lj|...>] Px[=0|pi|val] hist_min_x hist_max_x num_bins_x Py[=0|pi|val] hist_min_y hist_max_y num_bins_y tol temperature numpad metadatafile freefile use_mask\n"
double HIST_MAXx,HIST_MINx,BIN_WIDTHx;
double HIST_MAXy,HIST_MINy,BIN_WIDTHy;
double TOL;
@@ -35,7 +35,7 @@ int NUM_BINSx, NUM_BINSy;
int PERIODICx, PERIODICy;
double PERIODx, PERIODy;
double *data1,**num,***bias;
-
+double k_B = k_B_DEFAULT;
int main(int argc, char *argv[])
{
@@ -117,6 +118,61 @@ else
PERIODIC = 0;
@@ -76,6 +76,61 @@ for (i=0; i<argc; i++)
}
printf("\n");
+// set k_B according to LAMMPS units settings
+if (strcmp(argv[1],"units") == 0)
@ -135,9 +136,116 @@ index 487871b..526908e 100644
+ argv += 2;
+ }
+
// Parse command line arguments
if (argc != 9 && argc !=11)
PERIODICx = parse_periodic(argv[1], &PERIODx);
if (PERIODICx)
{
diff --git a/wham-2d/wham-2d.h b/wham-2d/wham-2d.h
index b17e4bd..5fc17ff 100644
--- a/wham-2d/wham-2d.h
+++ b/wham-2d/wham-2d.h
@@ -20,15 +20,15 @@ extern int NUM_BINSy;
extern int PERIODICx,PERIODICy; // flags to turn on periodicity
extern double PERIODx, PERIODy; // flags to control periodic interval
+extern double k_B;
// A couple of predefined periodic units
#define DEGREES 360.0
#define RADIANS 6.28318530717959
-#define k_B 0.001982923700 // Boltzmann's constant in kcal/mol K
-//#define k_B 0.0083144621 // Boltzmann's constant kJ/mol-K
-//#define k_B 1.0 // Boltzmann's constant in reduced units
-
+#define k_B_DEFAULT 0.001982923700 // Boltzmann's constant in kcal/mol K
+//#define k_B_DEFAULT 0.0083144621 // Boltzmann's constant kJ/mol-K
+//#define k_B_DEFAULT 1.0 // Boltzmann's constant in reduced units
// Value inserted for the free energy of masked values
diff --git a/wham/wham.c b/wham/wham.c
index 487871b..1496eed 100644
--- a/wham/wham.c
+++ b/wham/wham.c
@@ -21,7 +21,7 @@
#include "wham.h"
-#define COMMAND_LINE "Command line: wham [P|Ppi|Pval] hist_min hist_max num_bins tol temperature numpad metadatafile freefile [num_MC_trials randSeed]\n"
+#define COMMAND_LINE "Command line: wham [units <real|metal|lj|...>] [P|Ppi|Pval] hist_min hist_max num_bins tol temperature numpad metadatafile freefile [num_MC_trials randSeed]\n"
double HIST_MAX,HIST_MIN,BIN_WIDTH,TOL;
double *HISTOGRAM;
@@ -29,6 +29,7 @@ double kT;
int NUM_BINS;
int PERIODIC;
double PERIOD;
+double k_B = k_B_DEFAULT;
int main(int argc, char *argv[])
{
@@ -82,6 +83,61 @@ for (i=0; i<argc; i++)
}
printf("\n");
+// set k_B according to LAMMPS units settings
+if (strcmp(argv[1],"units") == 0)
+ {
+ if (argc < 3)
+ {
+ printf( COMMAND_LINE );
+ exit(-1);
+ }
+
+ if (strcmp(argv[2], "lj") == 0)
+ {
+ k_B = 1.0;
+ }
+ else if (strcmp(argv[2], "real") == 0)
+ {
+ k_B = 0.0019872067;
+ }
+ else if (strcmp(argv[2], "metal") == 0)
+ {
+ k_B = 8.617343e-5;
+ }
+ else if (strcmp(argv[2], "si") == 0)
+ {
+ k_B = 1.3806504e-23;
+ }
+ else if (strcmp(argv[2], "cgs") == 0)
+ {
+ k_B = 1.3806504e-16;
+ }
+ else if (strcmp(argv[2], "electron") == 0)
+ {
+ k_B = 3.16681534e-6;
+ }
+ else if (strcmp(argv[2], "micro") == 0)
+ {
+ k_B = 1.3806504e-8;
+ }
+ else if (strcmp(argv[2], "nano") == 0)
+ {
+ k_B = 0.013806504;
+ }
+ else if (strcmp(argv[2], "default") == 0)
+ {
+ k_B = k_B_DEFAULT;
+ }
+ else
+ {
+ printf("Unknown unit style: %s\n%s", argv[2], COMMAND_LINE);
+ exit(-1);
+ }
+ printf("# Setting value of k_B to = %.15g\n", k_B);
+ argc -= 2;
+ argv += 2;
+ }
+
if (toupper(argv[1][0]) == 'P')
{
PERIODIC = 1;
diff --git a/wham/wham.h b/wham/wham.h
index aacc1e8..7d509f2 100644
--- a/wham/wham.h