From da98363a25ec6cb6f18aed85fc621525fbf94fa7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 24 Dec 2022 22:49:42 -0500 Subject: [PATCH 1/8] implement angle_write command --- doc/src/Commands_all.rst | 1 + doc/src/angle_table.rst | 11 +- doc/src/angle_write.rst | 98 ++++++++ doc/src/bond_write.rst | 2 +- doc/src/commands_list.rst | 1 + doc/utils/sphinx-config/false_positives.txt | 1 + src/angle_write.cpp | 244 ++++++++++++++++++++ src/angle_write.h | 34 +++ 8 files changed, 387 insertions(+), 5 deletions(-) create mode 100644 doc/src/angle_write.rst create mode 100644 src/angle_write.cpp create mode 100644 src/angle_write.h 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 From bbdc6fd3ab1ef7ce093fcc6337342e680d0942d5 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 25 Dec 2022 04:35:03 -0500 Subject: [PATCH 2/8] fix file handle leak --- src/angle_write.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/angle_write.cpp b/src/angle_write.cpp index e0d6827a66..ca0118eda5 100644 --- a/src/angle_write.cpp +++ b/src/angle_write.cpp @@ -147,6 +147,7 @@ void AngleWrite::command(int narg, char **arg) fgets(line, MAXLINE, coeffs); writer->input->one(fmt::format("angle_coeff {}", line)); } + fclose(coeffs); platform::unlink(coeffs_file); // initialize system From ecf11f2f20d16de4a81bed70efc6e21090148dbb Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 25 Dec 2022 06:06:34 -0500 Subject: [PATCH 3/8] use macro for keeping repetitive code consistent --- src/angle_write.cpp | 55 +++++++++++++-------------------------------- 1 file changed, 15 insertions(+), 40 deletions(-) diff --git a/src/angle_write.cpp b/src/angle_write.cpp index ca0118eda5..554f03d18f 100644 --- a/src/angle_write.cpp +++ b/src/angle_write.cpp @@ -170,23 +170,19 @@ void AngleWrite::command(int narg, char **arg) force->angle_style, atype); fmt::print(fp, "\n{}\nN {} EQ {:.15g}\n\n", keyword, n, theta0); +#define GET_ENERGY(myphi, mytheta) \ + theta = mytheta; \ + atom3[0] = cos(theta * DEG2RAD); \ + atom3[1] = sin(theta * DEG2RAD); \ + myphi = angle->single(atype, i2, i1, i3) + 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); + GET_ENERGY(phi, 0.0); + GET_ENERGY(phi1, epsilon); + GET_ENERGY(phi2, 2.0 * epsilon); f = (1.5 * phi - 2.0 * phi1 + 0.5 * phi2) / epsilon; if (!std::isfinite(f)) f = 0.0; @@ -194,20 +190,10 @@ void AngleWrite::command(int narg, char **arg) 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); + GET_ENERGY(phi1, dtheta * static_cast(i) - epsilon); + GET_ENERGY(phi2, dtheta * static_cast(i) + epsilon); + GET_ENERGY(phi, dtheta * static_cast(i)); - 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 @@ -217,20 +203,9 @@ void AngleWrite::command(int narg, char **arg) } // 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); + GET_ENERGY(phi, 180.0); + GET_ENERGY(phi1, 180.0 - epsilon); + GET_ENERGY(phi2, 180.0 - 2.0 * epsilon); f = (2.0 * phi1 - 1.5 * phi - 0.5 * phi2) / epsilon; if (!std::isfinite(f)) f = 0.0; From c4f2befb1f535daf06f42c7bebe98ee3c4ba21c9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 25 Dec 2022 06:25:45 -0500 Subject: [PATCH 4/8] add sanity check on valid angle type --- src/angle_write.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/angle_write.cpp b/src/angle_write.cpp index 554f03d18f..3b7c7f2e24 100644 --- a/src/angle_write.cpp +++ b/src/angle_write.cpp @@ -56,6 +56,9 @@ void AngleWrite::command(int narg, char **arg) if (narg != 4) error->all(FLERR, "Illegal angle_write command"); int atype = utils::inumeric(FLERR, arg[0], false, lmp); + if ((atype <= 0) || (atype > atom->nangletypes)) + error->all(FLERR, "Invalid angle type {} in angle_write command", atype); + int n = utils::inumeric(FLERR, arg[1], false, lmp); std::string table_file = arg[2]; std::string keyword = arg[3]; From a4f8cb9a9230566d72d579e7bdb0a260130d3ef7 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 25 Dec 2022 06:54:40 -0500 Subject: [PATCH 5/8] explicitly disallow angle_write with angle_style class2 --- src/angle_write.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/angle_write.cpp b/src/angle_write.cpp index 3b7c7f2e24..e7c95be3cf 100644 --- a/src/angle_write.cpp +++ b/src/angle_write.cpp @@ -51,6 +51,10 @@ void AngleWrite::command(int narg, char **arg) if (angle && (force->angle->writedata == 0)) error->all(FLERR, "Angle style must support writing coeffs to data file for angle_write"); + if (angle && (utils::strmatch(force->angle_style, "^class2"))) + error->all(FLERR, "Angle_write command is not compatible with angle style {}", + force->angle_style); + // parse arguments if (narg != 4) error->all(FLERR, "Illegal angle_write command"); From bbfc7381fb81dca21ffc9b109e6761150cc4d63a Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 25 Dec 2022 06:54:53 -0500 Subject: [PATCH 6/8] updates and corrections for docs --- doc/src/angle_write.rst | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/doc/src/angle_write.rst b/doc/src/angle_write.rst index 16a4ee4e2c..8bc94b5191 100644 --- a/doc/src/angle_write.rst +++ b/doc/src/angle_write.rst @@ -28,10 +28,10 @@ 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 `. +force with respect to the angle, 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 @@ -69,7 +69,7 @@ separate, internally created, new LAMMPS instance with a dummy system of 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 +consequence of 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. @@ -77,20 +77,19 @@ following conditions must be met: :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 - `. + :doc:`angle style class2 ` 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 if the requirement is not met. It is thus recommended to make + careful tests for any created tables. Related commands """""""""""""""" :doc:`angle_style table `, :doc:`bond_write `, -:doc:`angle_style `, :doc:`angle_coeff ` +:doc:`dihedral_write `, :doc:`angle_style `, +:doc:`angle_coeff ` Default """"""" From 4ac830bf73a90820070b96c54057c36a829fd3c9 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 25 Dec 2022 06:57:38 -0500 Subject: [PATCH 7/8] add dihedral_write command --- doc/src/Commands_all.rst | 1 + doc/src/commands_list.rst | 1 + doc/src/dihedral_write.rst | 99 ++++++++++ doc/utils/sphinx-config/false_positives.txt | 1 + src/dihedral_write.cpp | 206 ++++++++++++++++++++ src/dihedral_write.h | 34 ++++ 6 files changed, 342 insertions(+) create mode 100644 doc/src/dihedral_write.rst create mode 100644 src/dihedral_write.cpp create mode 100644 src/dihedral_write.h diff --git a/doc/src/Commands_all.rst b/doc/src/Commands_all.rst index 2ba5cba51d..4a035f5b41 100644 --- a/doc/src/Commands_all.rst +++ b/doc/src/Commands_all.rst @@ -46,6 +46,7 @@ table above. * :doc:`dielectric ` * :doc:`dihedral_coeff ` * :doc:`dihedral_style ` + * :doc:`dihedral_write ` * :doc:`dimension ` * :doc:`displace_atoms ` * :doc:`dump ` diff --git a/doc/src/commands_list.rst b/doc/src/commands_list.rst index 5731cbb755..cea237f14e 100644 --- a/doc/src/commands_list.rst +++ b/doc/src/commands_list.rst @@ -28,6 +28,7 @@ Commands dielectric dihedral_coeff dihedral_style + dihedral_write dimension displace_atoms dump diff --git a/doc/src/dihedral_write.rst b/doc/src/dihedral_write.rst new file mode 100644 index 0000000000..d0deb1d2e8 --- /dev/null +++ b/doc/src/dihedral_write.rst @@ -0,0 +1,99 @@ +.. index:: dihedral_write + +dihedral_write command +====================== + +Syntax +"""""" + +.. code-block:: LAMMPS + + dihedral_write dtype N file keyword + +* dtype = dihedral 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 + + dihedral_write 1 500 table.txt Harmonic_1 + dihedral_write 3 1000 table.txt Harmonic_3 + +Description +""""""""""" + +Write energy and force values to a file as a function of the dihedral +angle for the currently defined dihedral potential. Force in this +context means the force with respect to the dihedral angle, 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:`dihedral 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 dihedrals ranging from 0 +degrees to 360 degrees for 4 interacting atoms forming an dihedral type +dtype, using the appropriate :doc:`dihedral_coeff ` +coefficients. N evenly spaced dihedrals are used. Since 0 and 360 +degrees are the same dihedral angle, the latter entry is skipped. + +For example, for N = 6, values would be computed at +:math:`\phi = 0, 60, 120, 180, 240, 300`. + +The file is written in the format used as input for the +:doc:`dihedral_style table ` option with *keyword* as +the section name. Each line written to the file lists an index number +(1-N), an dihedral 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 *dihedral_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 *dihedral_write* will refuse to add a +table to an existing file if the units are not the same. + +Restrictions +"""""""""""" + +All force field coefficients for dihedrals and other kinds of interactions +must be set before this command can be invoked. + +The table of the dihedral energy and force data data is created by using a +separate, internally created, new LAMMPS instance with a dummy system of +4 atoms for which the dihedral potential energy is computed after +transferring the dihedral style and coefficients and arranging the 4 atoms +into the corresponding geometries. The dihedral force is then determined +from the potential energies through numerical differentiation. As a +consequence of this approach, not all dihedral styles are compatible. The +following conditions must be met: + +- The dihedral style must be able to write its coefficients to a data file. + This condition excludes for example :doc:`dihedral style hybrid ` and + :doc:`dihedral style table `. +- The potential function must not have any terms that depend on geometry + properties other than the dihedral. This condition excludes for + example :doc:`dihedral style class2 `. Please note + that the *write_dihedral* command has no way of checking for this + condition. It will check the style name against an internal list of + known to be incompatible styles. The resulting tables may be bogus + for unlisted dihedral styles if the requirement is not met. It is + thus recommended to make careful tests for any created tables. + +Related commands +"""""""""""""""" + +:doc:`dihedral_style table `, :doc:`bond_write `, +:doc:`angle_write `, :doc:`dihedral_style `, +:doc:`dihedral_coeff ` + +Default +""""""" + +none diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 13cc96f993..855f306f84 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -829,6 +829,7 @@ dtemp dtgrow dtheta dtshrink +dtype du dU Ducastelle diff --git a/src/dihedral_write.cpp b/src/dihedral_write.cpp new file mode 100644 index 0000000000..eedcde0eb8 --- /dev/null +++ b/src/dihedral_write.cpp @@ -0,0 +1,206 @@ +/* ---------------------------------------------------------------------- + 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 "dihedral_write.h" + +#include "atom.h" +#include "atom_vec.h" +#include "comm.h" +#include "dihedral.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 DihedralWrite::command(int narg, char **arg) +{ + // sanity checks + + if (domain->box_exist == 0) + error->all(FLERR, "Dihedral_write command before simulation box is defined"); + if (atom->avec->dihedrals_allow == 0) + error->all(FLERR, "Dihedral_write command when no dihedrals allowed"); + auto dihedral = force->dihedral; + if (dihedral == nullptr) + error->all(FLERR, "Dihedral_write command before an dihedral_style is defined"); + if (dihedral && (force->dihedral->writedata == 0)) + error->all(FLERR, "Dihedral style must support writing coeffs to data file for dihedral_write"); + + if (dihedral && + (utils::strmatch(force->dihedral_style, "^charmm") || + utils::strmatch(force->dihedral_style, "^class2"))) + error->all(FLERR, "Dihedral_write command is not compatible with dihedral style {}", + force->dihedral_style); + + // parse arguments + + if (narg != 4) error->all(FLERR, "Illegal dihedral_write command"); + + int dtype = utils::inumeric(FLERR, arg[0], false, lmp); + if ((dtype <= 0) || (dtype > atom->ndihedraltypes)) + error->all(FLERR, "Invalid dihedral type {} in dihedral_write command", dtype); + + 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(); + + // write out all dihedral_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 dihedral_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->dihedral->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 dihedral_write\n", utils::current_date(), + update->unit_style); + } + if (fp == nullptr) + error->one(FLERR, "Cannot open dihedral_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 dihedral potential + // const char *args[] = {"DihedralWrite", "-nocite", "-echo", "none", + // "-log", "none", "-screen", "none"}; + const char *args[] = {"DihedralWrite", "-nocite", "-echo", "screen", "-log", "none"}; + char **argv = (char **) args; + int argc = sizeof(args) / sizeof(char *); + LAMMPS *writer = new LAMMPS(argc, argv, singlecomm); + + // create dummy system replicating dihedral style settings + writer->input->one(fmt::format("units {}", update->unit_style)); + writer->input->one("atom_style molecular"); + writer->input->one("atom_modify map array"); + writer->input->one("boundary f f f"); + writer->input->one("region box block -2 2 -2 2 -2 2"); + writer->input->one(fmt::format("create_box {} box dihedral/types {} " + "extra/dihedral/per/atom 1 " + "extra/special/per/atom 4", + atom->ntypes, atom->ndihedraltypes)); + writer->input->one("create_atoms 1 single 1.0 0.0 -1.0"); + writer->input->one("create_atoms 1 single 0.0 0.0 -1.0"); + writer->input->one("create_atoms 1 single 0.0 0.0 1.0"); + writer->input->one("create_atoms 1 single 1.0 0.0 1.0"); + writer->input->one(fmt::format("create_bonds single/dihedral {} 1 2 3 4", dtype)); + + writer->input->one("pair_style zero 10.0"); + writer->input->one("pair_coeff * *"); + writer->input->one("mass * 1.0"); + writer->input->one(fmt::format("dihedral_style {}", force->dihedral_style)); + FILE *coeffs; + char line[MAXLINE]; + coeffs = fopen(coeffs_file.c_str(), "r"); + for (int i = 0; i < atom->ndihedraltypes; ++i) { + fgets(line, MAXLINE, coeffs); + writer->input->one(fmt::format("dihedral_coeff {}", line)); + } + fclose(coeffs); + platform::unlink(coeffs_file); + + // must complete a full setup() to initialize system including neighbor and dihedral lists. + + writer->input->one("run 0 post no"); + + // move third atom to reproduce dihedrals + + double theta, phi, phi1, phi2, f; + dihedral = writer->force->dihedral; + auto atom4 = writer->atom->x[writer->atom->map(4)]; + + // evaluate energy and force at each of N distances + + fmt::print(fp, "# Dihedral potential {} for dihedral type {}: i,theta,energy,force\n", + force->dihedral_style, dtype); + fmt::print(fp, "\n{}\nN {} DEGREES\n\n", keyword, n); + +#define GET_ENERGY(myphi, mytheta) \ + theta = mytheta; \ + atom4[0] = cos(theta * DEG2RAD); \ + atom4[1] = sin(theta * DEG2RAD); \ + dihedral->energy = 0.0; \ + dihedral->compute(ENERGY_GLOBAL, 0); \ + myphi = dihedral->energy + + const double dtheta = 360.0 / static_cast(n); + for (int i = 0; i < n; i++) { + GET_ENERGY(phi1, dtheta * static_cast(i) - epsilon); + GET_ENERGY(phi2, dtheta * static_cast(i) + epsilon); + GET_ENERGY(phi, dtheta * static_cast(i)); + + 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); + } + + // clean up + delete writer; + fclose(fp); + } + MPI_Comm_free(&singlecomm); +} diff --git a/src/dihedral_write.h b/src/dihedral_write.h new file mode 100644 index 0000000000..b17f8d19da --- /dev/null +++ b/src/dihedral_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(dihedral_write,DihedralWrite); +// clang-format on +#else + +#ifndef LMP_DIHEDRAL_WRITE_H +#define LMP_DIHEDRAL_WRITE_H + +#include "command.h" + +namespace LAMMPS_NS { + +class DihedralWrite : public Command { + public: + DihedralWrite(class LAMMPS *lmp) : Command(lmp){}; + void command(int, char **) override; +}; +} // namespace LAMMPS_NS +#endif +#endif From 63ddb07c594fd7124a20df849f4e6b6ec05d2698 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 25 Dec 2022 07:04:20 -0500 Subject: [PATCH 8/8] add versionadded tags --- doc/src/angle_write.rst | 2 ++ doc/src/dihedral_write.rst | 2 ++ 2 files changed, 4 insertions(+) diff --git a/doc/src/angle_write.rst b/doc/src/angle_write.rst index 8bc94b5191..1541a7120a 100644 --- a/doc/src/angle_write.rst +++ b/doc/src/angle_write.rst @@ -26,6 +26,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + 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 force with respect to the angle, not the force on individual atoms. diff --git a/doc/src/dihedral_write.rst b/doc/src/dihedral_write.rst index d0deb1d2e8..94fa3da0fa 100644 --- a/doc/src/dihedral_write.rst +++ b/doc/src/dihedral_write.rst @@ -26,6 +26,8 @@ Examples Description """"""""""" +.. versionadded:: TBD + Write energy and force values to a file as a function of the dihedral angle for the currently defined dihedral potential. Force in this context means the force with respect to the dihedral angle, not the