diff --git a/doc/src/Commands_all.rst b/doc/src/Commands_all.rst index bd9d0ed09c..4a035f5b41 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 ` @@ -45,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/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..1541a7120a --- /dev/null +++ b/doc/src/angle_write.rst @@ -0,0 +1,99 @@ +.. 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 +""""""""""" + +.. 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. +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 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. + 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 + :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:`dihedral_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..cea237f14e 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 @@ -27,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..94fa3da0fa --- /dev/null +++ b/doc/src/dihedral_write.rst @@ -0,0 +1,101 @@ +.. 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 +""""""""""" + +.. 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 +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 7cb409d040..855f306f84 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 @@ -828,6 +829,7 @@ dtemp dtgrow dtheta dtshrink +dtype du dU Ducastelle diff --git a/src/angle_write.cpp b/src/angle_write.cpp new file mode 100644 index 0000000000..e7c95be3cf --- /dev/null +++ b/src/angle_write.cpp @@ -0,0 +1,227 @@ +/* ---------------------------------------------------------------------- + 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"); + + 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"); + + 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]; + + // 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)); + } + fclose(coeffs); + 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); + +#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 + + 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; + 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++) { + 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); + } + + // get force for divergent 180 degree angle from interpolation to the left + 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; + 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 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