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

new math constants: DEG2RAD and RAD2DEG

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-03-17 19:39:22 -04:00
parent 6e7eedf148
commit 58030009da
14 changed files with 106 additions and 110 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
src/MOLECULE/angle_charmm.cpp
View File

@ -29,7 +29,8 @@
#include <cmath>
using namespace LAMMPS_NS;
using MathConst::MY_PI;
using MathConst::DEG2RAD;
using MathConst::RAD2DEG;
static constexpr double SMALL = 0.001;
@ -211,7 +212,7 @@ void AngleCharmm::coeff(int narg, char **arg)
int count = 0;
for (int i = ilo; i <= ihi; i++) {
k[i] = k_one;
theta0[i] = theta0_one / 180.0 * MY_PI;
theta0[i] = DEG2RAD * theta0_one;
k_ub[i] = k_ub_one;
r_ub[i] = r_ub_one;
setflag[i] = 1;
@ -269,7 +270,7 @@ void AngleCharmm::read_restart(FILE *fp)
void AngleCharmm::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nangletypes; i++)
fprintf(fp, "%d %g %g %g %g\n", i, k[i], theta0[i] / MY_PI * 180.0, k_ub[i], r_ub[i]);
fprintf(fp, "%d %g %g %g %g\n", i, k[i], RAD2DEG * theta0[i], k_ub[i], r_ub[i]);
}
/* ---------------------------------------------------------------------- */

7
src/MOLECULE/angle_cosine_squared.cpp
View File

@ -29,7 +29,8 @@
#include <cmath>
using namespace LAMMPS_NS;
using MathConst::MY_PI;
using MathConst::DEG2RAD;
using MathConst::RAD2DEG;
/* ---------------------------------------------------------------------- */
@ -181,7 +182,7 @@ void AngleCosineSquared::coeff(int narg, char **arg)
int count = 0;
for (int i = ilo; i <= ihi; i++) {
k[i] = k_one;
theta0[i] = theta0_one / 180.0 * MY_PI;
theta0[i] = DEG2RAD * theta0_one;
setflag[i] = 1;
count++;
}
@ -231,7 +232,7 @@ void AngleCosineSquared::read_restart(FILE *fp)
void AngleCosineSquared::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nangletypes; i++)
fprintf(fp, "%d %g %g\n", i, k[i], theta0[i] / MY_PI * 180.0);
fprintf(fp, "%d %g %g\n", i, k[i], RAD2DEG * theta0[i]);
}
/* ---------------------------------------------------------------------- */

7
src/MOLECULE/angle_harmonic.cpp
View File

@ -25,7 +25,8 @@
#include <cmath>
using namespace LAMMPS_NS;
using namespace MathConst;
using MathConst::DEG2RAD;
using MathConst::RAD2DEG;
static constexpr double SMALL = 0.001;
@ -183,7 +184,7 @@ void AngleHarmonic::coeff(int narg, char **arg)
int count = 0;
for (int i = ilo; i <= ihi; i++) {
k[i] = k_one;
theta0[i] = theta0_one / 180.0 * MY_PI;
theta0[i] = DEG2RAD * theta0_one;
setflag[i] = 1;
count++;
}
@ -233,7 +234,7 @@ void AngleHarmonic::read_restart(FILE *fp)
void AngleHarmonic::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nangletypes; i++)
fprintf(fp, "%d %g %g\n", i, k[i], theta0[i] / MY_PI * 180.0);
fprintf(fp, "%d %g %g\n", i, k[i], RAD2DEG * theta0[i]);
}
/* ---------------------------------------------------------------------- */

12
src/MOLECULE/angle_table.cpp
View File

@ -32,7 +32,9 @@
#include <cstring>
using namespace LAMMPS_NS;
using MathConst::DEG2RAD;
using MathConst::MY_PI;
using MathConst::RAD2DEG;
enum { LINEAR, SPLINE };
@ -243,8 +245,8 @@ void AngleTable::coeff(int narg, char **arg)
// convert theta from degrees to radians
for (int i = 0; i < tb->ninput; i++) {
tb->afile[i] *= MY_PI / 180.0;
tb->ffile[i] *= 180.0 / MY_PI;
tb->afile[i] *= DEG2RAD;
tb->ffile[i] *= RAD2DEG;
}
// spline read-in and compute a,e,f vectors within table
@ -519,10 +521,10 @@ void AngleTable::param_extract(Table *tb, char *line)
tb->fpflag = 1;
tb->fplo = values.next_double();
tb->fphi = values.next_double();
tb->fplo *= (180.0 / MY_PI) * (180.0 / MY_PI);
tb->fphi *= (180.0 / MY_PI) * (180.0 / MY_PI);
tb->fplo *= RAD2DEG * RAD2DEG;
tb->fphi *= RAD2DEG * RAD2DEG;
} else if (word == "EQ") {
tb->theta0 = values.next_double() / 180.0 * MY_PI;
tb->theta0 = DEG2RAD * values.next_double();
} else {
error->one(FLERR, "Invalid keyword in angle table parameters");
}

36
src/angle_zero.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -20,18 +19,19 @@
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
using MathConst::DEG2RAD;
using MathConst::RAD2DEG;
/* ---------------------------------------------------------------------- */
AngleZero::AngleZero(LAMMPS *lmp) : Angle(lmp), coeffflag(1) {}
AngleZero::AngleZero(LAMMPS *_lmp) : Angle(_lmp), coeffflag(1) {}
/* ---------------------------------------------------------------------- */
@ -54,12 +54,13 @@ void AngleZero::compute(int eflag, int vflag)
void AngleZero::settings(int narg, char **arg)
{
if ((narg != 0) && (narg != 1))
error->all(FLERR,"Illegal angle_style command");
if ((narg != 0) && (narg != 1)) error->all(FLERR, "Illegal angle_style command");
if (narg == 1) {
if (strcmp("nocoeff",arg[0]) == 0) coeffflag=0;
else error->all(FLERR,"Illegal angle_style command");
if (strcmp("nocoeff", arg[0]) == 0)
coeffflag = 0;
else
error->all(FLERR, "Illegal angle_style command");
}
}
@ -68,11 +69,11 @@ void AngleZero::settings(int narg, char **arg)
void AngleZero::allocate()
{
allocated = 1;
int n = atom->nangletypes;
const int np1 = atom->nangletypes + 1;
memory->create(theta0,n+1,"angle:theta0");
memory->create(setflag,n+1,"angle:setflag");
for (int i = 1; i <= n; i++) setflag[i] = 0;
memory->create(theta0, np1, "angle:theta0");
memory->create(setflag, np1, "angle:setflag");
for (int i = 1; i < np1; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
@ -90,15 +91,14 @@ void AngleZero::coeff(int narg, char **arg)
utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
double theta0_one = 0.0;
if (coeffflag && (narg == 2))
theta0_one = utils::numeric(FLERR,arg[1],false,lmp);
if (coeffflag && (narg == 2)) theta0_one = utils::numeric(FLERR, arg[1], false, lmp);
// convert theta0 from degrees to radians
int count = 0;
for (int i = ilo; i <= ihi; i++) {
setflag[i] = 1;
theta0[i] = theta0_one/180.0 * MY_PI;
theta0[i] = DEG2RAD * theta0_one;
count++;
}
@ -116,7 +116,8 @@ double AngleZero::equilibrium_angle(int i)
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void AngleZero::write_restart(FILE *fp) {
void AngleZero::write_restart(FILE *fp)
{
fwrite(&theta0[1], sizeof(double), atom->nangletypes, fp);
}
@ -141,8 +142,7 @@ void AngleZero::read_restart(FILE *fp)
void AngleZero::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nangletypes; i++)
fprintf(fp,"%d %g\n",i,theta0[i]/MY_PI*180.0);
for (int i = 1; i <= atom->nangletypes; i++) fprintf(fp, "%d %g\n", i, RAD2DEG * theta0[i]);
}
/* ---------------------------------------------------------------------- */

7
src/displace_atoms.cpp
View File

@ -37,13 +37,14 @@
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
using MathConst::DEG2RAD;
using MathConst::MY_2PI;
enum{MOVE,RAMP,RANDOM,ROTATE};
/* ---------------------------------------------------------------------- */
DisplaceAtoms::DisplaceAtoms(LAMMPS *lmp) : Command(lmp)
DisplaceAtoms::DisplaceAtoms(LAMMPS *_lmp) : Command(_lmp)
{
mvec = nullptr;
}
@ -235,7 +236,7 @@ void DisplaceAtoms::command(int narg, char **arg)
runit[1] = axis[1]/len;
runit[2] = axis[2]/len;
double angle = MY_PI*theta/180.0;
double angle = DEG2RAD*theta;
double cosine = cos(angle);
double sine = sin(angle);

22
src/dump_image.cpp
View File

@ -40,7 +40,7 @@
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
using MathConst::DEG2RAD;
#define BIG 1.0e20
@ -228,18 +228,15 @@ DumpImage::DumpImage(LAMMPS *lmp, int narg, char **arg) :
if (utils::strmatch(arg[iarg+1],"^v_")) {
thetastr = utils::strdup(arg[iarg+1]+2);
} else {
double theta = utils::numeric(FLERR,arg[iarg+1],false,lmp);
const double theta = utils::numeric(FLERR,arg[iarg+1],false,lmp);
if (theta < 0.0 || theta > 180.0)
error->all(FLERR,"Invalid dump image theta value");
theta *= MY_PI/180.0;
image->theta = theta;
image->theta = DEG2RAD * theta;
}
if (utils::strmatch(arg[iarg+2],"^v_")) {
phistr = utils::strdup(arg[iarg+2]+2);
} else {
double phi = utils::numeric(FLERR,arg[iarg+2],false,lmp);
phi *= MY_PI/180.0;
image->phi = phi;
image->phi = DEG2RAD * utils::numeric(FLERR,arg[iarg+2],false,lmp);
}
iarg += 3;
@ -652,18 +649,13 @@ void DumpImage::view_params()
// view direction theta and phi
if (thetastr) {
double theta = input->variable->compute_equal(thetavar);
const double theta = input->variable->compute_equal(thetavar);
if (theta < 0.0 || theta > 180.0)
error->all(FLERR,"Invalid dump image theta value");
theta *= MY_PI/180.0;
image->theta = theta;
image->theta = DEG2RAD * theta;
}
if (phistr) {
double phi = input->variable->compute_equal(phivar);
phi *= MY_PI/180.0;
image->phi = phi;
}
if (phistr) image->phi = DEG2RAD * input->variable->compute_equal(phivar);
// up vector

17
src/fix_gravity.cpp
View File

@ -29,14 +29,14 @@
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
using MathConst::DEG2RAD;
enum{CHUTE,SPHERICAL,VECTOR};
/* ---------------------------------------------------------------------- */
FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
FixGravity::FixGravity(LAMMPS *_lmp, int narg, char **arg) :
Fix(_lmp, narg, arg),
mstr(nullptr), vstr(nullptr), pstr(nullptr), tstr(nullptr),
xstr(nullptr), ystr(nullptr), zstr(nullptr)
{
@ -135,7 +135,6 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) :
// initializations
degree2rad = MY_PI/180.0;
time_origin = update->ntimestep;
eflag = 0;
@ -328,12 +327,12 @@ void FixGravity::set_acceleration()
theta = 180.0 - vert;
}
if (domain->dimension == 3) {
xgrav = sin(degree2rad * theta) * cos(degree2rad * phi);
ygrav = sin(degree2rad * theta) * sin(degree2rad * phi);
zgrav = cos(degree2rad * theta);
xgrav = sin(DEG2RAD * theta) * cos(DEG2RAD * phi);
ygrav = sin(DEG2RAD * theta) * sin(DEG2RAD * phi);
zgrav = cos(DEG2RAD * theta);
} else {
xgrav = sin(degree2rad * theta);
ygrav = cos(degree2rad * theta);
xgrav = sin(DEG2RAD * theta);
ygrav = cos(DEG2RAD * theta);
zgrav = 0.0;
}
} else if (style == VECTOR) {

1
src/fix_gravity.h
View File

@ -45,7 +45,6 @@ class FixGravity : public Fix {
double vert, phi, theta;
double xdir, ydir, zdir;
double xgrav, ygrav, zgrav, xacc, yacc, zacc;
double degree2rad;
int ilevel_respa;
int time_origin;
double gvec[3];

8
src/fix_restrain.cpp
View File

@ -33,7 +33,7 @@
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
using MathConst::DEG2RAD;
enum{BOND,LBOUND,ANGLE,DIHEDRAL};
@ -119,8 +119,7 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
ids[nrestrain][2] = utils::tnumeric(FLERR,arg[iarg+3],false,lmp);
kstart[nrestrain] = utils::numeric(FLERR,arg[iarg+4],false,lmp);
kstop[nrestrain] = utils::numeric(FLERR,arg[iarg+5],false,lmp);
target[nrestrain] = utils::numeric(FLERR,arg[iarg+6],false,lmp);
target[nrestrain] *= MY_PI / 180.0;
target[nrestrain] = DEG2RAD * utils::numeric(FLERR,arg[iarg+6],false,lmp);
iarg += 7;
} else if (strcmp(arg[iarg],"dihedral") == 0) {
if (iarg+8 > narg) error->all(FLERR,"Illegal fix restrain command");
@ -132,8 +131,7 @@ FixRestrain::FixRestrain(LAMMPS *lmp, int narg, char **arg) :
ids[nrestrain][3] = utils::tnumeric(FLERR,arg[iarg+4],false,lmp);
kstart[nrestrain] = utils::numeric(FLERR,arg[iarg+5],false,lmp);
kstop[nrestrain] = utils::numeric(FLERR,arg[iarg+6],false,lmp);
target[nrestrain] = utils::numeric(FLERR,arg[iarg+7],false,lmp);
target[nrestrain] *= MY_PI / 180.0;
target[nrestrain] = DEG2RAD * utils::numeric(FLERR,arg[iarg+7],false,lmp);
cos_target[nrestrain] = cos(target[nrestrain]);
sin_target[nrestrain] = sin(target[nrestrain]);
iarg += 8;

8
src/image.cpp
View File

@ -40,7 +40,9 @@
#endif
using namespace LAMMPS_NS;
using namespace MathConst;
using MathConst::DEG2RAD;
using MathConst::MY_PI;
using MathConst::MY_PI4;
#define NCOLORS 140
#define NELEMENTS 109
@ -61,8 +63,8 @@ Image::Image(LAMMPS *lmp, int nmap_caller) : Pointers(lmp)
// defaults for 3d viz
width = height = 512;
theta = 60.0 * MY_PI/180.0;
phi = 30.0 * MY_PI/180.0;
theta = 60.0 * DEG2RAD;
phi = 30.0 * DEG2RAD;
zoom = 1.0;
persp = 0.0;
shiny = 1.0;

2
src/math_const.h
View File

@ -30,6 +30,8 @@ namespace MathConst {
static constexpr double MY_SQRT2 = 1.41421356237309504880; // sqrt(2)
static constexpr double MY_ISQRT2 = 0.707106781186547524401; // 1/sqrt(2)
static constexpr double MY_CUBEROOT2 = 1.25992104989487316476; // 2*(1/3)
static constexpr double DEG2RAD = MY_PI / 180.0; // degree to radians
static constexpr double RAD2DEG = 180.0 / MY_PI; // radians to degree
} // namespace MathConst
} // namespace LAMMPS_NS

5
src/set.cpp
View File

@ -321,10 +321,7 @@ void Set::command(int narg, char **arg)
} else if (strcmp(arg[iarg],"theta") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal set command");
if (utils::strmatch(arg[iarg+1],"^v_")) varparse(arg[iarg+1],1);
else {
dvalue = utils::numeric(FLERR,arg[iarg+1],false,lmp);
dvalue *= MY_PI/180.0;
}
else dvalue = DEG2RAD * utils::numeric(FLERR,arg[iarg+1],false,lmp);
if (!atom->line_flag)
error->all(FLERR,"Cannot set this attribute for this atom style");
set(THETA);

59
unittest/force-styles/tests/fix-timestep-spring_rg.yaml
View File

@ -1,7 +1,8 @@
---
lammps_version: 10 Feb 2021
date_generated: Fri Feb 26 23:09:58 2021
lammps_version: 17 Feb 2022
date_generated: Thu Mar 17 19:43:17 2022
epsilon: 2e-14
skip_tests:
prerequisites: ! |
atom full
fix spring/rg
@ -11,9 +12,9 @@ post_commands: ! |
fix test solute spring/rg 10.0 NULL
input_file: in.fourmol
natoms: 29
global_scalar: 2.28284405113276
global_scalar: 2.2828440511327632
run_pos: ! |2
1 -2.7045538630662386e-01 2.4912156114977786e+00 -1.6695860623999764e-01
1 -2.7045538630662386e-01 2.4912156114977786e+00 -1.6695860623999767e-01
2 3.1004042556254785e-01 2.9612350151324356e+00 -8.5466362489112224e-01
3 -7.0398524408175600e-01 1.2305507880027069e+00 -6.2777529566908619e-01
4 -1.5818155384081236e+00 1.4837405418512759e+00 -1.2538710203031622e+00
@ -43,33 +44,33 @@ run_pos: ! |2
28 -2.7473562684513411e+00 -4.0200819932379330e+00 1.5830052163433954e+00
29 -1.3126000191359855e+00 -3.5962518039482929e+00 2.2746342468737835e+00
run_vel: ! |2
1 8.1709752174462542e-03 1.6515688596142518e-02 4.7900602931841710e-03
2 5.4503971185775359e-03 5.1783216918025263e-03 -1.4372851752161298e-03
3 -8.2293181444702025e-03 -1.2926946755956580e-02 -4.0984684087231526e-03
4 -3.7691559672827545e-03 -6.5727436126565286e-03 -1.1183447702477240e-03
5 -1.1021393058885439e-02 -9.8909848370648083e-03 -2.8413942770682352e-03
1 8.1709752174462542e-03 1.6515688596142514e-02 4.7900602931841684e-03
2 5.4503971185775411e-03 5.1783216918025367e-03 -1.4372851752161293e-03
3 -8.2293181444702043e-03 -1.2926946755956573e-02 -4.0984684087231457e-03
4 -3.7691559672827506e-03 -6.5727436126565338e-03 -1.1183447702477318e-03
5 -1.1021393058885434e-02 -9.8909848370648135e-03 -2.8413942770682383e-03
6 -3.9676421867697674e-02 4.6816949894224760e-02 3.7148621014876211e-02
7 9.1057372865591753e-04 -1.0128576591025795e-02 -5.1567805812137879e-02
8 7.9064819523862306e-03 -3.3506478635237690e-03 3.4557066595965268e-02
8 7.9064819523862306e-03 -3.3506478635237699e-03 3.4557066595965268e-02
9 1.5643732414084173e-03 3.7365718597172125e-03 1.5047128287408378e-02
10 2.9201230978839865e-02 -2.9249253288638165e-02 -1.5018033840725952e-02
10 2.9201230978839872e-02 -2.9249253288638168e-02 -1.5018033840725952e-02
11 -4.7837548086858278e-03 -3.7476576867801988e-03 -2.3461210747735017e-03
12 2.2691604899250286e-03 -3.4766287463205091e-04 -3.0638559638336716e-03
13 2.7524108097519047e-03 5.8172920420346171e-03 -7.9444854604175251e-04
14 3.5242391705623136e-03 -5.7939413457956208e-03 -3.9473455124488044e-03
15 -1.8552706489784355e-03 -5.8556732093022392e-03 6.2939676447533247e-03
16 1.8681112813127045e-02 -1.3261876488936336e-02 -4.5638880990403745e-02
17 -1.2896548250722183e-02 9.7532806974487893e-03 3.7296012357338384e-02
18 -8.0065794849152119e-04 -8.6270473216451689e-04 -1.4483040697105679e-03
19 1.2452390836177711e-03 -2.5061097118792130e-03 7.2998631009750618e-03
20 3.5930060229590064e-03 3.6938860309242609e-03 3.2322732687903584e-03
21 -1.4689220370648721e-03 -2.7352129763334656e-04 7.0581624213162688e-04
22 -7.0694199254627624e-03 -4.2577148924886872e-03 2.8079117614129140e-04
23 6.0446963117387039e-03 -1.4000131614805551e-03 2.5819754847001223e-03
24 3.1926367903015851e-04 -9.9445664751400781e-04 1.4999996958302617e-04
25 1.3789754516173812e-04 -4.4335894884620919e-03 -8.1808136724237357e-04
26 2.0485904035224033e-03 2.7813358633828971e-03 4.3245727149209051e-03
27 4.5604120292954369e-04 -1.0305523026882344e-03 2.1188058381008324e-04
28 -6.2544520861855133e-03 1.4127711176149574e-03 -1.8429821884795954e-03
29 6.4110631534415174e-04 3.1273432719600572e-03 3.7253671105653696e-03
12 2.2691604899250325e-03 -3.4766287463204712e-04 -3.0638559638336751e-03
13 2.7524108097519012e-03 5.8172920420346101e-03 -7.9444854604174731e-04
14 3.5242391705623084e-03 -5.7939413457956243e-03 -3.9473455124488009e-03
15 -1.8552706489784417e-03 -5.8556732093022409e-03 6.2939676447533325e-03
16 1.8681112813127041e-02 -1.3261876488936333e-02 -4.5638880990403745e-02
17 -1.2896548250722183e-02 9.7532806974487875e-03 3.7296012357338384e-02
18 -8.0065794849152119e-04 -8.6270473216451592e-04 -1.4483040697105707e-03
19 1.2452390836177679e-03 -2.5061097118792169e-03 7.2998631009750653e-03
20 3.5930060229590077e-03 3.6938860309242627e-03 3.2322732687903649e-03
21 -1.4689220370648717e-03 -2.7352129763334574e-04 7.0581624213162438e-04
22 -7.0694199254627719e-03 -4.2577148924886890e-03 2.8079117614129563e-04
23 6.0446963117387047e-03 -1.4000131614805551e-03 2.5819754847001292e-03
24 3.1926367903015803e-04 -9.9445664751400564e-04 1.4999996958302473e-04
25 1.3789754516173503e-04 -4.4335894884621032e-03 -8.1808136724237216e-04
26 2.0485904035224068e-03 2.7813358633828954e-03 4.3245727149209120e-03
27 4.5604120292954331e-04 -1.0305523026882316e-03 2.1188058381008232e-04
28 -6.2544520861855151e-03 1.4127711176149509e-03 -1.8429821884795954e-03
29 6.4110631534415575e-04 3.1273432719600524e-03 3.7253671105653718e-03
...