enable and apply clang-format
This commit is contained in:
Axel Kohlmeyer
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -18,25 +17,25 @@
|
||||
|
||||
#include "angle_charmm.h"
|
||||
|
||||
#include <cmath>
|
||||
#include "atom.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
using MathConst::MY_PI;
|
||||
|
||||
#define SMALL 0.001
|
||||
static constexpr double SMALL = 0.001;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleCharmm::AngleCharmm(LAMMPS *lmp) : Angle(lmp) {}
|
||||
AngleCharmm::AngleCharmm(LAMMPS *_lmp) : Angle(_lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -55,15 +54,15 @@ AngleCharmm::~AngleCharmm()
|
||||
|
||||
void AngleCharmm::compute(int eflag, int vflag)
|
||||
{
|
||||
int i1,i2,i3,n,type;
|
||||
double delx1,dely1,delz1,delx2,dely2,delz2;
|
||||
double eangle,f1[3],f3[3];
|
||||
double dtheta,tk;
|
||||
double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
|
||||
double delxUB,delyUB,delzUB,rsqUB,rUB,dr,rk,forceUB;
|
||||
int i1, i2, i3, n, type;
|
||||
double delx1, dely1, delz1, delx2, dely2, delz2;
|
||||
double eangle, f1[3], f3[3];
|
||||
double dtheta, tk;
|
||||
double rsq1, rsq2, r1, r2, c, s, a, a11, a12, a22;
|
||||
double delxUB, delyUB, delzUB, rsqUB, rUB, dr, rk, forceUB;
|
||||
|
||||
eangle = 0.0;
|
||||
ev_init(eflag,vflag);
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
@ -84,7 +83,7 @@ void AngleCharmm::compute(int eflag, int vflag)
|
||||
dely1 = x[i1][1] - x[i2][1];
|
||||
delz1 = x[i1][2] - x[i2][2];
|
||||
|
||||
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
|
||||
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
|
||||
r1 = sqrt(rsq1);
|
||||
|
||||
// 2nd bond
|
||||
@ -93,7 +92,7 @@ void AngleCharmm::compute(int eflag, int vflag)
|
||||
dely2 = x[i3][1] - x[i2][1];
|
||||
delz2 = x[i3][2] - x[i2][2];
|
||||
|
||||
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
|
||||
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
|
||||
r2 = sqrt(rsq2);
|
||||
|
||||
// Urey-Bradley bond
|
||||
@ -102,7 +101,7 @@ void AngleCharmm::compute(int eflag, int vflag)
|
||||
delyUB = x[i3][1] - x[i1][1];
|
||||
delzUB = x[i3][2] - x[i1][2];
|
||||
|
||||
rsqUB = delxUB*delxUB + delyUB*delyUB + delzUB*delzUB;
|
||||
rsqUB = delxUB * delxUB + delyUB * delyUB + delzUB * delzUB;
|
||||
rUB = sqrt(rsqUB);
|
||||
|
||||
// Urey-Bradley force & energy
|
||||
@ -110,42 +109,44 @@ void AngleCharmm::compute(int eflag, int vflag)
|
||||
dr = rUB - r_ub[type];
|
||||
rk = k_ub[type] * dr;
|
||||
|
||||
if (rUB > 0.0) forceUB = -2.0*rk/rUB;
|
||||
else forceUB = 0.0;
|
||||
if (rUB > 0.0)
|
||||
forceUB = -2.0 * rk / rUB;
|
||||
else
|
||||
forceUB = 0.0;
|
||||
|
||||
if (eflag) eangle = rk*dr;
|
||||
if (eflag) eangle = rk * dr;
|
||||
|
||||
// angle (cos and sin)
|
||||
|
||||
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
s = sqrt(1.0 - c*c);
|
||||
s = sqrt(1.0 - c * c);
|
||||
if (s < SMALL) s = SMALL;
|
||||
s = 1.0/s;
|
||||
s = 1.0 / s;
|
||||
|
||||
// harmonic force & energy
|
||||
|
||||
dtheta = acos(c) - theta0[type];
|
||||
tk = k[type] * dtheta;
|
||||
|
||||
if (eflag) eangle += tk*dtheta;
|
||||
if (eflag) eangle += tk * dtheta;
|
||||
|
||||
a = -2.0 * tk * s;
|
||||
a11 = a*c / rsq1;
|
||||
a12 = -a / (r1*r2);
|
||||
a22 = a*c / rsq2;
|
||||
a11 = a * c / rsq1;
|
||||
a12 = -a / (r1 * r2);
|
||||
a22 = a * c / rsq2;
|
||||
|
||||
f1[0] = a11*delx1 + a12*delx2 - delxUB*forceUB;
|
||||
f1[1] = a11*dely1 + a12*dely2 - delyUB*forceUB;
|
||||
f1[2] = a11*delz1 + a12*delz2 - delzUB*forceUB;
|
||||
f1[0] = a11 * delx1 + a12 * delx2 - delxUB * forceUB;
|
||||
f1[1] = a11 * dely1 + a12 * dely2 - delyUB * forceUB;
|
||||
f1[2] = a11 * delz1 + a12 * delz2 - delzUB * forceUB;
|
||||
|
||||
f3[0] = a22*delx2 + a12*delx1 + delxUB*forceUB;
|
||||
f3[1] = a22*dely2 + a12*dely1 + delyUB*forceUB;
|
||||
f3[2] = a22*delz2 + a12*delz1 + delzUB*forceUB;
|
||||
f3[0] = a22 * delx2 + a12 * delx1 + delxUB * forceUB;
|
||||
f3[1] = a22 * dely2 + a12 * dely1 + delyUB * forceUB;
|
||||
f3[2] = a22 * delz2 + a12 * delz1 + delzUB * forceUB;
|
||||
|
||||
// apply force to each of 3 atoms
|
||||
|
||||
@ -167,8 +168,9 @@ void AngleCharmm::compute(int eflag, int vflag)
|
||||
f[i3][2] += f3[2];
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
|
||||
delx1,dely1,delz1,delx2,dely2,delz2);
|
||||
if (evflag)
|
||||
ev_tally(i1, i2, i3, nlocal, newton_bond, eangle, f1, f3, delx1, dely1, delz1, delx2, dely2,
|
||||
delz2);
|
||||
}
|
||||
}
|
||||
|
||||
@ -177,14 +179,14 @@ void AngleCharmm::compute(int eflag, int vflag)
|
||||
void AngleCharmm::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nangletypes;
|
||||
const int np1 = atom->nangletypes + 1;
|
||||
|
||||
memory->create(k,n+1,"angle:k");
|
||||
memory->create(theta0,n+1,"angle:theta0");
|
||||
memory->create(k_ub,n+1,"angle:k_ub");
|
||||
memory->create(r_ub,n+1,"angle:r_ub");
|
||||
memory->create(setflag,n+1,"angle:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
memory->create(k, np1, "angle:k");
|
||||
memory->create(theta0, np1, "angle:theta0");
|
||||
memory->create(k_ub, np1, "angle:k_ub");
|
||||
memory->create(r_ub, np1, "angle:r_ub");
|
||||
memory->create(setflag, np1, "angle:setflag");
|
||||
for (int i = 1; i < np1; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -193,30 +195,30 @@ void AngleCharmm::allocate()
|
||||
|
||||
void AngleCharmm::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 5) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (narg != 5) error->all(FLERR, "Incorrect args for angle coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi;
|
||||
utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
|
||||
int ilo, ihi;
|
||||
utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
|
||||
|
||||
double k_one = utils::numeric(FLERR,arg[1],false,lmp);
|
||||
double theta0_one = utils::numeric(FLERR,arg[2],false,lmp);
|
||||
double k_ub_one = utils::numeric(FLERR,arg[3],false,lmp);
|
||||
double r_ub_one = utils::numeric(FLERR,arg[4],false,lmp);
|
||||
double k_one = utils::numeric(FLERR, arg[1], false, lmp);
|
||||
double theta0_one = utils::numeric(FLERR, arg[2], false, lmp);
|
||||
double k_ub_one = utils::numeric(FLERR, arg[3], false, lmp);
|
||||
double r_ub_one = utils::numeric(FLERR, arg[4], false, lmp);
|
||||
|
||||
// convert theta0 from degrees to radians
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
k[i] = k_one;
|
||||
theta0[i] = theta0_one/180.0 * MY_PI;
|
||||
theta0[i] = theta0_one / 180.0 * MY_PI;
|
||||
k_ub[i] = k_ub_one;
|
||||
r_ub[i] = r_ub_one;
|
||||
setflag[i] = 1;
|
||||
count++;
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (count == 0) error->all(FLERR, "Incorrect args for angle coefficients");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -232,10 +234,10 @@ double AngleCharmm::equilibrium_angle(int i)
|
||||
|
||||
void AngleCharmm::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k_ub[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&r_ub[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k[1], sizeof(double), atom->nangletypes, fp);
|
||||
fwrite(&theta0[1], sizeof(double), atom->nangletypes, fp);
|
||||
fwrite(&k_ub[1], sizeof(double), atom->nangletypes, fp);
|
||||
fwrite(&r_ub[1], sizeof(double), atom->nangletypes, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -247,15 +249,15 @@ void AngleCharmm::read_restart(FILE *fp)
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) {
|
||||
utils::sfread(FLERR,&k[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR,&theta0[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR,&k_ub[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR,&r_ub[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR, &k[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &theta0[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &k_ub[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &r_ub[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&k_ub[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&r_ub[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&k[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&theta0[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&k_ub[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&r_ub[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
|
||||
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
|
||||
}
|
||||
@ -267,8 +269,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], theta0[i] / MY_PI * 180.0, k_ub[i], r_ub[i]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -280,23 +281,23 @@ double AngleCharmm::single(int type, int i1, int i2, int i3)
|
||||
double delx1 = x[i1][0] - x[i2][0];
|
||||
double dely1 = x[i1][1] - x[i2][1];
|
||||
double delz1 = x[i1][2] - x[i2][2];
|
||||
domain->minimum_image(delx1,dely1,delz1);
|
||||
double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
|
||||
domain->minimum_image(delx1, dely1, delz1);
|
||||
double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1);
|
||||
|
||||
double delx2 = x[i3][0] - x[i2][0];
|
||||
double dely2 = x[i3][1] - x[i2][1];
|
||||
double delz2 = x[i3][2] - x[i2][2];
|
||||
domain->minimum_image(delx2,dely2,delz2);
|
||||
double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
|
||||
domain->minimum_image(delx2, dely2, delz2);
|
||||
double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2);
|
||||
|
||||
double delxUB = x[i3][0] - x[i1][0];
|
||||
double delyUB = x[i3][1] - x[i1][1];
|
||||
double delzUB = x[i3][2] - x[i1][2];
|
||||
domain->minimum_image(delxUB,delyUB,delzUB);
|
||||
double rUB = sqrt(delxUB*delxUB + delyUB*delyUB + delzUB*delzUB);
|
||||
domain->minimum_image(delxUB, delyUB, delzUB);
|
||||
double rUB = sqrt(delxUB * delxUB + delyUB * delyUB + delzUB * delzUB);
|
||||
|
||||
double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
@ -305,5 +306,5 @@ double AngleCharmm::single(int type, int i1, int i2, int i3)
|
||||
double dr = rUB - r_ub[type];
|
||||
double rk = k_ub[type] * dr;
|
||||
|
||||
return (tk*dtheta + rk*dr);
|
||||
return (tk * dtheta + rk * dr);
|
||||
}
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -14,25 +13,23 @@
|
||||
|
||||
#include "angle_cosine.h"
|
||||
|
||||
#include <cmath>
|
||||
#include "atom.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
#define SMALL 0.001
|
||||
using MathConst::MY_PI;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleCosine::AngleCosine(LAMMPS *lmp) : Angle(lmp) {}
|
||||
AngleCosine::AngleCosine(LAMMPS *_lmp) : Angle(_lmp) {}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
@ -48,13 +45,13 @@ AngleCosine::~AngleCosine()
|
||||
|
||||
void AngleCosine::compute(int eflag, int vflag)
|
||||
{
|
||||
int i1,i2,i3,n,type;
|
||||
double delx1,dely1,delz1,delx2,dely2,delz2;
|
||||
double eangle,f1[3],f3[3];
|
||||
double rsq1,rsq2,r1,r2,c,a,a11,a12,a22;
|
||||
int i1, i2, i3, n, type;
|
||||
double delx1, dely1, delz1, delx2, dely2, delz2;
|
||||
double eangle, f1[3], f3[3];
|
||||
double rsq1, rsq2, r1, r2, c, a, a11, a12, a22;
|
||||
|
||||
eangle = 0.0;
|
||||
ev_init(eflag,vflag);
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
@ -75,7 +72,7 @@ void AngleCosine::compute(int eflag, int vflag)
|
||||
dely1 = x[i1][1] - x[i2][1];
|
||||
delz1 = x[i1][2] - x[i2][2];
|
||||
|
||||
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
|
||||
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
|
||||
r1 = sqrt(rsq1);
|
||||
|
||||
// 2nd bond
|
||||
@ -84,31 +81,31 @@ void AngleCosine::compute(int eflag, int vflag)
|
||||
dely2 = x[i3][1] - x[i2][1];
|
||||
delz2 = x[i3][2] - x[i2][2];
|
||||
|
||||
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
|
||||
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
|
||||
r2 = sqrt(rsq2);
|
||||
|
||||
// c = cosine of angle
|
||||
|
||||
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
// force & energy
|
||||
|
||||
if (eflag) eangle = k[type]*(1.0+c);
|
||||
if (eflag) eangle = k[type] * (1.0 + c);
|
||||
|
||||
a = k[type];
|
||||
a11 = a*c / rsq1;
|
||||
a12 = -a / (r1*r2);
|
||||
a22 = a*c / rsq2;
|
||||
a11 = a * c / rsq1;
|
||||
a12 = -a / (r1 * r2);
|
||||
a22 = a * c / rsq2;
|
||||
|
||||
f1[0] = a11*delx1 + a12*delx2;
|
||||
f1[1] = a11*dely1 + a12*dely2;
|
||||
f1[2] = a11*delz1 + a12*delz2;
|
||||
f3[0] = a22*delx2 + a12*delx1;
|
||||
f3[1] = a22*dely2 + a12*dely1;
|
||||
f3[2] = a22*delz2 + a12*delz1;
|
||||
f1[0] = a11 * delx1 + a12 * delx2;
|
||||
f1[1] = a11 * dely1 + a12 * dely2;
|
||||
f1[2] = a11 * delz1 + a12 * delz2;
|
||||
f3[0] = a22 * delx2 + a12 * delx1;
|
||||
f3[1] = a22 * dely2 + a12 * dely1;
|
||||
f3[2] = a22 * delz2 + a12 * delz1;
|
||||
|
||||
// apply force to each of 3 atoms
|
||||
|
||||
@ -130,8 +127,9 @@ void AngleCosine::compute(int eflag, int vflag)
|
||||
f[i3][2] += f3[2];
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
|
||||
delx1,dely1,delz1,delx2,dely2,delz2);
|
||||
if (evflag)
|
||||
ev_tally(i1, i2, i3, nlocal, newton_bond, eangle, f1, f3, delx1, dely1, delz1, delx2, dely2,
|
||||
delz2);
|
||||
}
|
||||
}
|
||||
|
||||
@ -140,11 +138,11 @@ void AngleCosine::compute(int eflag, int vflag)
|
||||
void AngleCosine::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nangletypes;
|
||||
const int np1 = atom->nangletypes + 1;
|
||||
|
||||
memory->create(k,n+1,"angle:k");
|
||||
memory->create(setflag,n+1,"angle:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
memory->create(k, np1, "angle:k");
|
||||
memory->create(setflag, np1, "angle:setflag");
|
||||
for (int i = 1; i < np1; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -153,13 +151,13 @@ void AngleCosine::allocate()
|
||||
|
||||
void AngleCosine::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 2) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (narg != 2) error->all(FLERR, "Incorrect args for angle coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi;
|
||||
utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
|
||||
int ilo, ihi;
|
||||
utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
|
||||
|
||||
double k_one = utils::numeric(FLERR,arg[1],false,lmp);
|
||||
double k_one = utils::numeric(FLERR, arg[1], false, lmp);
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
@ -168,7 +166,7 @@ void AngleCosine::coeff(int narg, char **arg)
|
||||
count++;
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (count == 0) error->all(FLERR, "Incorrect args for angle coefficients");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -184,7 +182,7 @@ double AngleCosine::equilibrium_angle(int /*i*/)
|
||||
|
||||
void AngleCosine::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k[1], sizeof(double), atom->nangletypes, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -195,8 +193,9 @@ void AngleCosine::read_restart(FILE *fp)
|
||||
{
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) utils::sfread(FLERR,&k[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
if (comm->me == 0)
|
||||
utils::sfread(FLERR, &k[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
MPI_Bcast(&k[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
|
||||
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
|
||||
}
|
||||
@ -207,8 +206,7 @@ void AngleCosine::read_restart(FILE *fp)
|
||||
|
||||
void AngleCosine::write_data(FILE *fp)
|
||||
{
|
||||
for (int i = 1; i <= atom->nangletypes; i++)
|
||||
fprintf(fp,"%d %g\n",i,k[i]);
|
||||
for (int i = 1; i <= atom->nangletypes; i++) fprintf(fp, "%d %g\n", i, k[i]);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -220,19 +218,19 @@ double AngleCosine::single(int type, int i1, int i2, int i3)
|
||||
double delx1 = x[i1][0] - x[i2][0];
|
||||
double dely1 = x[i1][1] - x[i2][1];
|
||||
double delz1 = x[i1][2] - x[i2][2];
|
||||
domain->minimum_image(delx1,dely1,delz1);
|
||||
double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
|
||||
domain->minimum_image(delx1, dely1, delz1);
|
||||
double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1);
|
||||
|
||||
double delx2 = x[i3][0] - x[i2][0];
|
||||
double dely2 = x[i3][1] - x[i2][1];
|
||||
double delz2 = x[i3][2] - x[i2][2];
|
||||
domain->minimum_image(delx2,dely2,delz2);
|
||||
double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
|
||||
domain->minimum_image(delx2, dely2, delz2);
|
||||
double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2);
|
||||
|
||||
double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
return k[type]*(1.0+c);
|
||||
return k[type] * (1.0 + c);
|
||||
}
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -18,25 +17,23 @@
|
||||
|
||||
#include "angle_cosine_squared.h"
|
||||
|
||||
#include <cmath>
|
||||
#include "atom.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
#define SMALL 0.001
|
||||
using MathConst::MY_PI;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleCosineSquared::AngleCosineSquared(LAMMPS *lmp) : Angle(lmp)
|
||||
AngleCosineSquared::AngleCosineSquared(LAMMPS *_lmp) : Angle(_lmp)
|
||||
{
|
||||
k = nullptr;
|
||||
theta0 = nullptr;
|
||||
@ -57,14 +54,14 @@ AngleCosineSquared::~AngleCosineSquared()
|
||||
|
||||
void AngleCosineSquared::compute(int eflag, int vflag)
|
||||
{
|
||||
int i1,i2,i3,n,type;
|
||||
double delx1,dely1,delz1,delx2,dely2,delz2;
|
||||
double eangle,f1[3],f3[3];
|
||||
double dcostheta,tk;
|
||||
double rsq1,rsq2,r1,r2,c,a,a11,a12,a22;
|
||||
int i1, i2, i3, n, type;
|
||||
double delx1, dely1, delz1, delx2, dely2, delz2;
|
||||
double eangle, f1[3], f3[3];
|
||||
double dcostheta, tk;
|
||||
double rsq1, rsq2, r1, r2, c, a, a11, a12, a22;
|
||||
|
||||
eangle = 0.0;
|
||||
ev_init(eflag,vflag);
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
@ -85,7 +82,7 @@ void AngleCosineSquared::compute(int eflag, int vflag)
|
||||
dely1 = x[i1][1] - x[i2][1];
|
||||
delz1 = x[i1][2] - x[i2][2];
|
||||
|
||||
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
|
||||
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
|
||||
r1 = sqrt(rsq1);
|
||||
|
||||
// 2nd bond
|
||||
@ -94,13 +91,13 @@ void AngleCosineSquared::compute(int eflag, int vflag)
|
||||
dely2 = x[i3][1] - x[i2][1];
|
||||
delz2 = x[i3][2] - x[i2][2];
|
||||
|
||||
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
|
||||
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
|
||||
r2 = sqrt(rsq2);
|
||||
|
||||
// angle (cos and sin)
|
||||
|
||||
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
@ -110,19 +107,19 @@ void AngleCosineSquared::compute(int eflag, int vflag)
|
||||
dcostheta = c - cos(theta0[type]);
|
||||
tk = k[type] * dcostheta;
|
||||
|
||||
if (eflag) eangle = tk*dcostheta;
|
||||
if (eflag) eangle = tk * dcostheta;
|
||||
|
||||
a = 2.0 * tk;
|
||||
a11 = a*c / rsq1;
|
||||
a12 = -a / (r1*r2);
|
||||
a22 = a*c / rsq2;
|
||||
a11 = a * c / rsq1;
|
||||
a12 = -a / (r1 * r2);
|
||||
a22 = a * c / rsq2;
|
||||
|
||||
f1[0] = a11*delx1 + a12*delx2;
|
||||
f1[1] = a11*dely1 + a12*dely2;
|
||||
f1[2] = a11*delz1 + a12*delz2;
|
||||
f3[0] = a22*delx2 + a12*delx1;
|
||||
f3[1] = a22*dely2 + a12*dely1;
|
||||
f3[2] = a22*delz2 + a12*delz1;
|
||||
f1[0] = a11 * delx1 + a12 * delx2;
|
||||
f1[1] = a11 * dely1 + a12 * dely2;
|
||||
f1[2] = a11 * delz1 + a12 * delz2;
|
||||
f3[0] = a22 * delx2 + a12 * delx1;
|
||||
f3[1] = a22 * dely2 + a12 * dely1;
|
||||
f3[2] = a22 * delz2 + a12 * delz1;
|
||||
|
||||
// apply force to each of 3 atoms
|
||||
|
||||
@ -144,8 +141,9 @@ void AngleCosineSquared::compute(int eflag, int vflag)
|
||||
f[i3][2] += f3[2];
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
|
||||
delx1,dely1,delz1,delx2,dely2,delz2);
|
||||
if (evflag)
|
||||
ev_tally(i1, i2, i3, nlocal, newton_bond, eangle, f1, f3, delx1, dely1, delz1, delx2, dely2,
|
||||
delz2);
|
||||
}
|
||||
}
|
||||
|
||||
@ -154,13 +152,13 @@ void AngleCosineSquared::compute(int eflag, int vflag)
|
||||
void AngleCosineSquared::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nangletypes;
|
||||
const int np1 = atom->nangletypes + 1;
|
||||
|
||||
memory->create(k,n+1,"angle:k");
|
||||
memory->create(theta0,n+1,"angle:theta0");
|
||||
memory->create(k, np1, "angle:k");
|
||||
memory->create(theta0, np1, "angle:theta0");
|
||||
|
||||
memory->create(setflag,n+1,"angle:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
memory->create(setflag, np1, "angle:setflag");
|
||||
for (int i = 1; i < np1; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -169,26 +167,26 @@ void AngleCosineSquared::allocate()
|
||||
|
||||
void AngleCosineSquared::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 3) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (narg != 3) error->all(FLERR, "Incorrect args for angle coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi;
|
||||
utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
|
||||
int ilo, ihi;
|
||||
utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
|
||||
|
||||
double k_one = utils::numeric(FLERR,arg[1],false,lmp);
|
||||
double theta0_one = utils::numeric(FLERR,arg[2],false,lmp);
|
||||
double k_one = utils::numeric(FLERR, arg[1], false, lmp);
|
||||
double theta0_one = utils::numeric(FLERR, arg[2], false, lmp);
|
||||
|
||||
// convert theta0 from degrees to radians
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
k[i] = k_one;
|
||||
theta0[i] = theta0_one/180.0 * MY_PI;
|
||||
theta0[i] = theta0_one / 180.0 * MY_PI;
|
||||
setflag[i] = 1;
|
||||
count++;
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (count == 0) error->all(FLERR, "Incorrect args for angle coefficients");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -204,8 +202,8 @@ double AngleCosineSquared::equilibrium_angle(int i)
|
||||
|
||||
void AngleCosineSquared::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k[1], sizeof(double), atom->nangletypes, fp);
|
||||
fwrite(&theta0[1], sizeof(double), atom->nangletypes, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -217,11 +215,11 @@ void AngleCosineSquared::read_restart(FILE *fp)
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) {
|
||||
utils::sfread(FLERR,&k[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR,&theta0[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR, &k[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &theta0[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&k[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&theta0[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
|
||||
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
|
||||
}
|
||||
@ -233,7 +231,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], theta0[i] / MY_PI * 180.0);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -245,21 +243,21 @@ double AngleCosineSquared::single(int type, int i1, int i2, int i3)
|
||||
double delx1 = x[i1][0] - x[i2][0];
|
||||
double dely1 = x[i1][1] - x[i2][1];
|
||||
double delz1 = x[i1][2] - x[i2][2];
|
||||
domain->minimum_image(delx1,dely1,delz1);
|
||||
double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
|
||||
domain->minimum_image(delx1, dely1, delz1);
|
||||
double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1);
|
||||
|
||||
double delx2 = x[i3][0] - x[i2][0];
|
||||
double dely2 = x[i3][1] - x[i2][1];
|
||||
double delz2 = x[i3][2] - x[i2][2];
|
||||
domain->minimum_image(delx2,dely2,delz2);
|
||||
double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
|
||||
domain->minimum_image(delx2, dely2, delz2);
|
||||
double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2);
|
||||
|
||||
double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
double dcostheta = c - cos(theta0[type]);
|
||||
double tk = k[type] * dcostheta;
|
||||
return tk*dcostheta;
|
||||
return tk * dcostheta;
|
||||
}
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -14,25 +13,25 @@
|
||||
|
||||
#include "angle_harmonic.h"
|
||||
|
||||
#include <cmath>
|
||||
#include "atom.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "neighbor.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
|
||||
#define SMALL 0.001
|
||||
static constexpr double SMALL = 0.001;
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleHarmonic::AngleHarmonic(LAMMPS *lmp) : Angle(lmp)
|
||||
AngleHarmonic::AngleHarmonic(LAMMPS *_lmp) : Angle(_lmp)
|
||||
{
|
||||
k = nullptr;
|
||||
theta0 = nullptr;
|
||||
@ -53,14 +52,14 @@ AngleHarmonic::~AngleHarmonic()
|
||||
|
||||
void AngleHarmonic::compute(int eflag, int vflag)
|
||||
{
|
||||
int i1,i2,i3,n,type;
|
||||
double delx1,dely1,delz1,delx2,dely2,delz2;
|
||||
double eangle,f1[3],f3[3];
|
||||
double dtheta,tk;
|
||||
double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
|
||||
int i1, i2, i3, n, type;
|
||||
double delx1, dely1, delz1, delx2, dely2, delz2;
|
||||
double eangle, f1[3], f3[3];
|
||||
double dtheta, tk;
|
||||
double rsq1, rsq2, r1, r2, c, s, a, a11, a12, a22;
|
||||
|
||||
eangle = 0.0;
|
||||
ev_init(eflag,vflag);
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
@ -81,7 +80,7 @@ void AngleHarmonic::compute(int eflag, int vflag)
|
||||
dely1 = x[i1][1] - x[i2][1];
|
||||
delz1 = x[i1][2] - x[i2][2];
|
||||
|
||||
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
|
||||
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
|
||||
r1 = sqrt(rsq1);
|
||||
|
||||
// 2nd bond
|
||||
@ -90,39 +89,39 @@ void AngleHarmonic::compute(int eflag, int vflag)
|
||||
dely2 = x[i3][1] - x[i2][1];
|
||||
delz2 = x[i3][2] - x[i2][2];
|
||||
|
||||
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
|
||||
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
|
||||
r2 = sqrt(rsq2);
|
||||
|
||||
// angle (cos and sin)
|
||||
|
||||
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
s = sqrt(1.0 - c*c);
|
||||
s = sqrt(1.0 - c * c);
|
||||
if (s < SMALL) s = SMALL;
|
||||
s = 1.0/s;
|
||||
s = 1.0 / s;
|
||||
|
||||
// force & energy
|
||||
|
||||
dtheta = acos(c) - theta0[type];
|
||||
tk = k[type] * dtheta;
|
||||
|
||||
if (eflag) eangle = tk*dtheta;
|
||||
if (eflag) eangle = tk * dtheta;
|
||||
|
||||
a = -2.0 * tk * s;
|
||||
a11 = a*c / rsq1;
|
||||
a12 = -a / (r1*r2);
|
||||
a22 = a*c / rsq2;
|
||||
a11 = a * c / rsq1;
|
||||
a12 = -a / (r1 * r2);
|
||||
a22 = a * c / rsq2;
|
||||
|
||||
f1[0] = a11*delx1 + a12*delx2;
|
||||
f1[1] = a11*dely1 + a12*dely2;
|
||||
f1[2] = a11*delz1 + a12*delz2;
|
||||
f3[0] = a22*delx2 + a12*delx1;
|
||||
f3[1] = a22*dely2 + a12*dely1;
|
||||
f3[2] = a22*delz2 + a12*delz1;
|
||||
f1[0] = a11 * delx1 + a12 * delx2;
|
||||
f1[1] = a11 * dely1 + a12 * dely2;
|
||||
f1[2] = a11 * delz1 + a12 * delz2;
|
||||
f3[0] = a22 * delx2 + a12 * delx1;
|
||||
f3[1] = a22 * dely2 + a12 * dely1;
|
||||
f3[2] = a22 * delz2 + a12 * delz1;
|
||||
|
||||
// apply force to each of 3 atoms
|
||||
|
||||
@ -144,8 +143,9 @@ void AngleHarmonic::compute(int eflag, int vflag)
|
||||
f[i3][2] += f3[2];
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
|
||||
delx1,dely1,delz1,delx2,dely2,delz2);
|
||||
if (evflag)
|
||||
ev_tally(i1, i2, i3, nlocal, newton_bond, eangle, f1, f3, delx1, dely1, delz1, delx2, dely2,
|
||||
delz2);
|
||||
}
|
||||
}
|
||||
|
||||
@ -154,13 +154,13 @@ void AngleHarmonic::compute(int eflag, int vflag)
|
||||
void AngleHarmonic::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nangletypes;
|
||||
const int np1 = atom->nangletypes + 1;
|
||||
|
||||
memory->create(k,n+1,"angle:k");
|
||||
memory->create(theta0,n+1,"angle:theta0");
|
||||
memory->create(k, np1, "angle:k");
|
||||
memory->create(theta0, np1, "angle:theta0");
|
||||
|
||||
memory->create(setflag,n+1,"angle:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
memory->create(setflag, np1, "angle:setflag");
|
||||
for (int i = 1; i < np1; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -169,26 +169,26 @@ void AngleHarmonic::allocate()
|
||||
|
||||
void AngleHarmonic::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 3) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (narg != 3) error->all(FLERR, "Incorrect args for angle coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi;
|
||||
utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
|
||||
int ilo, ihi;
|
||||
utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
|
||||
|
||||
double k_one = utils::numeric(FLERR,arg[1],false,lmp);
|
||||
double theta0_one = utils::numeric(FLERR,arg[2],false,lmp);
|
||||
double k_one = utils::numeric(FLERR, arg[1], false, lmp);
|
||||
double theta0_one = utils::numeric(FLERR, arg[2], false, lmp);
|
||||
|
||||
// convert theta0 from degrees to radians
|
||||
|
||||
int count = 0;
|
||||
for (int i = ilo; i <= ihi; i++) {
|
||||
k[i] = k_one;
|
||||
theta0[i] = theta0_one/180.0 * MY_PI;
|
||||
theta0[i] = theta0_one / 180.0 * MY_PI;
|
||||
setflag[i] = 1;
|
||||
count++;
|
||||
}
|
||||
|
||||
if (count == 0) error->all(FLERR,"Incorrect args for angle coefficients");
|
||||
if (count == 0) error->all(FLERR, "Incorrect args for angle coefficients");
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -204,8 +204,8 @@ double AngleHarmonic::equilibrium_angle(int i)
|
||||
|
||||
void AngleHarmonic::write_restart(FILE *fp)
|
||||
{
|
||||
fwrite(&k[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&theta0[1],sizeof(double),atom->nangletypes,fp);
|
||||
fwrite(&k[1], sizeof(double), atom->nangletypes, fp);
|
||||
fwrite(&theta0[1], sizeof(double), atom->nangletypes, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -217,11 +217,11 @@ void AngleHarmonic::read_restart(FILE *fp)
|
||||
allocate();
|
||||
|
||||
if (comm->me == 0) {
|
||||
utils::sfread(FLERR,&k[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR,&theta0[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
|
||||
utils::sfread(FLERR, &k[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &theta0[1], sizeof(double), atom->nangletypes, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&k[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&theta0[1],atom->nangletypes,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&k[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&theta0[1], atom->nangletypes, MPI_DOUBLE, 0, world);
|
||||
|
||||
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
|
||||
}
|
||||
@ -233,7 +233,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], theta0[i] / MY_PI * 180.0);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -245,21 +245,21 @@ double AngleHarmonic::single(int type, int i1, int i2, int i3)
|
||||
double delx1 = x[i1][0] - x[i2][0];
|
||||
double dely1 = x[i1][1] - x[i2][1];
|
||||
double delz1 = x[i1][2] - x[i2][2];
|
||||
domain->minimum_image(delx1,dely1,delz1);
|
||||
double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
|
||||
domain->minimum_image(delx1, dely1, delz1);
|
||||
double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1);
|
||||
|
||||
double delx2 = x[i3][0] - x[i2][0];
|
||||
double dely2 = x[i3][1] - x[i2][1];
|
||||
double delz2 = x[i3][2] - x[i2][2];
|
||||
domain->minimum_image(delx2,dely2,delz2);
|
||||
double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
|
||||
domain->minimum_image(delx2, dely2, delz2);
|
||||
double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2);
|
||||
|
||||
double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
double dtheta = acos(c) - theta0[type];
|
||||
double tk = k[type] * dtheta;
|
||||
return tk*dtheta;
|
||||
return tk * dtheta;
|
||||
}
|
||||
|
||||
@ -1,4 +1,3 @@
|
||||
// clang-format off
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://www.lammps.org/, Sandia National Laboratories
|
||||
@ -18,33 +17,31 @@
|
||||
|
||||
#include "angle_table.h"
|
||||
|
||||
#include <cmath>
|
||||
|
||||
#include <cstring>
|
||||
#include "atom.h"
|
||||
#include "neighbor.h"
|
||||
#include "domain.h"
|
||||
#include "comm.h"
|
||||
#include "domain.h"
|
||||
#include "error.h"
|
||||
#include "force.h"
|
||||
#include "math_const.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
|
||||
#include "tokenizer.h"
|
||||
#include "neighbor.h"
|
||||
#include "table_file_reader.h"
|
||||
#include "tokenizer.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
using namespace MathConst;
|
||||
using MathConst::MY_PI;
|
||||
|
||||
enum{LINEAR,SPLINE};
|
||||
enum { LINEAR, SPLINE };
|
||||
|
||||
#define SMALL 0.001
|
||||
#define TINY 1.E-10
|
||||
#define TINY 1.E-10
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
AngleTable::AngleTable(LAMMPS *lmp) : Angle(lmp)
|
||||
AngleTable::AngleTable(LAMMPS *_lmp) : Angle(_lmp)
|
||||
{
|
||||
writedata = 0;
|
||||
ntables = 0;
|
||||
@ -69,14 +66,14 @@ AngleTable::~AngleTable()
|
||||
|
||||
void AngleTable::compute(int eflag, int vflag)
|
||||
{
|
||||
int i1,i2,i3,n,type;
|
||||
double eangle,f1[3],f3[3];
|
||||
double delx1,dely1,delz1,delx2,dely2,delz2;
|
||||
double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
|
||||
double theta,u,mdu; //mdu: minus du, -du/dx=f
|
||||
int i1, i2, i3, n, type;
|
||||
double eangle, f1[3], f3[3];
|
||||
double delx1, dely1, delz1, delx2, dely2, delz2;
|
||||
double rsq1, rsq2, r1, r2, c, s, a, a11, a12, a22;
|
||||
double theta, u, mdu; //mdu: minus du, -du/dx=f
|
||||
|
||||
eangle = 0.0;
|
||||
ev_init(eflag,vflag);
|
||||
ev_init(eflag, vflag);
|
||||
|
||||
double **x = atom->x;
|
||||
double **f = atom->f;
|
||||
@ -97,7 +94,7 @@ void AngleTable::compute(int eflag, int vflag)
|
||||
dely1 = x[i1][1] - x[i2][1];
|
||||
delz1 = x[i1][2] - x[i2][2];
|
||||
|
||||
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1;
|
||||
rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
|
||||
r1 = sqrt(rsq1);
|
||||
|
||||
// 2nd bond
|
||||
@ -106,39 +103,39 @@ void AngleTable::compute(int eflag, int vflag)
|
||||
dely2 = x[i3][1] - x[i2][1];
|
||||
delz2 = x[i3][2] - x[i2][2];
|
||||
|
||||
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2;
|
||||
rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
|
||||
r2 = sqrt(rsq2);
|
||||
|
||||
// angle (cos and sin)
|
||||
|
||||
c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
s = sqrt(1.0 - c*c);
|
||||
s = sqrt(1.0 - c * c);
|
||||
if (s < SMALL) s = SMALL;
|
||||
s = 1.0/s;
|
||||
s = 1.0 / s;
|
||||
|
||||
// tabulated force & energy
|
||||
|
||||
theta = acos(c);
|
||||
uf_lookup(type,theta,u,mdu);
|
||||
uf_lookup(type, theta, u, mdu);
|
||||
|
||||
if (eflag) eangle = u;
|
||||
|
||||
a = mdu * s;
|
||||
a11 = a*c / rsq1;
|
||||
a12 = -a / (r1*r2);
|
||||
a22 = a*c / rsq2;
|
||||
a11 = a * c / rsq1;
|
||||
a12 = -a / (r1 * r2);
|
||||
a22 = a * c / rsq2;
|
||||
|
||||
f1[0] = a11*delx1 + a12*delx2;
|
||||
f1[1] = a11*dely1 + a12*dely2;
|
||||
f1[2] = a11*delz1 + a12*delz2;
|
||||
f3[0] = a22*delx2 + a12*delx1;
|
||||
f3[1] = a22*dely2 + a12*dely1;
|
||||
f3[2] = a22*delz2 + a12*delz1;
|
||||
f1[0] = a11 * delx1 + a12 * delx2;
|
||||
f1[1] = a11 * dely1 + a12 * dely2;
|
||||
f1[2] = a11 * delz1 + a12 * delz2;
|
||||
f3[0] = a22 * delx2 + a12 * delx1;
|
||||
f3[1] = a22 * dely2 + a12 * dely1;
|
||||
f3[2] = a22 * delz2 + a12 * delz1;
|
||||
|
||||
// apply force to each of 3 atoms
|
||||
|
||||
@ -160,8 +157,9 @@ void AngleTable::compute(int eflag, int vflag)
|
||||
f[i3][2] += f3[2];
|
||||
}
|
||||
|
||||
if (evflag) ev_tally(i1,i2,i3,nlocal,newton_bond,eangle,f1,f3,
|
||||
delx1,dely1,delz1,delx2,dely2,delz2);
|
||||
if (evflag)
|
||||
ev_tally(i1, i2, i3, nlocal, newton_bond, eangle, f1, f3, delx1, dely1, delz1, delx2, dely2,
|
||||
delz2);
|
||||
}
|
||||
}
|
||||
|
||||
@ -170,13 +168,13 @@ void AngleTable::compute(int eflag, int vflag)
|
||||
void AngleTable::allocate()
|
||||
{
|
||||
allocated = 1;
|
||||
int n = atom->nangletypes;
|
||||
const int np1 = atom->nangletypes + 1;
|
||||
|
||||
memory->create(theta0,n+1,"angle:theta0");
|
||||
memory->create(tabindex,n+1,"angle:tabindex");
|
||||
memory->create(theta0, np1, "angle:theta0");
|
||||
memory->create(tabindex, np1, "angle:tabindex");
|
||||
|
||||
memory->create(setflag,n+1,"angle:setflag");
|
||||
for (int i = 1; i <= n; i++) setflag[i] = 0;
|
||||
memory->create(setflag, np1, "angle:setflag");
|
||||
for (int i = 1; i < np1; i++) setflag[i] = 0;
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -185,14 +183,17 @@ void AngleTable::allocate()
|
||||
|
||||
void AngleTable::settings(int narg, char **arg)
|
||||
{
|
||||
if (narg != 2) error->all(FLERR,"Illegal angle_style command");
|
||||
if (narg != 2) error->all(FLERR, "Illegal angle_style command");
|
||||
|
||||
if (strcmp(arg[0],"linear") == 0) tabstyle = LINEAR;
|
||||
else if (strcmp(arg[0],"spline") == 0) tabstyle = SPLINE;
|
||||
else error->all(FLERR,"Unknown table style in angle style table");
|
||||
if (strcmp(arg[0], "linear") == 0)
|
||||
tabstyle = LINEAR;
|
||||
else if (strcmp(arg[0], "spline") == 0)
|
||||
tabstyle = SPLINE;
|
||||
else
|
||||
error->all(FLERR, "Unknown table style in angle style table");
|
||||
|
||||
tablength = utils::inumeric(FLERR,arg[1],false,lmp);
|
||||
if (tablength < 2) error->all(FLERR,"Illegal number of angle table entries");
|
||||
tablength = utils::inumeric(FLERR, arg[1], false, lmp);
|
||||
if (tablength < 2) error->all(FLERR, "Illegal number of angle table entries");
|
||||
|
||||
// delete old tables, since cannot just change settings
|
||||
|
||||
@ -200,8 +201,8 @@ void AngleTable::settings(int narg, char **arg)
|
||||
memory->sfree(tables);
|
||||
|
||||
if (allocated) {
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(tabindex);
|
||||
memory->destroy(setflag);
|
||||
memory->destroy(tabindex);
|
||||
}
|
||||
allocated = 0;
|
||||
|
||||
@ -215,36 +216,35 @@ void AngleTable::settings(int narg, char **arg)
|
||||
|
||||
void AngleTable::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 3) error->all(FLERR,"Illegal angle_coeff command");
|
||||
if (narg != 3) error->all(FLERR, "Illegal angle_coeff command");
|
||||
if (!allocated) allocate();
|
||||
|
||||
int ilo,ihi;
|
||||
utils::bounds(FLERR,arg[0],1,atom->nangletypes,ilo,ihi,error);
|
||||
int ilo, ihi;
|
||||
utils::bounds(FLERR, arg[0], 1, atom->nangletypes, ilo, ihi, error);
|
||||
|
||||
int me;
|
||||
MPI_Comm_rank(world,&me);
|
||||
tables = (Table *)
|
||||
memory->srealloc(tables,(ntables+1)*sizeof(Table),"angle:tables");
|
||||
MPI_Comm_rank(world, &me);
|
||||
tables = (Table *) memory->srealloc(tables, (ntables + 1) * sizeof(Table), "angle:tables");
|
||||
Table *tb = &tables[ntables];
|
||||
null_table(tb);
|
||||
if (me == 0) read_table(tb,arg[1],arg[2]);
|
||||
if (me == 0) read_table(tb, arg[1], arg[2]);
|
||||
bcast_table(tb);
|
||||
|
||||
// error check on table parameters
|
||||
|
||||
if (tb->ninput <= 1) error->one(FLERR,"Invalid angle table length");
|
||||
if (tb->ninput <= 1) error->one(FLERR, "Invalid angle table length");
|
||||
|
||||
double alo,ahi;
|
||||
double alo, ahi;
|
||||
alo = tb->afile[0];
|
||||
ahi = tb->afile[tb->ninput-1];
|
||||
if (fabs(alo-0.0) > TINY || fabs(ahi-180.0) > TINY)
|
||||
error->all(FLERR,"Angle table must range from 0 to 180 degrees");
|
||||
ahi = tb->afile[tb->ninput - 1];
|
||||
if (fabs(alo - 0.0) > TINY || fabs(ahi - 180.0) > TINY)
|
||||
error->all(FLERR, "Angle table must range from 0 to 180 degrees");
|
||||
|
||||
// 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] *= MY_PI / 180.0;
|
||||
tb->ffile[i] *= 180.0 / MY_PI;
|
||||
}
|
||||
|
||||
// spline read-in and compute a,e,f vectors within table
|
||||
@ -263,7 +263,7 @@ void AngleTable::coeff(int narg, char **arg)
|
||||
}
|
||||
ntables++;
|
||||
|
||||
if (count == 0) error->all(FLERR,"Illegal angle_coeff command");
|
||||
if (count == 0) error->all(FLERR, "Illegal angle_coeff command");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -301,8 +301,8 @@ void AngleTable::read_restart(FILE *fp)
|
||||
|
||||
void AngleTable::write_restart_settings(FILE *fp)
|
||||
{
|
||||
fwrite(&tabstyle,sizeof(int),1,fp);
|
||||
fwrite(&tablength,sizeof(int),1,fp);
|
||||
fwrite(&tabstyle, sizeof(int), 1, fp);
|
||||
fwrite(&tablength, sizeof(int), 1, fp);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -312,11 +312,11 @@ void AngleTable::write_restart_settings(FILE *fp)
|
||||
void AngleTable::read_restart_settings(FILE *fp)
|
||||
{
|
||||
if (comm->me == 0) {
|
||||
utils::sfread(FLERR,&tabstyle,sizeof(int),1,fp,nullptr,error);
|
||||
utils::sfread(FLERR,&tablength,sizeof(int),1,fp,nullptr,error);
|
||||
utils::sfread(FLERR, &tabstyle, sizeof(int), 1, fp, nullptr, error);
|
||||
utils::sfread(FLERR, &tablength, sizeof(int), 1, fp, nullptr, error);
|
||||
}
|
||||
MPI_Bcast(&tabstyle,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&tablength,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&tabstyle, 1, MPI_INT, 0, world);
|
||||
MPI_Bcast(&tablength, 1, MPI_INT, 0, world);
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -328,23 +328,23 @@ double AngleTable::single(int type, int i1, int i2, int i3)
|
||||
double delx1 = x[i1][0] - x[i2][0];
|
||||
double dely1 = x[i1][1] - x[i2][1];
|
||||
double delz1 = x[i1][2] - x[i2][2];
|
||||
domain->minimum_image(delx1,dely1,delz1);
|
||||
double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1);
|
||||
domain->minimum_image(delx1, dely1, delz1);
|
||||
double r1 = sqrt(delx1 * delx1 + dely1 * dely1 + delz1 * delz1);
|
||||
|
||||
double delx2 = x[i3][0] - x[i2][0];
|
||||
double dely2 = x[i3][1] - x[i2][1];
|
||||
double delz2 = x[i3][2] - x[i2][2];
|
||||
domain->minimum_image(delx2,dely2,delz2);
|
||||
double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2);
|
||||
domain->minimum_image(delx2, dely2, delz2);
|
||||
double r2 = sqrt(delx2 * delx2 + dely2 * dely2 + delz2 * delz2);
|
||||
|
||||
double c = delx1*delx2 + dely1*dely2 + delz1*delz2;
|
||||
c /= r1*r2;
|
||||
double c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
|
||||
c /= r1 * r2;
|
||||
if (c > 1.0) c = 1.0;
|
||||
if (c < -1.0) c = -1.0;
|
||||
|
||||
double theta = acos(c);
|
||||
double u=0.0;
|
||||
u_lookup(type,theta,u);
|
||||
double u = 0.0;
|
||||
u_lookup(type, theta, u);
|
||||
return u;
|
||||
}
|
||||
|
||||
@ -385,11 +385,9 @@ void AngleTable::read_table(Table *tb, char *file, char *keyword)
|
||||
{
|
||||
TableFileReader reader(lmp, file, "angle");
|
||||
|
||||
char * line = reader.find_section_start(keyword);
|
||||
char *line = reader.find_section_start(keyword);
|
||||
|
||||
if (!line) {
|
||||
error->one(FLERR,"Did not find keyword in table file");
|
||||
}
|
||||
if (!line) { error->one(FLERR, "Did not find keyword in table file"); }
|
||||
|
||||
// read args on 2nd line of section
|
||||
// allocate table arrays for file values
|
||||
@ -419,10 +417,9 @@ void AngleTable::read_table(Table *tb, char *file, char *keyword)
|
||||
|
||||
// warn if data was read incompletely, e.g. columns were missing
|
||||
|
||||
if (cerror) {
|
||||
std::string str = fmt::format("{} of {} lines in table were incomplete or could not be parsed completely", cerror, tb->ninput);
|
||||
error->warning(FLERR,str.c_str());
|
||||
}
|
||||
if (cerror)
|
||||
error->warning(FLERR, "{} of {} lines in table incomplete or could not be parsed", cerror,
|
||||
tb->ninput);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -432,22 +429,22 @@ void AngleTable::read_table(Table *tb, char *file, char *keyword)
|
||||
|
||||
void AngleTable::spline_table(Table *tb)
|
||||
{
|
||||
memory->create(tb->e2file,tb->ninput,"angle:e2file");
|
||||
memory->create(tb->f2file,tb->ninput,"angle:f2file");
|
||||
memory->create(tb->e2file, tb->ninput, "angle:e2file");
|
||||
memory->create(tb->f2file, tb->ninput, "angle:f2file");
|
||||
|
||||
double ep0 = - tb->ffile[0];
|
||||
double epn = - tb->ffile[tb->ninput-1];
|
||||
spline(tb->afile,tb->efile,tb->ninput,ep0,epn,tb->e2file);
|
||||
double ep0 = -tb->ffile[0];
|
||||
double epn = -tb->ffile[tb->ninput - 1];
|
||||
spline(tb->afile, tb->efile, tb->ninput, ep0, epn, tb->e2file);
|
||||
|
||||
if (tb->fpflag == 0) {
|
||||
tb->fplo = (tb->ffile[1] - tb->ffile[0]) / (tb->afile[1] - tb->afile[0]);
|
||||
tb->fphi = (tb->ffile[tb->ninput-1] - tb->ffile[tb->ninput-2]) /
|
||||
(tb->afile[tb->ninput-1] - tb->afile[tb->ninput-2]);
|
||||
tb->fphi = (tb->ffile[tb->ninput - 1] - tb->ffile[tb->ninput - 2]) /
|
||||
(tb->afile[tb->ninput - 1] - tb->afile[tb->ninput - 2]);
|
||||
}
|
||||
|
||||
double fp0 = tb->fplo;
|
||||
double fpn = tb->fphi;
|
||||
spline(tb->afile,tb->ffile,tb->ninput,fp0,fpn,tb->f2file);
|
||||
spline(tb->afile, tb->ffile, tb->ninput, fp0, fpn, tb->f2file);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -458,44 +455,44 @@ void AngleTable::compute_table(Table *tb)
|
||||
{
|
||||
// delta = table spacing in angle for N-1 bins
|
||||
|
||||
int tlm1 = tablength-1;
|
||||
int tlm1 = tablength - 1;
|
||||
tb->delta = MY_PI / tlm1;
|
||||
tb->invdelta = 1.0/tb->delta;
|
||||
tb->deltasq6 = tb->delta*tb->delta / 6.0;
|
||||
tb->invdelta = 1.0 / tb->delta;
|
||||
tb->deltasq6 = tb->delta * tb->delta / 6.0;
|
||||
|
||||
// N-1 evenly spaced bins in angle from 0 to PI
|
||||
// ang,e,f = value at lower edge of bin
|
||||
// de,df values = delta values of e,f
|
||||
// ang,e,f are N in length so de,df arrays can compute difference
|
||||
|
||||
memory->create(tb->ang,tablength,"angle:ang");
|
||||
memory->create(tb->e,tablength,"angle:e");
|
||||
memory->create(tb->de,tablength,"angle:de");
|
||||
memory->create(tb->f,tablength,"angle:f");
|
||||
memory->create(tb->df,tablength,"angle:df");
|
||||
memory->create(tb->e2,tablength,"angle:e2");
|
||||
memory->create(tb->f2,tablength,"angle:f2");
|
||||
memory->create(tb->ang, tablength, "angle:ang");
|
||||
memory->create(tb->e, tablength, "angle:e");
|
||||
memory->create(tb->de, tablength, "angle:de");
|
||||
memory->create(tb->f, tablength, "angle:f");
|
||||
memory->create(tb->df, tablength, "angle:df");
|
||||
memory->create(tb->e2, tablength, "angle:e2");
|
||||
memory->create(tb->f2, tablength, "angle:f2");
|
||||
|
||||
double a;
|
||||
for (int i = 0; i < tablength; i++) {
|
||||
a = i*tb->delta;
|
||||
a = i * tb->delta;
|
||||
tb->ang[i] = a;
|
||||
tb->e[i] = splint(tb->afile,tb->efile,tb->e2file,tb->ninput,a);
|
||||
tb->f[i] = splint(tb->afile,tb->ffile,tb->f2file,tb->ninput,a);
|
||||
tb->e[i] = splint(tb->afile, tb->efile, tb->e2file, tb->ninput, a);
|
||||
tb->f[i] = splint(tb->afile, tb->ffile, tb->f2file, tb->ninput, a);
|
||||
}
|
||||
|
||||
for (int i = 0; i < tlm1; i++) {
|
||||
tb->de[i] = tb->e[i+1] - tb->e[i];
|
||||
tb->df[i] = tb->f[i+1] - tb->f[i];
|
||||
tb->de[i] = tb->e[i + 1] - tb->e[i];
|
||||
tb->df[i] = tb->f[i + 1] - tb->f[i];
|
||||
}
|
||||
// get final elements from linear extrapolation
|
||||
tb->de[tlm1] = 2.0*tb->de[tlm1-1] - tb->de[tlm1-2];
|
||||
tb->df[tlm1] = 2.0*tb->df[tlm1-1] - tb->df[tlm1-2];
|
||||
tb->de[tlm1] = 2.0 * tb->de[tlm1 - 1] - tb->de[tlm1 - 2];
|
||||
tb->df[tlm1] = 2.0 * tb->df[tlm1 - 1] - tb->df[tlm1 - 2];
|
||||
|
||||
double ep0 = - tb->f[0];
|
||||
double epn = - tb->f[tlm1];
|
||||
spline(tb->ang,tb->e,tablength,ep0,epn,tb->e2);
|
||||
spline(tb->ang,tb->f,tablength,tb->fplo,tb->fphi,tb->f2);
|
||||
double ep0 = -tb->f[0];
|
||||
double epn = -tb->f[tlm1];
|
||||
spline(tb->ang, tb->e, tablength, ep0, epn, tb->e2);
|
||||
spline(tb->ang, tb->f, tablength, tb->fplo, tb->fphi, tb->f2);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -522,19 +519,19 @@ 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 *= (180.0 / MY_PI) * (180.0 / MY_PI);
|
||||
tb->fphi *= (180.0 / MY_PI) * (180.0 / MY_PI);
|
||||
} else if (word == "EQ") {
|
||||
tb->theta0 = values.next_double()/180.0*MY_PI;
|
||||
tb->theta0 = values.next_double() / 180.0 * MY_PI;
|
||||
} else {
|
||||
error->one(FLERR,"Invalid keyword in angle table parameters");
|
||||
error->one(FLERR, "Invalid keyword in angle table parameters");
|
||||
}
|
||||
}
|
||||
} catch(TokenizerException &e) {
|
||||
} catch (TokenizerException &e) {
|
||||
error->one(FLERR, e.what());
|
||||
}
|
||||
|
||||
if (tb->ninput == 0) error->one(FLERR,"Angle table parameters did not set N");
|
||||
if (tb->ninput == 0) error->one(FLERR, "Angle table parameters did not set N");
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
@ -545,80 +542,84 @@ void AngleTable::param_extract(Table *tb, char *line)
|
||||
|
||||
void AngleTable::bcast_table(Table *tb)
|
||||
{
|
||||
MPI_Bcast(&tb->ninput,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&tb->ninput, 1, MPI_INT, 0, world);
|
||||
|
||||
int me;
|
||||
MPI_Comm_rank(world,&me);
|
||||
MPI_Comm_rank(world, &me);
|
||||
if (me > 0) {
|
||||
memory->create(tb->afile,tb->ninput,"angle:afile");
|
||||
memory->create(tb->efile,tb->ninput,"angle:efile");
|
||||
memory->create(tb->ffile,tb->ninput,"angle:ffile");
|
||||
memory->create(tb->afile, tb->ninput, "angle:afile");
|
||||
memory->create(tb->efile, tb->ninput, "angle:efile");
|
||||
memory->create(tb->ffile, tb->ninput, "angle:ffile");
|
||||
}
|
||||
|
||||
MPI_Bcast(tb->afile,tb->ninput,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(tb->efile,tb->ninput,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(tb->ffile,tb->ninput,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(tb->afile, tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->efile, tb->ninput, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(tb->ffile, tb->ninput, MPI_DOUBLE, 0, world);
|
||||
|
||||
MPI_Bcast(&tb->fpflag,1,MPI_INT,0,world);
|
||||
MPI_Bcast(&tb->fpflag, 1, MPI_INT, 0, world);
|
||||
if (tb->fpflag) {
|
||||
MPI_Bcast(&tb->fplo,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&tb->fphi,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&tb->fplo, 1, MPI_DOUBLE, 0, world);
|
||||
MPI_Bcast(&tb->fphi, 1, MPI_DOUBLE, 0, world);
|
||||
}
|
||||
MPI_Bcast(&tb->theta0,1,MPI_DOUBLE,0,world);
|
||||
MPI_Bcast(&tb->theta0, 1, MPI_DOUBLE, 0, world);
|
||||
}
|
||||
|
||||
/* ----------------------------------------------------------------------
|
||||
spline and splint routines modified from Numerical Recipes
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
void AngleTable::spline(double *x, double *y, int n,
|
||||
double yp1, double ypn, double *y2)
|
||||
void AngleTable::spline(double *x, double *y, int n, double yp1, double ypn, double *y2)
|
||||
{
|
||||
int i,k;
|
||||
double p,qn,sig,un;
|
||||
int i, k;
|
||||
double p, qn, sig, un;
|
||||
double *u = new double[n];
|
||||
|
||||
if (yp1 > 0.99e300) y2[0] = u[0] = 0.0;
|
||||
if (yp1 > 0.99e300)
|
||||
y2[0] = u[0] = 0.0;
|
||||
else {
|
||||
y2[0] = -0.5;
|
||||
u[0] = (3.0/(x[1]-x[0])) * ((y[1]-y[0]) / (x[1]-x[0]) - yp1);
|
||||
u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - yp1);
|
||||
}
|
||||
for (i = 1; i < n-1; i++) {
|
||||
sig = (x[i]-x[i-1]) / (x[i+1]-x[i-1]);
|
||||
p = sig*y2[i-1] + 2.0;
|
||||
y2[i] = (sig-1.0) / p;
|
||||
u[i] = (y[i+1]-y[i]) / (x[i+1]-x[i]) - (y[i]-y[i-1]) / (x[i]-x[i-1]);
|
||||
u[i] = (6.0*u[i] / (x[i+1]-x[i-1]) - sig*u[i-1]) / p;
|
||||
for (i = 1; i < n - 1; i++) {
|
||||
sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
|
||||
p = sig * y2[i - 1] + 2.0;
|
||||
y2[i] = (sig - 1.0) / p;
|
||||
u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
|
||||
u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
|
||||
}
|
||||
if (ypn > 0.99e300) qn = un = 0.0;
|
||||
if (ypn > 0.99e300)
|
||||
qn = un = 0.0;
|
||||
else {
|
||||
qn = 0.5;
|
||||
un = (3.0/(x[n-1]-x[n-2])) * (ypn - (y[n-1]-y[n-2]) / (x[n-1]-x[n-2]));
|
||||
un = (3.0 / (x[n - 1] - x[n - 2])) * (ypn - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
|
||||
}
|
||||
y2[n-1] = (un-qn*u[n-2]) / (qn*y2[n-2] + 1.0);
|
||||
for (k = n-2; k >= 0; k--) y2[k] = y2[k]*y2[k+1] + u[k];
|
||||
y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
|
||||
for (k = n - 2; k >= 0; k--) y2[k] = y2[k] * y2[k + 1] + u[k];
|
||||
|
||||
delete [] u;
|
||||
delete[] u;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double AngleTable::splint(double *xa, double *ya, double *y2a, int n, double x)
|
||||
{
|
||||
int klo,khi,k;
|
||||
double h,b,a,y;
|
||||
int klo, khi, k;
|
||||
double h, b, a, y;
|
||||
|
||||
klo = 0;
|
||||
khi = n-1;
|
||||
while (khi-klo > 1) {
|
||||
k = (khi+klo) >> 1;
|
||||
if (xa[k] > x) khi = k;
|
||||
else klo = k;
|
||||
khi = n - 1;
|
||||
while (khi - klo > 1) {
|
||||
k = (khi + klo) >> 1;
|
||||
if (xa[k] > x)
|
||||
khi = k;
|
||||
else
|
||||
klo = k;
|
||||
}
|
||||
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;
|
||||
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;
|
||||
return y;
|
||||
}
|
||||
|
||||
@ -628,31 +629,29 @@ double AngleTable::splint(double *xa, double *ya, double *y2a, int n, double x)
|
||||
|
||||
void AngleTable::uf_lookup(int type, double x, double &u, double &f)
|
||||
{
|
||||
if (!std::isfinite(x)) {
|
||||
error->one(FLERR,"Illegal angle in angle style table");
|
||||
}
|
||||
if (!std::isfinite(x)) { error->one(FLERR, "Illegal angle in angle style table"); }
|
||||
|
||||
double fraction,a,b;
|
||||
double fraction, a, b;
|
||||
const Table *tb = &tables[tabindex[type]];
|
||||
|
||||
// invdelta is based on tablength-1
|
||||
int itable = static_cast<int> (x * tb->invdelta);
|
||||
int itable = static_cast<int>(x * tb->invdelta);
|
||||
if (itable < 0) itable = 0;
|
||||
if (itable >= tablength) itable = tablength-1;
|
||||
if (itable >= tablength) itable = tablength - 1;
|
||||
|
||||
if (tabstyle == LINEAR) {
|
||||
fraction = (x - tb->ang[itable]) * tb->invdelta;
|
||||
u = tb->e[itable] + fraction*tb->de[itable];
|
||||
f = tb->f[itable] + fraction*tb->df[itable];
|
||||
u = tb->e[itable] + fraction * tb->de[itable];
|
||||
f = tb->f[itable] + fraction * tb->df[itable];
|
||||
} else if (tabstyle == SPLINE) {
|
||||
fraction = (x - tb->ang[itable]) * tb->invdelta;
|
||||
|
||||
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;
|
||||
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;
|
||||
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;
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
||||
@ -662,26 +661,24 @@ void AngleTable::uf_lookup(int type, double x, double &u, double &f)
|
||||
|
||||
void AngleTable::u_lookup(int type, double x, double &u)
|
||||
{
|
||||
if (!std::isfinite(x)) {
|
||||
error->one(FLERR,"Illegal angle in angle style table");
|
||||
}
|
||||
if (!std::isfinite(x)) { error->one(FLERR, "Illegal angle in angle style table"); }
|
||||
|
||||
double fraction,a,b;
|
||||
double fraction, a, b;
|
||||
const Table *tb = &tables[tabindex[type]];
|
||||
int itable = static_cast<int> ( x * tb->invdelta);
|
||||
int itable = static_cast<int>(x * tb->invdelta);
|
||||
|
||||
if (itable < 0) itable = 0;
|
||||
if (itable >= tablength) itable = tablength-1;
|
||||
if (itable >= tablength) itable = tablength - 1;
|
||||
|
||||
if (tabstyle == LINEAR) {
|
||||
fraction = (x - tb->ang[itable]) * tb->invdelta;
|
||||
u = tb->e[itable] + fraction*tb->de[itable];
|
||||
u = tb->e[itable] + fraction * tb->de[itable];
|
||||
} else if (tabstyle == SPLINE) {
|
||||
fraction = (x - tb->ang[itable]) * tb->invdelta;
|
||||
|
||||
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;
|
||||
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;
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user