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

cosmetic and whitespace changes

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-01-04 23:08:32 -05:00
parent 40abc0886c
commit 4aba9e9bb6
6 changed files with 38 additions and 74 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
src/EXTRA-MOLECULE/dihedral_table_cut.cpp
View File

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

15
src/MOLECULE/angle_table.cpp
View File

@ -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<int> (x * tb->invdelta);
// invdelta is based on tablength-1
int itable = static_cast<int> (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;
}
}

15
src/MOLECULE/bond_table.cpp
View File

@ -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<int> ((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;
}
}

48
src/MOLECULE/dihedral_table.cpp
View File

@ -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());

5
src/OPENMP/angle_table_omp.cpp
View File

@ -16,15 +16,16 @@
Contributing author: Axel Kohlmeyer (Temple U)
------------------------------------------------------------------------- */
#include "omp_compat.h"
#include "angle_table_omp.h"
#include <cmath>
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include <cmath>
#include "omp_compat.h"
#include "suffix.h"
using namespace LAMMPS_NS;

4
src/OPENMP/bond_table_omp.cpp
View File

@ -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 <cmath>
#include "omp_compat.h"
#include "suffix.h"
using namespace LAMMPS_NS;