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

apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-10-01 00:58:38 -04:00
parent 01fb33cb5d
commit 434c170097
2 changed files with 180 additions and 181 deletions
Download Patch File Download Diff File Expand all files Collapse all files

189
src/EXTRA-MOLECULE/angle_gaussian.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
@ -34,9 +33,9 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
AngleGaussian::AngleGaussian(LAMMPS *lmp)
: Angle(lmp), nterms(nullptr), angle_temperature(nullptr),
alpha(nullptr), width(nullptr), theta0(nullptr)
AngleGaussian::AngleGaussian(LAMMPS *lmp) :
Angle(lmp), nterms(nullptr), angle_temperature(nullptr), alpha(nullptr), width(nullptr),
theta0(nullptr)
{
}
@ -49,13 +48,13 @@ AngleGaussian::~AngleGaussian()
memory->destroy(nterms);
memory->destroy(angle_temperature);
for (int i = 1; i <= atom->nangletypes; i++) {
delete [] alpha[i];
delete [] width[i];
delete [] theta0[i];
delete[] alpha[i];
delete[] width[i];
delete[] theta0[i];
}
delete [] alpha;
delete [] width;
delete [] theta0;
delete[] alpha;
delete[] width;
delete[] theta0;
}
}
@ -63,15 +62,15 @@ AngleGaussian::~AngleGaussian()
void AngleGaussian::compute(int eflag, int vflag)
{
int i1,i2,i3,n,type;
double delx1,dely1,delz1,delx2,dely2,delz2;
double eangle,f1[3],f3[3];
int i1, i2, i3, n, type;
double delx1, dely1, delz1, delx2, dely2, delz2;
double eangle, f1[3], f3[3];
double dtheta;
double rsq1,rsq2,r1,r2,c,s,a,a11,a12,a22;
double rsq1, rsq2, r1, r2, c, s, a, a11, a12, a22;
double prefactor, exponent, g_i, sum_g_i, sum_numerator;
eangle = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
@ -92,7 +91,7 @@ void AngleGaussian::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
@ -101,20 +100,20 @@ void AngleGaussian::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 < SMAL) s = SMAL;
s = 1.0/s;
s = 1.0 / s;
// force & energy
double theta = acos(c);
@ -123,28 +122,28 @@ void AngleGaussian::compute(int eflag, int vflag)
sum_numerator = 0.0;
for (int i = 0; i < nterms[type]; i++) {
dtheta = theta - theta0[type][i];
prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2)));
exponent = -2*dtheta*dtheta/(width[type][i]*width[type][i]);
g_i = prefactor*exp(exponent);
prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2)));
exponent = -2 * dtheta * dtheta / (width[type][i] * width[type][i]);
g_i = prefactor * exp(exponent);
sum_g_i += g_i;
sum_numerator += g_i*dtheta/(width[type][i]*width[type][i]);
sum_numerator += g_i * dtheta / (width[type][i] * width[type][i]);
}
if (sum_g_i < SMALL) sum_g_i = SMALL;
if (eflag) eangle = -(force->boltz*angle_temperature[type])*log(sum_g_i);
if (eflag) eangle = -(force->boltz * angle_temperature[type]) * log(sum_g_i);
// I should check about the sign of this expression
a = -4.0*(force->boltz*angle_temperature[type])*(sum_numerator/sum_g_i)*s;
a11 = a*c / rsq1;
a12 = -a / (r1*r2);
a22 = a*c / rsq2;
a = -4.0 * (force->boltz * angle_temperature[type]) * (sum_numerator / sum_g_i) * s;
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
@ -166,8 +165,9 @@ void AngleGaussian::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);
}
}
@ -176,20 +176,20 @@ void AngleGaussian::compute(int eflag, int vflag)
void AngleGaussian::allocate()
{
allocated = 1;
int n = atom->nangletypes+1;
int n = atom->nangletypes + 1;
memory->create(nterms,n,"angle:nterms");
memory->create(angle_temperature,n,"angle:angle_temperature");
memory->create(nterms, n, "angle:nterms");
memory->create(angle_temperature, n, "angle:angle_temperature");
alpha = new double*[n];
width = new double*[n];
theta0 = new double*[n];
memset(alpha,0,sizeof(double *)*n);
memset(width,0,sizeof(double *)*n);
memset(theta0,0,sizeof(double *)*n);
alpha = new double *[n];
width = new double *[n];
theta0 = new double *[n];
memset(alpha, 0, sizeof(double *) * n);
memset(width, 0, sizeof(double *) * n);
memset(theta0, 0, sizeof(double *) * n);
memory->create(setflag,n,"angle:setflag");
memset(setflag,0,sizeof(int)*n);
memory->create(setflag, n, "angle:setflag");
memset(setflag, 0, sizeof(int) * n);
}
/* ----------------------------------------------------------------------
@ -198,15 +198,14 @@ void AngleGaussian::allocate()
void AngleGaussian::coeff(int narg, char **arg)
{
if (narg < 6) error->all(FLERR,"Incorrect args for angle coefficients");
if (narg < 6) error->all(FLERR, "Incorrect args for angle coefficients");
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 angle_temperature_one = utils::numeric(FLERR,arg[1],false,lmp);
int n = utils::inumeric(FLERR,arg[2],false,lmp);
if (narg != 3*n + 3)
error->all(FLERR,"Incorrect args for angle coefficients");
double angle_temperature_one = utils::numeric(FLERR, arg[1], false, lmp);
int n = utils::inumeric(FLERR, arg[2], false, lmp);
if (narg != 3 * n + 3) error->all(FLERR, "Incorrect args for angle coefficients");
if (!allocated) allocate();
@ -217,21 +216,21 @@ void AngleGaussian::coeff(int narg, char **arg)
angle_temperature[i] = angle_temperature_one;
nterms[i] = n;
delete[] alpha[i];
alpha[i] = new double [n];
alpha[i] = new double[n];
delete[] width[i];
width[i] = new double [n];
width[i] = new double[n];
delete[] theta0[i];
theta0[i] = new double [n];
theta0[i] = new double[n];
for (int j = 0; j < n; j++) {
alpha[i][j] = utils::numeric(FLERR,arg[3+3*j],false,lmp);
width[i][j] = utils::numeric(FLERR,arg[4+3*j],false,lmp);
theta0[i][j] = utils::numeric(FLERR,arg[5+3*j],false,lmp)* MY_PI / 180.0;
alpha[i][j] = utils::numeric(FLERR, arg[3 + 3 * j], false, lmp);
width[i][j] = utils::numeric(FLERR, arg[4 + 3 * j], false, lmp);
theta0[i][j] = utils::numeric(FLERR, arg[5 + 3 * j], false, lmp) * MY_PI / 180.0;
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");
}
/* ---------------------------------------------------------------------- */
@ -247,12 +246,12 @@ double AngleGaussian::equilibrium_angle(int i)
void AngleGaussian::write_restart(FILE *fp)
{
fwrite(&angle_temperature[1],sizeof(double),atom->nangletypes,fp);
fwrite(&nterms[1],sizeof(int),atom->nangletypes,fp);
fwrite(&angle_temperature[1], sizeof(double), atom->nangletypes, fp);
fwrite(&nterms[1], sizeof(int), atom->nangletypes, fp);
for (int i = 1; i <= atom->nangletypes; i++) {
fwrite(alpha[i],sizeof(double),nterms[i],fp);
fwrite(width[i],sizeof(double),nterms[i],fp);
fwrite(theta0[i],sizeof(double),nterms[i],fp);
fwrite(alpha[i], sizeof(double), nterms[i], fp);
fwrite(width[i], sizeof(double), nterms[i], fp);
fwrite(theta0[i], sizeof(double), nterms[i], fp);
}
}
@ -265,31 +264,32 @@ void AngleGaussian::read_restart(FILE *fp)
allocate();
if (comm->me == 0) {
utils::sfread(FLERR,&angle_temperature[1],sizeof(double),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR,&nterms[1],sizeof(int),atom->nangletypes,fp,nullptr,error);
utils::sfread(FLERR, &angle_temperature[1], sizeof(double), atom->nangletypes, fp, nullptr,
error);
utils::sfread(FLERR, &nterms[1], sizeof(int), atom->nangletypes, fp, nullptr, error);
}
MPI_Bcast(&angle_temperature[1],atom->nangletypes,MPI_DOUBLE,0,world);
MPI_Bcast(&nterms[1],atom->nangletypes,MPI_INT,0,world);
MPI_Bcast(&angle_temperature[1], atom->nangletypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&nterms[1], atom->nangletypes, MPI_INT, 0, world);
// allocate
for (int i = 1; i <= atom->nangletypes; i++) {
alpha[i] = new double [nterms[i]];
width[i] = new double [nterms[i]];
theta0[i] = new double [nterms[i]];
alpha[i] = new double[nterms[i]];
width[i] = new double[nterms[i]];
theta0[i] = new double[nterms[i]];
}
if (comm->me == 0) {
for (int i = 1; i <= atom->nangletypes; i++) {
utils::sfread(FLERR,alpha[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR,width[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR,theta0[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR, alpha[i], sizeof(double), nterms[i], fp, nullptr, error);
utils::sfread(FLERR, width[i], sizeof(double), nterms[i], fp, nullptr, error);
utils::sfread(FLERR, theta0[i], sizeof(double), nterms[i], fp, nullptr, error);
}
}
for (int i = 1; i <= atom->nangletypes; i++) {
MPI_Bcast(alpha[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(width[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(theta0[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(alpha[i], nterms[i], MPI_DOUBLE, 0, world);
MPI_Bcast(width[i], nterms[i], MPI_DOUBLE, 0, world);
MPI_Bcast(theta0[i], nterms[i], MPI_DOUBLE, 0, world);
}
for (int i = 1; i <= atom->nangletypes; i++) setflag[i] = 1;
@ -302,13 +302,12 @@ void AngleGaussian::read_restart(FILE *fp)
void AngleGaussian::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nangletypes; i++) {
fprintf(fp,"%d %g %d",i,angle_temperature[i],nterms[i]);
fprintf(fp, "%d %g %d", i, angle_temperature[i], nterms[i]);
for (int j = 0; j < nterms[i]; j++) {
fprintf(fp," %g %g %g",alpha[i][j],width[i][j],(theta0[i][j]/MY_PI)*180.0);
fprintf(fp, " %g %g %g", alpha[i][j], width[i][j], (theta0[i][j] / MY_PI) * 180.0);
}
fprintf(fp, "\n");
}
}
/* ---------------------------------------------------------------------- */
@ -320,30 +319,30 @@ double AngleGaussian::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 theta = acos(c);
double sum_g_i = 0.0;
for (int i = 0; i < nterms[type]; i++) {
double dtheta = theta - theta0[type][i];
double prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2)));
double exponent = -2*dtheta*dtheta/(width[type][i]*width[type][i]);
double g_i = prefactor*exp(exponent);
double prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2)));
double exponent = -2 * dtheta * dtheta / (width[type][i] * width[type][i]);
double g_i = prefactor * exp(exponent);
sum_g_i += g_i;
}
if (sum_g_i < SMALL) sum_g_i = SMALL;
return -(force->boltz*angle_temperature[type])*log(sum_g_i);
return -(force->boltz * angle_temperature[type]) * log(sum_g_i);
}

172
src/EXTRA-MOLECULE/bond_gaussian.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
@ -15,12 +14,12 @@
#include "bond_gaussian.h"
#include "atom.h"
#include "neighbor.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
@ -32,9 +31,9 @@ using namespace MathConst;
/* ---------------------------------------------------------------------- */
BondGaussian::BondGaussian(LAMMPS *lmp)
: Bond(lmp), nterms(nullptr), bond_temperature(nullptr),
alpha(nullptr), width(nullptr), r0(nullptr)
BondGaussian::BondGaussian(LAMMPS *lmp) :
Bond(lmp), nterms(nullptr), bond_temperature(nullptr), alpha(nullptr), width(nullptr),
r0(nullptr)
{
reinitflag = 1;
}
@ -48,13 +47,13 @@ BondGaussian::~BondGaussian()
memory->destroy(nterms);
memory->destroy(bond_temperature);
for (int i = 1; i <= atom->nbondtypes; i++) {
delete [] alpha[i];
delete [] width[i];
delete [] r0[i];
delete[] alpha[i];
delete[] width[i];
delete[] r0[i];
}
delete [] alpha;
delete [] width;
delete [] r0;
delete[] alpha;
delete[] width;
delete[] r0;
}
}
@ -62,13 +61,13 @@ BondGaussian::~BondGaussian()
void BondGaussian::compute(int eflag, int vflag)
{
int i1,i2,n,type;
double delx,dely,delz,ebond,fbond;
double rsq,r,dr;
int i1, i2, n, type;
double delx, dely, delz, ebond, fbond;
double rsq, r, dr;
double prefactor, exponent, g_i, sum_g_i, sum_numerator;
ebond = 0.0;
ev_init(eflag,vflag);
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
@ -86,43 +85,45 @@ void BondGaussian::compute(int eflag, int vflag)
dely = x[i1][1] - x[i2][1];
delz = x[i1][2] - x[i2][2];
rsq = delx*delx + dely*dely + delz*delz;
rsq = delx * delx + dely * dely + delz * delz;
r = sqrt(rsq);
sum_g_i = 0.0;
sum_numerator = 0.0;
for (int i = 0; i < nterms[type]; i++) {
dr = r - r0[type][i];
prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2)));
exponent = -2*dr*dr/(width[type][i]*width[type][i]);
g_i = prefactor*exp(exponent);
prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2)));
exponent = -2 * dr * dr / (width[type][i] * width[type][i]);
g_i = prefactor * exp(exponent);
sum_g_i += g_i;
sum_numerator += g_i*dr/(width[type][i]*width[type][i]);
sum_numerator += g_i * dr / (width[type][i] * width[type][i]);
}
// force & energy
if (sum_g_i < SMALL) sum_g_i = SMALL;
if (r > 0.0) fbond = -4.0*(force->boltz*bond_temperature[type])*(sum_numerator/sum_g_i)/r;
else fbond = 0.0;
if (r > 0.0)
fbond = -4.0 * (force->boltz * bond_temperature[type]) * (sum_numerator / sum_g_i) / r;
else
fbond = 0.0;
if (eflag) ebond = -(force->boltz*bond_temperature[type])*log(sum_g_i);
if (eflag) ebond = -(force->boltz * bond_temperature[type]) * log(sum_g_i);
// apply force to each of 2 atoms
if (newton_bond || i1 < nlocal) {
f[i1][0] += delx*fbond;
f[i1][1] += dely*fbond;
f[i1][2] += delz*fbond;
f[i1][0] += delx * fbond;
f[i1][1] += dely * fbond;
f[i1][2] += delz * fbond;
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= delx*fbond;
f[i2][1] -= dely*fbond;
f[i2][2] -= delz*fbond;
f[i2][0] -= delx * fbond;
f[i2][1] -= dely * fbond;
f[i2][2] -= delz * fbond;
}
if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz);
if (evflag) ev_tally(i1, i2, nlocal, newton_bond, ebond, fbond, delx, dely, delz);
}
}
@ -131,20 +132,20 @@ void BondGaussian::compute(int eflag, int vflag)
void BondGaussian::allocate()
{
allocated = 1;
int n = atom->nbondtypes+1;
int n = atom->nbondtypes + 1;
memory->create(nterms,n,"bond:nterms");
memory->create(bond_temperature,n,"bond:bond_temperature");
memory->create(nterms, n, "bond:nterms");
memory->create(bond_temperature, n, "bond:bond_temperature");
alpha = new double *[n];
width = new double *[n];
r0 = new double *[n];
memset(alpha,0,sizeof(double *)*n);
memset(width,0,sizeof(double *)*n);
memset(r0,0,sizeof(double *)*n);
memset(alpha, 0, sizeof(double *) * n);
memset(width, 0, sizeof(double *) * n);
memset(r0, 0, sizeof(double *) * n);
memory->create(setflag,n,"bond:setflag");
memset(setflag,0,sizeof(int)*n);
memory->create(setflag, n, "bond:setflag");
memset(setflag, 0, sizeof(int) * n);
}
/* ----------------------------------------------------------------------
@ -153,15 +154,14 @@ void BondGaussian::allocate()
void BondGaussian::coeff(int narg, char **arg)
{
if (narg < 6) error->all(FLERR,"Incorrect args for bond coefficients");
if (narg < 6) error->all(FLERR, "Incorrect args for bond coefficients");
int ilo,ihi;
utils::bounds(FLERR,arg[0],1,atom->nbondtypes,ilo,ihi,error);
int ilo, ihi;
utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error);
double bond_temp_one = utils::numeric(FLERR,arg[1],false,lmp);
int n = utils::inumeric(FLERR,arg[2],false,lmp);
if (narg != 3*n + 3)
error->all(FLERR,"Incorrect args for bond coefficients");
double bond_temp_one = utils::numeric(FLERR, arg[1], false, lmp);
int n = utils::inumeric(FLERR, arg[2], false, lmp);
if (narg != 3 * n + 3) error->all(FLERR, "Incorrect args for bond coefficients");
if (!allocated) allocate();
@ -170,21 +170,21 @@ void BondGaussian::coeff(int narg, char **arg)
bond_temperature[i] = bond_temp_one;
nterms[i] = n;
delete[] alpha[i];
alpha[i] = new double [n];
alpha[i] = new double[n];
delete[] width[i];
width[i] = new double [n];
width[i] = new double[n];
delete[] r0[i];
r0[i] = new double [n];
r0[i] = new double[n];
for (int j = 0; j < n; j++) {
alpha[i][j] = utils::numeric(FLERR,arg[3+3*j],false,lmp);
width[i][j] = utils::numeric(FLERR,arg[4+3*j],false,lmp);
r0[i][j] = utils::numeric(FLERR,arg[5+3*j],false,lmp);
alpha[i][j] = utils::numeric(FLERR, arg[3 + 3 * j], false, lmp);
width[i][j] = utils::numeric(FLERR, arg[4 + 3 * j], false, lmp);
r0[i][j] = utils::numeric(FLERR, arg[5 + 3 * j], false, lmp);
setflag[i] = 1;
}
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for bond coefficients");
if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients");
}
/* ----------------------------------------------------------------------
@ -202,12 +202,12 @@ double BondGaussian::equilibrium_distance(int i)
void BondGaussian::write_restart(FILE *fp)
{
fwrite(&bond_temperature[1],sizeof(double),atom->nbondtypes,fp);
fwrite(&nterms[1],sizeof(int),atom->nbondtypes,fp);
fwrite(&bond_temperature[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&nterms[1], sizeof(int), atom->nbondtypes, fp);
for (int i = 1; i <= atom->nbondtypes; i++) {
fwrite(alpha[i],sizeof(double),nterms[i],fp);
fwrite(width[i],sizeof(double),nterms[i],fp);
fwrite(r0[i],sizeof(double),nterms[i],fp);
fwrite(alpha[i], sizeof(double), nterms[i], fp);
fwrite(width[i], sizeof(double), nterms[i], fp);
fwrite(r0[i], sizeof(double), nterms[i], fp);
}
}
@ -220,31 +220,32 @@ void BondGaussian::read_restart(FILE *fp)
allocate();
if (comm->me == 0) {
utils::sfread(FLERR,&bond_temperature[1],sizeof(double),atom->nbondtypes,fp,nullptr,error);
utils::sfread(FLERR,&nterms[1],sizeof(int),atom->nbondtypes,fp,nullptr,error);
utils::sfread(FLERR, &bond_temperature[1], sizeof(double), atom->nbondtypes, fp, nullptr,
error);
utils::sfread(FLERR, &nterms[1], sizeof(int), atom->nbondtypes, fp, nullptr, error);
}
MPI_Bcast(&bond_temperature[1],atom->nbondtypes,MPI_DOUBLE,0,world);
MPI_Bcast(&nterms[1],atom->nbondtypes,MPI_INT,0,world);
MPI_Bcast(&bond_temperature[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&nterms[1], atom->nbondtypes, MPI_INT, 0, world);
// allocate
for (int i = 1; i <= atom->nbondtypes; i++) {
alpha[i] = new double [nterms[i]];
width[i] = new double [nterms[i]];
r0[i] = new double [nterms[i]];
alpha[i] = new double[nterms[i]];
width[i] = new double[nterms[i]];
r0[i] = new double[nterms[i]];
}
if (comm->me == 0) {
for (int i = 1; i <= atom->nbondtypes; i++) {
utils::sfread(FLERR,alpha[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR,width[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR,r0[i],sizeof(double),nterms[i],fp,nullptr,error);
utils::sfread(FLERR, alpha[i], sizeof(double), nterms[i], fp, nullptr, error);
utils::sfread(FLERR, width[i], sizeof(double), nterms[i], fp, nullptr, error);
utils::sfread(FLERR, r0[i], sizeof(double), nterms[i], fp, nullptr, error);
}
}
for (int i = 1; i <= atom->nbondtypes; i++) {
MPI_Bcast(alpha[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(width[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(r0[i],nterms[i],MPI_DOUBLE,0,world);
MPI_Bcast(alpha[i], nterms[i], MPI_DOUBLE, 0, world);
MPI_Bcast(width[i], nterms[i], MPI_DOUBLE, 0, world);
MPI_Bcast(r0[i], nterms[i], MPI_DOUBLE, 0, world);
}
for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
@ -257,9 +258,9 @@ void BondGaussian::read_restart(FILE *fp)
void BondGaussian::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nbondtypes; i++) {
fprintf(fp,"%d %g %d",i,bond_temperature[i],nterms[i]);
fprintf(fp, "%d %g %d", i, bond_temperature[i], nterms[i]);
for (int j = 0; j < nterms[i]; j++) {
fprintf(fp," %g %g %g",alpha[i][j],width[i][j],r0[i][j]);
fprintf(fp, " %g %g %g", alpha[i][j], width[i][j], r0[i][j]);
}
fprintf(fp, "\n");
}
@ -267,26 +268,25 @@ void BondGaussian::write_data(FILE *fp)
/* ---------------------------------------------------------------------- */
double BondGaussian::single(int type, double rsq, int /*i*/, int /*j*/,
double &fforce)
double BondGaussian::single(int type, double rsq, int /*i*/, int /*j*/, double &fforce)
{
double r = sqrt(rsq);
fforce = 0;
double sum_g_i = 0.0;
double sum_numerator = 0.0;
for (int i = 0; i < nterms[type]; i++) {
double dr = r - r0[type][i];
double prefactor = (alpha[type][i]/(width[type][i]*sqrt(MY_PI2)));
double exponent = -2*dr*dr/(width[type][i]*width[type][i]);
double g_i = prefactor*exp(exponent);
sum_g_i += g_i;
sum_numerator += g_i*dr/(width[type][i]*width[type][i]);
}
for (int i = 0; i < nterms[type]; i++) {
double dr = r - r0[type][i];
double prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2)));
double exponent = -2 * dr * dr / (width[type][i] * width[type][i]);
double g_i = prefactor * exp(exponent);
sum_g_i += g_i;
sum_numerator += g_i * dr / (width[type][i] * width[type][i]);
}
if (sum_g_i < SMALL) sum_g_i = SMALL;
if (r > 0.0) fforce = -4.0*(force->boltz*bond_temperature[type])*(sum_numerator/sum_g_i)/r;
if (r > 0.0)
fforce = -4.0 * (force->boltz * bond_temperature[type]) * (sum_numerator / sum_g_i) / r;
return -(force->boltz*bond_temperature[type])*log(sum_g_i);
return -(force->boltz * bond_temperature[type]) * log(sum_g_i);
}