From 4aba9e9bb61be24e7514abbf924844e586008e52 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 Jan 2022 23:08:32 -0500 Subject: [PATCH] cosmetic and whitespace changes --- src/EXTRA-MOLECULE/dihedral_table_cut.cpp | 25 ++++-------- src/MOLECULE/angle_table.cpp | 15 +++---- src/MOLECULE/bond_table.cpp | 15 +++---- src/MOLECULE/dihedral_table.cpp | 48 +++++++---------------- src/OPENMP/angle_table_omp.cpp | 5 ++- src/OPENMP/bond_table_omp.cpp | 4 +- 6 files changed, 38 insertions(+), 74 deletions(-) diff --git a/src/EXTRA-MOLECULE/dihedral_table_cut.cpp b/src/EXTRA-MOLECULE/dihedral_table_cut.cpp index b132104bdd..522eff8626 100644 --- a/src/EXTRA-MOLECULE/dihedral_table_cut.cpp +++ b/src/EXTRA-MOLECULE/dihedral_table_cut.cpp @@ -82,20 +82,14 @@ enum { //GSL status return codes. GSL_EBADLEN = 19 }; - - // cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x, // with n control points at xa[], ya[], with parameters y2a[]. // The xa[] must be monotonically increasing and their // range should not exceed period (ie xa[n-1] < xa[0] + period). // x must lie in the range: [(xa[n-1]-period), (xa[0]+period)] // "period" is typically 2*PI. -static double cyc_splintD(double const *xa, - double const *ya, - double const *y2a, - int n, - double period, - double x) +static double cyc_splintD(double const *xa, double const *ya, double const *y2a, + int n, double period, double x) { int klo = -1; int khi = n; // (not n-1) @@ -490,8 +484,7 @@ void DihedralTableCut::coeff(int narg, char **arg) if (tb->ninput < 2) error->all(FLERR,"Invalid dihedral table length: {}",arg[5]); else if ((tb->ninput == 2) && (tabstyle == SPLINE)) - error->all(FLERR,"Invalid dihedral spline table length: {} " - "(Try linear)",arg[5]); + error->all(FLERR,"Invalid dihedral spline table length: {} (Try linear)",arg[5]); // check for monotonicity for (int i=0; i < tb->ninput-1; i++) { @@ -509,12 +502,10 @@ void DihedralTableCut::coeff(int narg, char **arg) double phihi = tb->phifile[tb->ninput-1]; if (tb->use_degrees) { if ((phihi - philo) >= 360) - error->all(FLERR,"Dihedral table angle range must be < 360 " - "degrees ({})",arg[5]); + error->all(FLERR,"Dihedral table angle range must be < 360 degrees ({})",arg[5]); } else { if ((phihi - philo) >= MY_2PI) - error->all(FLERR,"Dihedral table angle range must be < 2*PI " - "radians ({})",arg[5]); + error->all(FLERR,"Dihedral table angle range must be < 2*PI radians ({})",arg[5]); } // convert phi from degrees to radians @@ -532,9 +523,9 @@ void DihedralTableCut::coeff(int narg, char **arg) // We also want the angles to be sorted in increasing order. // This messy code fixes these problems with the user's data: { - double *phifile_tmp = new double [tb->ninput]; //temporary arrays - double *ffile_tmp = new double [tb->ninput]; //used for sorting - double *efile_tmp = new double [tb->ninput]; + double *phifile_tmp = new double[tb->ninput]; //temporary arrays + double *ffile_tmp = new double[tb->ninput]; //used for sorting + double *efile_tmp = new double[tb->ninput]; // After re-imaging, does the range of angles cross the 0 or 2*PI boundary? // If so, find the discontinuity: diff --git a/src/MOLECULE/angle_table.cpp b/src/MOLECULE/angle_table.cpp index 84a1dafb27..2d7356db01 100644 --- a/src/MOLECULE/angle_table.cpp +++ b/src/MOLECULE/angle_table.cpp @@ -618,8 +618,7 @@ double AngleTable::splint(double *xa, double *ya, double *y2a, int n, double x) h = xa[khi]-xa[klo]; a = (xa[khi]-x) / h; b = (x-xa[klo]) / h; - y = a*ya[klo] + b*ya[khi] + - ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0; + y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0; return y; } @@ -635,8 +634,9 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f) double fraction,a,b; const Table *tb = &tables[tabindex[type]]; - int itable = static_cast (x * tb->invdelta); + // invdelta is based on tablength-1 + int itable = static_cast (x * tb->invdelta); if (itable < 0) itable = 0; if (itable >= tablength) itable = tablength-1; @@ -650,11 +650,9 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f) b = (x - tb->ang[itable]) * tb->invdelta; a = 1.0 - b; u = a * tb->e[itable] + b * tb->e[itable+1] + - ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * - tb->deltasq6; + ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6; f = a * tb->f[itable] + b * tb->f[itable+1] + - ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) * - tb->deltasq6; + ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) * tb->deltasq6; } } @@ -684,7 +682,6 @@ void AngleTable::u_lookup(int type, double x, double &u) b = (x - tb->ang[itable]) * tb->invdelta; a = 1.0 - b; u = a * tb->e[itable] + b * tb->e[itable+1] + - ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * - tb->deltasq6; + ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6; } } diff --git a/src/MOLECULE/bond_table.cpp b/src/MOLECULE/bond_table.cpp index 58c07b5aea..48d7377682 100644 --- a/src/MOLECULE/bond_table.cpp +++ b/src/MOLECULE/bond_table.cpp @@ -581,8 +581,7 @@ double BondTable::splint(double *xa, double *ya, double *y2a, int n, double x) h = xa[khi]-xa[klo]; a = (xa[khi]-x) / h; b = (x-xa[klo]) / h; - y = a*ya[klo] + b*ya[khi] + - ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0; + y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0; return y; } @@ -601,11 +600,9 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f) const Table *tb = &tables[tabindex[type]]; const int itable = static_cast ((x - tb->lo) * tb->invdelta); if (itable < 0) - error->one(FLERR,"Bond length < table inner cutoff: " - "type {} length {:.8}",type,x); + error->one(FLERR,"Bond length < table inner cutoff: type {} length {:.8}",type,x); else if (itable >= tablength) - error->one(FLERR,"Bond length > table outer cutoff: " - "type {} length {:.8}",type,x); + error->one(FLERR,"Bond length > table outer cutoff: type {} length {:.8}",type,x); if (tabstyle == LINEAR) { fraction = (x - tb->r[itable]) * tb->invdelta; @@ -617,10 +614,8 @@ void BondTable::uf_lookup(int type, double x, double &u, double &f) b = (x - tb->r[itable]) * tb->invdelta; a = 1.0 - b; u = a * tb->e[itable] + b * tb->e[itable+1] + - ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * - tb->deltasq6; + ((a*a*a-a)*tb->e2[itable] + (b*b*b-b)*tb->e2[itable+1]) * tb->deltasq6; f = a * tb->f[itable] + b * tb->f[itable+1] + - ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) * - tb->deltasq6; + ((a*a*a-a)*tb->f2[itable] + (b*b*b-b)*tb->f2[itable+1]) * tb->deltasq6; } } diff --git a/src/MOLECULE/dihedral_table.cpp b/src/MOLECULE/dihedral_table.cpp index a91324dd98..dbca4a85c1 100644 --- a/src/MOLECULE/dihedral_table.cpp +++ b/src/MOLECULE/dihedral_table.cpp @@ -189,11 +189,8 @@ static int solve_cyc_tridiag( const double diag[], size_t d_stride, spline and splint routines modified from Numerical Recipes ------------------------------------------------------------------------- */ -static int cyc_spline(double const *xa, - double const *ya, - int n, - double period, - double *y2a, bool warn) +static int cyc_spline(double const *xa, double const *ya, int n, + double period, double *y2a, bool warn) { double *diag = new double[n]; @@ -264,12 +261,8 @@ static int cyc_spline(double const *xa, // range should not exceed period (ie xa[n-1] < xa[0] + period). // x must lie in the range: [(xa[n-1]-period), (xa[0]+period)] // "period" is typically 2*PI. -static double cyc_splint(double const *xa, - double const *ya, - double const *y2a, - int n, - double period, - double x) +static double cyc_splint(double const *xa, double const *ya, double const *y2a, + int n, double period, double x) { int klo = -1; int khi = n; @@ -302,11 +295,8 @@ static double cyc_splint(double const *xa, } // cyc_splint() -static double cyc_lin(double const *xa, - double const *ya, - int n, - double period, - double x) +static double cyc_lin(double const *xa, double const *ya, + int n, double period, double x) { int klo = -1; int khi = n; @@ -337,21 +327,14 @@ static double cyc_lin(double const *xa, } // cyc_lin() - - - // cyc_splintD(): Evaluate the deriviative of a cyclic spline at position x, // with n control points at xa[], ya[], with parameters y2a[]. // The xa[] must be monotonically increasing and their // range should not exceed period (ie xa[n-1] < xa[0] + period). // x must lie in the range: [(xa[n-1]-period), (xa[0]+period)] // "period" is typically 2*PI. -static double cyc_splintD(double const *xa, - double const *ya, - double const *y2a, - int n, - double period, - double x) +static double cyc_splintD(double const *xa, double const *ya, double const *y2a, + int n, double period, double x) { int klo = -1; int khi = n; // (not n-1) @@ -829,9 +812,9 @@ void DihedralTable::coeff(int narg, char **arg) // We also want the angles to be sorted in increasing order. // This messy code fixes these problems with the user's data: { - double *phifile_tmp = new double [tb->ninput]; //temporary arrays - double *ffile_tmp = new double [tb->ninput]; //used for sorting - double *efile_tmp = new double [tb->ninput]; + double *phifile_tmp = new double[tb->ninput]; //temporary arrays + double *ffile_tmp = new double[tb->ninput]; //used for sorting + double *efile_tmp = new double[tb->ninput]; // After re-imaging, does the range of angles cross the 0 or 2*PI boundary? // If so, find the discontinuity: @@ -1184,8 +1167,7 @@ void DihedralTable::compute_table(Table *tb) if (! tb->f_unspecified) tb->f[i] = cyc_splint(tb->phifile,tb->ffile,tb->f2file,tb->ninput,MY_2PI,phi); } - } // if (tabstyle == SPLINE) - else if (tabstyle == LINEAR) { + } else if (tabstyle == LINEAR) { if (! tb->f_unspecified) { for (int i = 0; i < tablength; i++) { double phi = i*tb->delta; @@ -1193,8 +1175,7 @@ void DihedralTable::compute_table(Table *tb) tb->e[i]= cyc_lin(tb->phifile,tb->efile,tb->ninput,MY_2PI,phi); tb->f[i]= cyc_lin(tb->phifile,tb->ffile,tb->ninput,MY_2PI,phi); } - } - else { + } else { for (int i = 0; i < tablength; i++) { double phi = i*tb->delta; tb->phi[i] = phi; @@ -1269,8 +1250,7 @@ void DihedralTable::param_extract(Table *tb, char *line) //else if (word == "EQ") { // tb->theta0 = values.next_double(); //} - else error->one(FLERR,"Invalid keyword in dihedral angle " - "table parameters ({})", word); + else error->one(FLERR,"Invalid keyword in dihedral angle table parameters ({})", word); } } catch (TokenizerException &e) { error->one(FLERR, e.what()); diff --git a/src/OPENMP/angle_table_omp.cpp b/src/OPENMP/angle_table_omp.cpp index 892f9295a5..cca34a67f7 100644 --- a/src/OPENMP/angle_table_omp.cpp +++ b/src/OPENMP/angle_table_omp.cpp @@ -16,15 +16,16 @@ Contributing author: Axel Kohlmeyer (Temple U) ------------------------------------------------------------------------- */ -#include "omp_compat.h" #include "angle_table_omp.h" -#include + #include "atom.h" #include "comm.h" #include "force.h" #include "neighbor.h" +#include +#include "omp_compat.h" #include "suffix.h" using namespace LAMMPS_NS; diff --git a/src/OPENMP/bond_table_omp.cpp b/src/OPENMP/bond_table_omp.cpp index faadca456a..dcc13c85c9 100644 --- a/src/OPENMP/bond_table_omp.cpp +++ b/src/OPENMP/bond_table_omp.cpp @@ -16,16 +16,16 @@ Contributing author: Axel Kohlmeyer (Temple U) ------------------------------------------------------------------------- */ -#include "omp_compat.h" #include "bond_table_omp.h" + #include "atom.h" #include "comm.h" #include "force.h" #include "neighbor.h" - #include +#include "omp_compat.h" #include "suffix.h" using namespace LAMMPS_NS;