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

implement angle_write command

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-12-24 22:49:42 -05:00
parent 07fe2fa29d
commit da98363a25
8 changed files with 387 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/Commands_all.rst
View File

@ -24,6 +24,7 @@ table above.
* :doc:`angle_coeff <angle_coeff>` * :doc:`angle_coeff <angle_coeff>`
* :doc:`angle_style <angle_style>` * :doc:`angle_style <angle_style>`
* :doc:`angle_write <angle_write>`
* :doc:`atom_modify <atom_modify>` * :doc:`atom_modify <atom_modify>`
* :doc:`atom_style <atom_style>` * :doc:`atom_style <atom_style>`
* :doc:`balance <balance>` * :doc:`balance <balance>`

11
doc/src/angle_table.rst
View File

@ -59,9 +59,12 @@ format of this file is described below.
---------- ----------
Suitable tables for use with this angle style can be created using the Suitable tables for use with this angle style can be created by LAMMPS
Python code in the ``tools/tabulate`` folder of the LAMMPS source code itself from existing angle styles using the :doc:`angle_write
distribution. <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 The format of a tabulated file is as follows (without the
parenthesized comments): parenthesized comments):
@ -154,7 +157,7 @@ for more info.
Related commands Related commands
"""""""""""""""" """"""""""""""""
:doc:`angle_coeff <angle_coeff>` :doc:`angle_coeff <angle_coeff>`, :doc:`angle_write <angle_write>`
Default Default
""""""" """""""

98
doc/src/angle_write.rst Normal file
View File

@ -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 <angle_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 <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 <angle_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 <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 <angle_hybrid>` and
:doc:`angle style table <angle_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 <angle_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
<angle_class2>`.
Related commands
""""""""""""""""
:doc:`angle_style table <angle_table>`, :doc:`bond_write <bond_write>`,
:doc:`angle_style <angle_style>`, :doc:`angle_coeff <angle_coeff>`
Default
"""""""
none

2
doc/src/bond_write.rst
View File

@ -67,7 +67,7 @@ be specified even if the potential has a finite value at r = 0.0.
Related commands Related commands
"""""""""""""""" """"""""""""""""
:doc:`bond_style table <bond_table>`, :doc:`bond_style table <bond_table>`, `angle_write <angle_write>`,
:doc:`bond_style <bond_style>`, :doc:`bond_coeff <bond_coeff>` :doc:`bond_style <bond_style>`, :doc:`bond_coeff <bond_coeff>`
Default Default

1
doc/src/commands_list.rst
View File

@ -6,6 +6,7 @@ Commands
angle_coeff angle_coeff
angle_style angle_style
angle_write
atom_modify atom_modify
atom_style atom_style
balance balance

1
doc/utils/sphinx-config/false_positives.txt
View File

@ -174,6 +174,7 @@ attrac
Atw Atw
Atwater Atwater
atwt atwt
atype
augt augt
AuO AuO
automagically automagically

244
src/angle_write.cpp Normal file
View File

@ -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 <cmath>
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<double>(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<double>(i) - epsilon;
atom3[0] = cos(theta * DEG2RAD);
atom3[1] = sin(theta * DEG2RAD);
phi1 = angle->single(atype, i2, i1, i3);
theta = dtheta * static_cast<double>(i) + epsilon;
atom3[0] = cos(theta * DEG2RAD);
atom3[1] = sin(theta * DEG2RAD);
phi2 = angle->single(atype, i2, i1, i3);
theta = dtheta * static_cast<double>(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);
}

34
src/angle_write.h Normal file
View File

@ -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