diff --git a/doc/src/Commands_all.rst b/doc/src/Commands_all.rst index bd9d0ed09c..2ba5cba51d 100644 --- a/doc/src/Commands_all.rst +++ b/doc/src/Commands_all.rst @@ -24,6 +24,7 @@ table above. * :doc:`angle_coeff ` * :doc:`angle_style ` + * :doc:`angle_write ` * :doc:`atom_modify ` * :doc:`atom_style ` * :doc:`balance ` diff --git a/doc/src/angle_table.rst b/doc/src/angle_table.rst index c0faa3f046..fba41db045 100644 --- a/doc/src/angle_table.rst +++ b/doc/src/angle_table.rst @@ -59,9 +59,12 @@ format of this file is described below. ---------- -Suitable tables for use with this angle style can be created using the -Python code in the ``tools/tabulate`` folder of the LAMMPS source code -distribution. +Suitable tables for use with this angle style can be created by LAMMPS +itself from existing angle styles using the :doc:`angle_write +` command. This can be useful to have a template file for +testing the angle style settings and to build a compatible custom file. +Another option to generate tables is the Python code in the +``tools/tabulate`` folder of the LAMMPS source code distribution. The format of a tabulated file is as follows (without the parenthesized comments): @@ -154,7 +157,7 @@ for more info. Related commands """""""""""""""" -:doc:`angle_coeff ` +:doc:`angle_coeff `, :doc:`angle_write ` Default """"""" diff --git a/doc/src/angle_write.rst b/doc/src/angle_write.rst new file mode 100644 index 0000000000..16a4ee4e2c --- /dev/null +++ b/doc/src/angle_write.rst @@ -0,0 +1,98 @@ +.. index:: angle_write + +angle_write command +=================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + angle_write atype N file keyword + +* atype = angle type +* N = # of values +* file = name of file to write values to +* keyword = section name in file for this set of tabulated values + +Examples +"""""""" + +.. code-block:: LAMMPS + + angle_write 1 500 table.txt Harmonic_1 + angle_write 3 1000 table.txt Harmonic_3 + +Description +""""""""""" + +Write energy and force values to a file as a function of angle for the +currently defined angle potential. Force in this context means the +angle force, not the force on individual atoms. This is useful for +plotting the potential function or otherwise debugging its values. The +resulting file can also be used as input for use with :doc:`angle style +table `. + +If the file already exists, the table of values is appended to the end +of the file to allow multiple tables of energy and force to be included +in one file. The individual sections may be identified by the *keyword*. + +The energy and force values are computed for angles ranging from 0 +degrees to 180 degrees for 3 interacting atoms forming an angle type +atype, using the appropriate :doc:`angle_coeff ` +coefficients. N evenly spaced angles are used. + +For example, for N = 6, values are computed at :math:`\theta = 0, 36, +72, 108, 144, 180`. + +The file is written in the format used as input for the +:doc:`angle_style table ` option with *keyword* as the +section name. Each line written to the file lists an index number +(1-N), an angle (in degrees), an energy (in energy units), and a force +(in force units per radians^2). In case a new file is created, the +first line will be a comment with a "DATE:" and "UNITS:" tag with the +current date and :doc:`units ` settings. For subsequent +invocations of the *angle_write* command for the same file, data will be +appended and the current units settings will be compared to the data +from the header, if present. The *angle_write* will refuse to add a +table to an existing file if the units are not the same. + +Restrictions +"""""""""""" + +All force field coefficients for angle and other kinds of interactions +must be set before this command can be invoked. + +The table of the angle energy and force data data is created by using a +separate, internally created, new LAMMPS instance with a dummy system of +3 atoms for which the angle potential energy is computed after +transferring the angle style and coefficients and arranging the 3 atoms +into the corresponding geometries. The angle force is then determined +from the potential energies through numerical differentiation. As a +consequence or this approach, not all angle styles are compatible. The +following conditions must be met: + +- The angle style must be able to write its coefficients to a data file. + This condition excludes for example :doc:`angle style hybrid ` and + :doc:`angle style table `. +- The potential function must not have any terms that depend on geometry + properties other than the angle. This condition excludes for example + all angle types for :doc:`angle style charmm ` that have + non-zero Urey-Bradley terms. Please note that the *write_angle* + command has no way of checking for this condition, so the resulting tables + may be bogus. It is strongly recommended to make careful tests for any + created tables. +- There must not be multiple angle interactions defined for the same + triple of atoms. This applies for example to :doc:`angle_style class2 + `. + +Related commands +"""""""""""""""" + +:doc:`angle_style table `, :doc:`bond_write `, +:doc:`angle_style `, :doc:`angle_coeff ` + +Default +""""""" + +none diff --git a/doc/src/bond_write.rst b/doc/src/bond_write.rst index 4170029239..51e0653f51 100644 --- a/doc/src/bond_write.rst +++ b/doc/src/bond_write.rst @@ -67,7 +67,7 @@ be specified even if the potential has a finite value at r = 0.0. Related commands """""""""""""""" -:doc:`bond_style table `, +:doc:`bond_style table `, `angle_write `, :doc:`bond_style `, :doc:`bond_coeff ` Default diff --git a/doc/src/commands_list.rst b/doc/src/commands_list.rst index ff2502189d..5731cbb755 100644 --- a/doc/src/commands_list.rst +++ b/doc/src/commands_list.rst @@ -6,6 +6,7 @@ Commands angle_coeff angle_style + angle_write atom_modify atom_style balance diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 7cb409d040..13cc96f993 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -174,6 +174,7 @@ attrac Atw Atwater atwt +atype augt AuO automagically diff --git a/src/angle_write.cpp b/src/angle_write.cpp new file mode 100644 index 0000000000..e0d6827a66 --- /dev/null +++ b/src/angle_write.cpp @@ -0,0 +1,244 @@ +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing authors: Axel Kohlmeyer (Temple U) +------------------------------------------------------------------------- */ + +#include "angle_write.h" + +#include "angle.h" +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "domain.h" +#include "error.h" +#include "force.h" +#include "input.h" +#include "lammps.h" +#include "math_const.h" +#include "update.h" + +#include +using namespace LAMMPS_NS; +using MathConst::DEG2RAD; +using MathConst::RAD2DEG; + +static constexpr double epsilon = 6.5e-6; +#define MAXLINE 1024 +/* ---------------------------------------------------------------------- */ + +void AngleWrite::command(int narg, char **arg) +{ + // sanity checks + + if (domain->box_exist == 0) + error->all(FLERR, "Angle_write command before simulation box is defined"); + if (atom->avec->angles_allow == 0) + error->all(FLERR, "Angle_write command when no angles allowed"); + auto angle = force->angle; + if (angle == nullptr) error->all(FLERR, "Angle_write command before an angle_style is defined"); + if (angle && (force->angle->writedata == 0)) + error->all(FLERR, "Angle style must support writing coeffs to data file for angle_write"); + + // parse arguments + + if (narg != 4) error->all(FLERR, "Illegal angle_write command"); + + int atype = utils::inumeric(FLERR, arg[0], false, lmp); + int n = utils::inumeric(FLERR, arg[1], false, lmp); + std::string table_file = arg[2]; + std::string keyword = arg[3]; + + // make sure system is initialized before calling any functions + + lmp->init(); + + double theta0 = angle->equilibrium_angle(atype) * RAD2DEG; + + // write out all angle_coeff settings to file. use function from write_data. + // open table file in append mode if it already exists + // add line with DATE: and UNITS: tag when creating new file + // otherwise make certain that units are consistent + // print header in format used by angle_style table + + FILE *fp = nullptr; + std::string coeffs_file = table_file + ".tmp.coeffs"; + if (comm->me == 0) { + + fp = fopen(coeffs_file.c_str(), "w"); + force->angle->write_data(fp); + fclose(fp); + + // units sanity check: + // - if this is the first time we write to this potential file, + // write out a line with "DATE:" and "UNITS:" tags + // - if the file already exists, print a message about appending + // while printing the date and check that units are consistent. + if (platform::file_is_readable(table_file)) { + std::string units = utils::get_potential_units(table_file, "table"); + if (!units.empty() && (units != update->unit_style)) { + error->one(FLERR, "Trying to append to a table file with UNITS: {} while units are {}", + units, update->unit_style); + } + std::string date = utils::get_potential_date(table_file, "table"); + utils::logmesg(lmp, "Appending to table file {} with DATE: {}\n", table_file, date); + fp = fopen(table_file.c_str(), "a"); + } else { + utils::logmesg(lmp, "Creating table file {} with DATE: {}\n", table_file, + utils::current_date()); + fp = fopen(table_file.c_str(), "w"); + if (fp) + fmt::print(fp, "# DATE: {} UNITS: {} Created by angle_write\n", utils::current_date(), + update->unit_style); + } + if (fp == nullptr) + error->one(FLERR, "Cannot open angle_write file {}: {}", table_file, utils::getsyserror()); + } + + // split communicator so that we can run a new LAMMPS class instance only on comm->me == 0 + + MPI_Comm singlecomm; + int color = (comm->me == 0) ? 1 : MPI_UNDEFINED; + int key = comm->me; + MPI_Comm_split(world, color, key, &singlecomm); + + if (comm->me == 0) { + // set up new LAMMPS instance with dummy system to evaluate angle potential + const char *args[] = {"AngleWrite", "-nocite", "-echo", "none", + "-log", "none", "-screen", "none"}; + char **argv = (char **) args; + int argc = sizeof(args) / sizeof(char *); + LAMMPS *writer = new LAMMPS(argc, argv, singlecomm); + + // create dummy system replicating angle style settings + writer->input->one(fmt::format("units {}", update->unit_style)); + writer->input->one("atom_style angle"); + writer->input->one("atom_modify map array"); + writer->input->one("boundary f f f"); + writer->input->one("region box block -2 2 -2 2 -1 1"); + writer->input->one(fmt::format("create_box {} box angle/types {} " + "extra/angle/per/atom 1 " + "extra/special/per/atom 4", + atom->ntypes, atom->nangletypes)); + writer->input->one("create_atoms 1 single 0.0 0.0 0.0"); + writer->input->one("create_atoms 1 single 1.0 0.0 0.0"); + writer->input->one("create_atoms 1 single -1.0 0.0 0.0"); + writer->input->one(fmt::format("create_bonds single/angle {} 2 1 3", atype)); + + writer->input->one("pair_style zero 10.0"); + writer->input->one("pair_coeff * *"); + writer->input->one("mass * 1.0"); + writer->input->one(fmt::format("angle_style {}", force->angle_style)); + FILE *coeffs; + char line[MAXLINE]; + coeffs = fopen(coeffs_file.c_str(), "r"); + for (int i = 0; i < atom->nangletypes; ++i) { + fgets(line, MAXLINE, coeffs); + writer->input->one(fmt::format("angle_coeff {}", line)); + } + platform::unlink(coeffs_file); + + // initialize system + + writer->init(); + + // move third atom to reproduce angles + + double theta, phi, phi1, phi2, f; + angle = writer->force->angle; + int i1, i2, i3; + i1 = writer->atom->map(1); + i2 = writer->atom->map(2); + i3 = writer->atom->map(3); + auto atom3 = writer->atom->x[i3]; + + // evaluate energy and force at each of N distances + + fmt::print(fp, "# Angle potential {} for angle type {}: i,theta,energy,force\n", + force->angle_style, atype); + fmt::print(fp, "\n{}\nN {} EQ {:.15g}\n\n", keyword, n, theta0); + + const double dtheta = 180.0 / static_cast(n - 1); + + // get force for divergent 0 degree angle from interpolation to the right + theta = 0.0; + atom3[0] = cos(theta * DEG2RAD); + atom3[1] = sin(theta * DEG2RAD); + phi = angle->single(atype, i2, i1, i3); + + theta = epsilon; + atom3[0] = cos(theta * DEG2RAD); + atom3[1] = sin(theta * DEG2RAD); + phi1 = angle->single(atype, i2, i1, i3); + + theta = 2.0 * epsilon; + atom3[0] = cos(theta * DEG2RAD); + atom3[1] = sin(theta * DEG2RAD); + phi2 = angle->single(atype, i2, i1, i3); + + f = (1.5 * phi - 2.0 * phi1 + 0.5 * phi2) / epsilon; + if (!std::isfinite(f)) f = 0.0; + if (!std::isfinite(phi)) phi = 0.0; + fprintf(fp, "%8d %- 22.15g %- 22.15g %- 22.15g\n", 1, 0.0, phi, f); + + for (int i = 1; i < n - 1; i++) { + theta = dtheta * static_cast(i) - epsilon; + atom3[0] = cos(theta * DEG2RAD); + atom3[1] = sin(theta * DEG2RAD); + phi1 = angle->single(atype, i2, i1, i3); + + theta = dtheta * static_cast(i) + epsilon; + atom3[0] = cos(theta * DEG2RAD); + atom3[1] = sin(theta * DEG2RAD); + phi2 = angle->single(atype, i2, i1, i3); + + theta = dtheta * static_cast(i); + atom3[0] = cos(theta * DEG2RAD); + atom3[1] = sin(theta * DEG2RAD); + phi = angle->single(atype, i2, i1, i3); + if (!std::isfinite(phi)) phi = 0.0; + + // get force from numerical differentiation + f = -0.5 * (phi2 - phi1) / epsilon; + if (!std::isfinite(f)) f = 0.0; + fprintf(fp, "%8d %- 22.15g %- 22.15g %- 22.15g\n", i + 1, theta, phi, f); + } + + // get force for divergent 180 degree angle from interpolation to the left + theta = 180.0; + atom3[0] = cos(theta * DEG2RAD); + atom3[1] = sin(theta * DEG2RAD); + phi = angle->single(atype, i2, i1, i3); + + theta = 180.0 - epsilon; + atom3[0] = cos(theta * DEG2RAD); + atom3[1] = sin(theta * DEG2RAD); + phi1 = angle->single(atype, i2, i1, i3); + + theta = 180.0 - 2.0 * epsilon; + atom3[0] = cos(theta * DEG2RAD); + atom3[1] = sin(theta * DEG2RAD); + phi2 = angle->single(atype, i2, i1, i3); + + f = (2.0 * phi1 - 1.5 * phi - 0.5 * phi2) / epsilon; + if (!std::isfinite(f)) f = 0.0; + if (!std::isfinite(phi)) phi = 0.0; + fprintf(fp, "%8d %- 22.15g %- 22.15g %- 22.15g\n", 1, 180.0, phi, f); + + // clean up + delete writer; + fclose(fp); + } + MPI_Comm_free(&singlecomm); +} diff --git a/src/angle_write.h b/src/angle_write.h new file mode 100644 index 0000000000..5c56cd8efc --- /dev/null +++ b/src/angle_write.h @@ -0,0 +1,34 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + Copyright (2003) Sandia Corporation. Under the terms of Contract + DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains + certain rights in this software. This software is distributed under + the GNU General Public License. + + See the README file in the top-level LAMMPS directory. +------------------------------------------------------------------------- */ + +#ifdef COMMAND_CLASS +// clang-format off +CommandStyle(angle_write,AngleWrite); +// clang-format on +#else + +#ifndef LMP_ANGLE_WRITE_H +#define LMP_ANGLE_WRITE_H + +#include "command.h" + +namespace LAMMPS_NS { + +class AngleWrite : public Command { + public: + AngleWrite(class LAMMPS *lmp) : Command(lmp){}; + void command(int, char **) override; +}; +} // namespace LAMMPS_NS +#endif +#endif