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

Merge branch 'master' into fix-property-array

Browse Source
This commit is contained in:
Steve Plimpton
2021-08-17 12:19:27 -06:00
parent e6af1eecb8 527240e941
commit 96280dc185
23 changed files with 971 additions and 885 deletions
Download Patch File Download Diff File Expand all files Collapse all files

9
doc/src/Developer_utils.rst
View File

@ -334,10 +334,11 @@ arguments of commands in LAMMPS are parsed and to make abstractions of
repetitive tasks.
The :cpp:class:`LAMMPS_NS::ArgInfo` class provides an abstraction
for parsing references to compute or fix styles or variables. These
would start with a "c\_", "f\_", "v\_" followed by the ID or name of
than instance and may be postfixed with one or two array indices
"[<number>]" with numbers > 0.
for parsing references to compute or fix styles, variables or custom
integer or double properties handled by :doc:`fix property/atom <fix_property_atom>`.
These would start with a "c\_", "f\_", "v\_", "d\_", "d2\_", "i\_", or "i2_"
followed by the ID or name of than instance and may be postfixed with
one or two array indices "[<number>]" with numbers > 0.
A typical code segment would look like this:

216
src/ASPHERE/pair_resquared.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
@ -18,25 +17,24 @@
#include "pair_resquared.h"
#include <cmath>
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "force.h"
#include "math_extra.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairRESquared::PairRESquared(LAMMPS *lmp) : Pair(lmp),
cr60(pow(60.0,1.0/3.0)),
b_alpha(45.0/56.0)
PairRESquared::PairRESquared(LAMMPS *lmp) :
Pair(lmp), cr60(pow(60.0, 1.0 / 3.0)), b_alpha(45.0 / 56.0)
{
single_enable = 0;
@ -137,9 +135,10 @@ void PairRESquared::compute(int eflag, int vflag)
r6inv = r2inv * r2inv * r2inv;
forcelj = r6inv * (lj1[itype][jtype] * r6inv - lj2[itype][jtype]);
forcelj *= -r2inv;
if (eflag) one_eng =
r6inv*(r6inv*lj3[itype][jtype]-lj4[itype][jtype]) -
offset[itype][jtype];
if (eflag) {
one_eng = r6inv * (r6inv * lj3[itype][jtype] - lj4[itype][jtype]);
one_eng -= offset[itype][jtype];
}
fforce[0] = r12[0] * forcelj;
fforce[1] = r12[1] * forcelj;
fforce[2] = r12[2] * forcelj;
@ -192,8 +191,8 @@ void PairRESquared::compute(int eflag, int vflag)
if (eflag) evdwl = factor_lj * one_eng;
if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
evdwl,0.0,fforce[0],fforce[1],fforce[2],
if (evflag)
ev_tally_xyz(i, j, nlocal, newton_pair, evdwl, 0.0, fforce[0], fforce[1], fforce[2],
-r12[0], -r12[1], -r12[2]);
}
}
@ -209,31 +208,30 @@ void PairRESquared::compute(int eflag, int vflag)
void PairRESquared::allocate()
{
allocated = 1;
int n = atom->ntypes;
const int n = atom->ntypes + 1;
memory->create(setflag,n+1,n+1,"pair:setflag");
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
memory->create(setflag, n, n, "pair:setflag");
for (int i = 1; i < n; i++)
for (int j = i; j < n; j++) setflag[i][j] = 0;
memory->create(cutsq,n+1,n+1,"pair:cutsq");
memory->create(cutsq, n, n, "pair:cutsq");
memory->create(form,n+1,n+1,"pair:form");
memory->create(epsilon,n+1,n+1,"pair:epsilon");
memory->create(sigma,n+1,n+1,"pair:sigma");
memory->create(shape1,n+1,3,"pair:shape1");
memory->create(shape2,n+1,3,"pair:shape2");
memory->create(well,n+1,3,"pair:well");
memory->create(cut,n+1,n+1,"pair:cut");
memory->create(lj1,n+1,n+1,"pair:lj1");
memory->create(lj2,n+1,n+1,"pair:lj2");
memory->create(lj3,n+1,n+1,"pair:lj3");
memory->create(lj4,n+1,n+1,"pair:lj4");
memory->create(offset,n+1,n+1,"pair:offset");
memory->create(form, n, n, "pair:form");
memory->create(epsilon, n, n, "pair:epsilon");
memory->create(sigma, n, n, "pair:sigma");
memory->create(shape1, n, 3, "pair:shape1");
memory->create(shape2, n, 3, "pair:shape2");
memory->create(well, n, 3, "pair:well");
memory->create(cut, n, n, "pair:cut");
memory->create(lj1, n, n, "pair:lj1");
memory->create(lj2, n, n, "pair:lj2");
memory->create(lj3, n, n, "pair:lj3");
memory->create(lj4, n, n, "pair:lj4");
memory->create(offset, n, n, "pair:offset");
lshape = new double[n+1];
setwell = new int[n+1];
for (int i = 1; i <= n; i++) setwell[i] = 0;
lshape = new double[n];
setwell = new int[n];
for (int i = 1; i < n; i++) setwell[i] = 0;
}
/* ----------------------------------------------------------------------
@ -262,8 +260,7 @@ void PairRESquared::settings(int narg, char **arg)
void PairRESquared::coeff(int narg, char **arg)
{
if (narg < 10 || narg > 11)
error->all(FLERR,"Incorrect args for pair coefficients");
if (narg < 10 || narg > 11) error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo, ihi, jlo, jhi;
@ -292,15 +289,19 @@ void PairRESquared::coeff(int narg, char **arg)
well[i][0] = eia_one;
well[i][1] = eib_one;
well[i][2] = eic_one;
if (eia_one == 1.0 && eib_one == 1.0 && eic_one == 1.0) setwell[i] = 2;
else setwell[i] = 1;
if (eia_one == 1.0 && eib_one == 1.0 && eic_one == 1.0)
setwell[i] = 2;
else
setwell[i] = 1;
}
if (eja_one != 0.0 || ejb_one != 0.0 || ejc_one != 0.0) {
well[j][0] = eja_one;
well[j][1] = ejb_one;
well[j][2] = ejc_one;
if (eja_one == 1.0 && ejb_one == 1.0 && ejc_one == 1.0) setwell[j] = 2;
else setwell[j] = 1;
if (eja_one == 1.0 && ejb_one == 1.0 && ejc_one == 1.0)
setwell[j] = 2;
else
setwell[j] = 1;
}
setflag[i][j] = 1;
count++;
@ -346,11 +347,9 @@ double PairRESquared::init_one(int i, int j)
error->all(FLERR, "Pair resquared epsilon a,b,c coeffs are not all set");
int ishape = 0;
if (shape1[i][0] != 0.0 && shape1[i][1] != 0.0 && shape1[i][2] != 0.0)
ishape = 1;
if (shape1[i][0] != 0.0 && shape1[i][1] != 0.0 && shape1[i][2] != 0.0) ishape = 1;
int jshape = 0;
if (shape1[j][0] != 0.0 && shape1[j][1] != 0.0 && shape1[j][2] != 0.0)
jshape = 1;
if (shape1[j][0] != 0.0 && shape1[j][1] != 0.0 && shape1[j][2] != 0.0) jshape = 1;
if (ishape == 0 && jshape == 0) {
form[i][j] = SPHERE_SPHERE;
@ -371,13 +370,11 @@ double PairRESquared::init_one(int i, int j)
if (setflag[i][j] == 0) {
if (setflag[j][i] == 0) {
if (ishape == 0 && jshape == 0) {
epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j],
sigma[i][i],sigma[j][j]);
epsilon[i][j] = mix_energy(epsilon[i][i], epsilon[j][j], sigma[i][i], sigma[j][j]);
sigma[i][j] = mix_distance(sigma[i][i], sigma[j][j]);
cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
} else
error->all(FLERR,
"Pair resquared epsilon and sigma coeffs are not all set");
error->all(FLERR, "Pair resquared epsilon and sigma coeffs are not all set");
}
epsilon[i][j] = epsilon[j][i];
sigma[i][j] = sigma[j][i];
@ -392,7 +389,8 @@ double PairRESquared::init_one(int i, int j)
if (offset_flag && (cut[i][j] > 0.0)) {
double ratio = sigma[i][j] / cut[i][j];
offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio, 12.0) - pow(ratio, 6.0));
} else offset[i][j] = 0.0;
} else
offset[i][j] = 0.0;
epsilon[j][i] = epsilon[i][j];
sigma[j][i] = sigma[i][j];
@ -505,10 +503,10 @@ void PairRESquared::precompute_i(const int i,RE2Vars &ws)
MathExtra::rotation_generator_x(ws.A, ws.lA[0]);
MathExtra::rotation_generator_y(ws.A, ws.lA[1]);
MathExtra::rotation_generator_z(ws.A, ws.lA[2]);
for (int i=0; i<3; i++) {
MathExtra::times3(aTs,ws.lA[i],ws.lAtwo[i]);
MathExtra::transpose_times3(ws.lA[i],ws.sa,ws.lAsa[i]);
MathExtra::plus3(ws.lAsa[i],ws.lAtwo[i],ws.lAsa[i]);
for (int m = 0; m < 3; m++) {
MathExtra::times3(aTs, ws.lA[m], ws.lAtwo[m]);
MathExtra::transpose_times3(ws.lA[m], ws.sa, ws.lAsa[m]);
MathExtra::plus3(ws.lAsa[m], ws.lAtwo[m], ws.lAsa[m]);
}
}
@ -517,6 +515,7 @@ void PairRESquared::precompute_i(const int i,RE2Vars &ws)
derivative of m (m2)
------------------------------------------------------------------------- */
// clang-format off
double PairRESquared::det_prime(const double m[3][3], const double m2[3][3])
{
double ans;
@ -531,16 +530,15 @@ double PairRESquared::det_prime(const double m[3][3], const double m2[3][3])
m2[2][0]*m[0][1]*m[1][2] - m2[2][0]*m[0][2]*m[1][1];
return ans;
}
// clang-format on
/* ----------------------------------------------------------------------
Compute the energy, force, torque for a pair (INTEGRATED-INTEGRATED)
------------------------------------------------------------------------- */
double PairRESquared::resquared_analytic(const int i, const int j,
const RE2Vars &wi, const RE2Vars &wj,
const double *r, const double rsq,
double *fforce, double *ttor,
double *rtor)
double PairRESquared::resquared_analytic(const int i, const int j, const RE2Vars &wi,
const RE2Vars &wj, const double *r, const double rsq,
double *fforce, double *ttor, double *rtor)
{
int *type = atom->type;
@ -649,16 +647,14 @@ double PairRESquared::resquared_analytic(const int i, const int j,
tprod = eta * chi * sigh;
double stemp = h12 / 2.0;
Ua = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)*
(shape1[type[i]][2]+stemp)*(shape1[type[j]][0]+stemp)*
(shape1[type[j]][1]+stemp)*(shape1[type[j]][2]+stemp);
Ua = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) *
(shape1[type[j]][0] + stemp) * (shape1[type[j]][1] + stemp) * (shape1[type[j]][2] + stemp);
Ua = (1.0 + 3.0 * tprod) * sprod / Ua;
Ua = epsilon[type[i]][type[j]] * Ua / -36.0;
stemp = h12 / cr60;
Ur = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)*
(shape1[type[i]][2]+stemp)*(shape1[type[j]][0]+stemp)*
(shape1[type[j]][1]+stemp)*(shape1[type[j]][2]+stemp);
Ur = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) *
(shape1[type[j]][0] + stemp) * (shape1[type[j]][1] + stemp) * (shape1[type[j]][2] + stemp);
Ur = (1.0 + b_alpha * tprod) * sprod / Ur;
Ur = epsilon[type[i]][type[j]] * Ur * pow(sigh, 6.0) / 2025.0;
@ -704,31 +700,25 @@ double PairRESquared::resquared_analytic(const int i, const int j,
spr[1] = 0.5 * sigma12p3 * s[1];
spr[2] = 0.5 * sigma12p3 * s[2];
stemp = 1.0/(shape1[type[i]][0]*2.0+h12)+
1.0/(shape1[type[i]][1]*2.0+h12)+
1.0/(shape1[type[i]][2]*2.0+h12)+
1.0/(shape1[type[j]][0]*2.0+h12)+
1.0/(shape1[type[j]][1]*2.0+h12)+
1.0/(shape1[type[j]][2]*2.0+h12);
stemp = 1.0 / (shape1[type[i]][0] * 2.0 + h12) + 1.0 / (shape1[type[i]][1] * 2.0 + h12) +
1.0 / (shape1[type[i]][2] * 2.0 + h12) + 1.0 / (shape1[type[j]][0] * 2.0 + h12) +
1.0 / (shape1[type[j]][1] * 2.0 + h12) + 1.0 / (shape1[type[j]][2] * 2.0 + h12);
hsec = h12 + 3.0 * sec;
dspu = 1.0 / h12 - 1.0 / hsec + stemp;
pbsu = 3.0 * sigma[type[i]][type[j]] / hsec;
stemp = 1.0/(shape1[type[i]][0]*cr60+h12)+
1.0/(shape1[type[i]][1]*cr60+h12)+
1.0/(shape1[type[i]][2]*cr60+h12)+
1.0/(shape1[type[j]][0]*cr60+h12)+
1.0/(shape1[type[j]][1]*cr60+h12)+
1.0/(shape1[type[j]][2]*cr60+h12);
stemp = 1.0 / (shape1[type[i]][0] * cr60 + h12) + 1.0 / (shape1[type[i]][1] * cr60 + h12) +
1.0 / (shape1[type[i]][2] * cr60 + h12) + 1.0 / (shape1[type[j]][0] * cr60 + h12) +
1.0 / (shape1[type[j]][1] * cr60 + h12) + 1.0 / (shape1[type[j]][2] * cr60 + h12);
hsec = h12 + b_alpha * sec;
dspr = 7.0 / h12 - 1.0 / hsec + stemp;
pbsr = b_alpha * sigma[type[i]][type[j]] / hsec;
for (int i=0; i<3; i++) {
u[0] = -rhat[i]*rhat[0];
u[1] = -rhat[i]*rhat[1];
u[2] = -rhat[i]*rhat[2];
u[i] += 1.0;
for (int m = 0; m < 3; m++) {
u[0] = -rhat[m] * rhat[0];
u[1] = -rhat[m] * rhat[1];
u[2] = -rhat[m] * rhat[2];
u[m] += 1.0;
u[0] /= rnorm;
u[1] /= rnorm;
u[2] /= rnorm;
@ -750,10 +740,10 @@ double PairRESquared::resquared_analytic(const int i, const int j,
deta -= ddH * tdH;
deta -= dsigma1 * teta1 + dsigma2 * teta2;
dchi = MathExtra::dot3(u, fourw);
dh12 = rhat[i]+MathExtra::dot3(u,spr);
dh12 = rhat[m] + MathExtra::dot3(u, spr);
dUa = pbsu * (eta * dchi + deta * chi) - dh12 * dspu;
dUr = pbsr * (eta * dchi + deta * chi) - dh12 * dspr;
fforce[i]=dUr*Ur+dUa*Ua;
fforce[m] = dUr * Ur + dUa * Ua;
}
// torque on i
@ -788,8 +778,7 @@ double PairRESquared::resquared_analytic(const int i, const int j,
// torque on j
if (!(force->newton_pair || j < atom->nlocal))
return Ua+Ur;
if (!(force->newton_pair || j < atom->nlocal)) return Ua + Ur;
MathExtra::vecmat(fourw, wj.aTe, fwae);
@ -826,10 +815,8 @@ double PairRESquared::resquared_analytic(const int i, const int j,
Compute the energy, force, torque for a pair (INTEGRATED-LJ)
------------------------------------------------------------------------- */
double PairRESquared::resquared_lj(const int i, const int j,
const RE2Vars &wi, const double *r,
const double rsq, double *fforce,
double *ttor, bool calc_torque)
double PairRESquared::resquared_lj(const int i, const int j, const RE2Vars &wi, const double *r,
const double rsq, double *fforce, double *ttor, bool calc_torque)
{
int *type = atom->type;
@ -882,8 +869,7 @@ double PairRESquared::resquared_lj(const int i, const int j,
scorrect[2] = scorrect[2] * scorrect[2] / 2.0;
MathExtra::transpose_diag3(wi.A, scorrect, aTs);
MathExtra::times3(aTs, wi.A, gamma);
for (int ii=0; ii<3; ii++)
MathExtra::times3(aTs,wi.lA[ii],lAtwo[ii]);
for (int ii = 0; ii < 3; ii++) MathExtra::times3(aTs, wi.lA[ii], lAtwo[ii]);
rnorm = sqrt(rsq);
rhat[0] = r[0] / rnorm;
@ -912,14 +898,14 @@ double PairRESquared::resquared_lj(const int i, const int j,
h12p3 = pow(h12, 3.0);
double sigmap3 = pow(sigma[type[i]][type[j]], 3.0);
double stemp = h12 / 2.0;
Ua = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)*
(shape1[type[i]][2]+stemp)*h12p3/8.0;
Ua = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) *
h12p3 / 8.0;
Ua = (1.0 + 3.0 * tprod) * lshape[type[i]] / Ua;
Ua = epsilon[type[i]][type[j]] * Ua * sigmap3 * solv_f_a;
stemp = h12 / cr60;
Ur = (shape1[type[i]][0]+stemp)*(shape1[type[i]][1]+stemp)*
(shape1[type[i]][2]+stemp)*h12p3/60.0;
Ur = (shape1[type[i]][0] + stemp) * (shape1[type[i]][1] + stemp) * (shape1[type[i]][2] + stemp) *
h12p3 / 60.0;
Ur = (1.0 + b_alpha * tprod) * lshape[type[i]] / Ur;
Ur = epsilon[type[i]][type[j]] * Ur * sigmap3 * pow(sigh, 6.0) * solv_f_r;
@ -934,35 +920,31 @@ double PairRESquared::resquared_lj(const int i, const int j,
spr[1] = 0.5 * sigma12p3 * s[1];
spr[2] = 0.5 * sigma12p3 * s[2];
stemp = 1.0/(shape1[type[i]][0]*2.0+h12)+
1.0/(shape1[type[i]][1]*2.0+h12)+
1.0/(shape1[type[i]][2]*2.0+h12)+
3.0/h12;
stemp = 1.0 / (shape1[type[i]][0] * 2.0 + h12) + 1.0 / (shape1[type[i]][1] * 2.0 + h12) +
1.0 / (shape1[type[i]][2] * 2.0 + h12) + 3.0 / h12;
hsec = h12 + 3.0 * sec;
dspu = 1.0 / h12 - 1.0 / hsec + stemp;
pbsu = 3.0 * sigma[type[i]][type[j]] / hsec;
stemp = 1.0/(shape1[type[i]][0]*cr60+h12)+
1.0/(shape1[type[i]][1]*cr60+h12)+
1.0/(shape1[type[i]][2]*cr60+h12)+
3.0/h12;
stemp = 1.0 / (shape1[type[i]][0] * cr60 + h12) + 1.0 / (shape1[type[i]][1] * cr60 + h12) +
1.0 / (shape1[type[i]][2] * cr60 + h12) + 3.0 / h12;
hsec = h12 + b_alpha * sec;
dspr = 7.0 / h12 - 1.0 / hsec + stemp;
pbsr = b_alpha * sigma[type[i]][type[j]] / hsec;
for (int i=0; i<3; i++) {
u[0] = -rhat[i]*rhat[0];
u[1] = -rhat[i]*rhat[1];
u[2] = -rhat[i]*rhat[2];
u[i] += 1.0;
for (int m = 0; m < 3; m++) {
u[0] = -rhat[m] * rhat[0];
u[1] = -rhat[m] * rhat[1];
u[2] = -rhat[m] * rhat[2];
u[m] += 1.0;
u[0] /= rnorm;
u[1] /= rnorm;
u[2] /= rnorm;
dchi = MathExtra::dot3(u, fourw);
dh12 = rhat[i]+MathExtra::dot3(u,spr);
dh12 = rhat[m] + MathExtra::dot3(u, spr);
dUa = pbsu * dchi - dh12 * dspu;
dUr = pbsr * dchi - dh12 * dspr;
fforce[i]=dUr*Ur+dUa*Ua;
fforce[m] = dUr * Ur + dUa * Ua;
}
// torque on i
@ -970,17 +952,17 @@ double PairRESquared::resquared_lj(const int i, const int j,
if (calc_torque) {
MathExtra::vecmat(fourw, wi.aTe, fwae);
for (int i=0; i<3; i++) {
MathExtra::matvec(wi.lA[i],rhat,p);
for (int m = 0; m < 3; m++) {
MathExtra::matvec(wi.lA[m], rhat, p);
double tempv[3];
MathExtra::matvec(wi.lA[i],w,tempv);
MathExtra::matvec(wi.lA[m], w, tempv);
dchi = -MathExtra::dot3(fwae, tempv);
MathExtra::matvec(lAtwo[i],spr,tempv);
MathExtra::matvec(lAtwo[m], spr, tempv);
dh12 = -MathExtra::dot3(s, tempv);
dUa = pbsu * dchi - dh12 * dspu;
dUr = pbsr * dchi - dh12 * dspr;
ttor[i] = -(dUa*Ua+dUr*Ur);
ttor[m] = -(dUa * Ua + dUr * Ur);
}
}

18
src/BODY/pair_body_rounded_polygon.cpp
View File

@ -21,21 +21,21 @@
#include "pair_body_rounded_polygon.h"
#include <cmath>
#include <cstring>
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_body.h"
#include "body_rounded_polygon.h"
#include "comm.h"
#include "force.h"
#include "fix.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "math_extra.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;

21
src/BODY/pair_body_rounded_polyhedron.cpp
View File

@ -23,22 +23,22 @@
#include "pair_body_rounded_polyhedron.h"
#include <cmath>
#include <cstring>
#include "atom.h"
#include "atom_vec_body.h"
#include "body_rounded_polyhedron.h"
#include "comm.h"
#include "force.h"
#include "fix.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
#include "math_extra.h"
#include "fix.h"
#include "force.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "modify.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace MathConst;
@ -1212,10 +1212,7 @@ int PairBodyRoundedPolyhedron::interaction_edge_to_edge(int ibody,
contact_list[num_contacts].unique = 1;
num_contacts++;
}
} else {
}
return interact;
}

10
src/DPD-REACT/pair_exp6_rx.cpp
View File

@ -617,8 +617,7 @@ void PairExp6rx::coeff(int narg, char **arg)
if (strcmp(site1,"1fluid") == 0)
isite1 = oneFluidApproxParameter;
else
{
else {
int isp;
for (isp = 0; isp < nspecies; isp++)
if (strcmp(site1, &atom->dvname[isp][0]) == 0) break;
@ -631,8 +630,7 @@ void PairExp6rx::coeff(int narg, char **arg)
if (strcmp(site2,"1fluid") == 0)
isite2 = oneFluidApproxParameter;
else
{
else {
int isp;
for (isp = 0; isp < nspecies; isp++)
if (strcmp(site2, &atom->dvname[isp][0]) == 0) break;
@ -644,8 +642,7 @@ void PairExp6rx::coeff(int narg, char **arg)
}
// Set the interaction potential type to the enumerated type.
for (int iparam = 0; iparam < nparams; ++iparam)
{
for (int iparam = 0; iparam < nparams; ++iparam) {
if (strcmp( params[iparam].potential, "exp6") == 0)
params[iparam].potentialType = exp6PotentialType;
else
@ -663,7 +660,6 @@ void PairExp6rx::coeff(int narg, char **arg)
scalingFlag = EXPONENT;
exponentR = utils::numeric(FLERR,arg[6],false,lmp);
exponentEpsilon = utils::numeric(FLERR,arg[7],false,lmp);
if (narg > 9) error->all(FLERR,"Incorrect args for pair coefficients");
if (narg == 9) cut_one = utils::numeric(FLERR,arg[8],false,lmp);
} else if (strcmp(arg[5],"polynomial") == 0) {
scalingFlag = POLYNOMIAL;

2
src/EXTRA-FIX/fix_gld.cpp
View File

@ -132,6 +132,7 @@ FixGLD::FixGLD(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR, "Illegal fix gld command");
}
if (strcmp(arg[iarg+1],"no") == 0) {
zeroflag = 0;
} else if (strcmp(arg[iarg+1],"yes") == 0) {
zeroflag = 1;
} else {
@ -144,6 +145,7 @@ FixGLD::FixGLD(LAMMPS *lmp, int narg, char **arg) :
error->all(FLERR, "Illegal fix gld command");
}
if (strcmp(arg[iarg+1],"no") == 0) {
freezeflag = 0;
} else if (strcmp(arg[iarg+1],"yes") == 0) {
freezeflag = 1;
for (int i = 0; i < atom->nlocal; i++) {

2
src/GPU/pppm_gpu.cpp
View File

@ -106,6 +106,8 @@ PPPMGPU::PPPMGPU(LAMMPS *lmp) : PPPM(lmp)
PPPMGPU::~PPPMGPU()
{
PPPM_GPU_API(clear)(poisson_time);
destroy_3d_offset(density_brick_gpu,nzlo_out,nylo_out);
destroy_3d_offset(vd_brick,nzlo_out,nylo_out);
}
/* ----------------------------------------------------------------------

5
src/GRANULAR/fix_wall_gran.cpp
View File

@ -437,7 +437,7 @@ FixWallGran::FixWallGran(LAMMPS *lmp, int narg, char **arg) :
// register with Atom class
history_one = nullptr;
grow_arrays(atom->nmax);
FixWallGran::grow_arrays(atom->nmax);
atom->add_callback(Atom::GROW);
atom->add_callback(Atom::RESTART);
@ -1555,8 +1555,7 @@ double FixWallGran::memory_usage()
void FixWallGran::grow_arrays(int nmax)
{
if (use_history) memory->grow(history_one,nmax,size_history,
"fix_wall_gran:history_one");
if (use_history) memory->grow(history_one,nmax,size_history,"fix_wall_gran:history_one");
if (peratom_flag) {
memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom");
}

5
src/GRANULAR/fix_wall_gran_region.cpp
View File

@ -62,7 +62,7 @@ FixWallGranRegion::FixWallGranRegion(LAMMPS *lmp, int narg, char **arg) :
ncontact = nullptr;
walls = nullptr;
history_many = nullptr;
grow_arrays(atom->nmax);
FixWallGranRegion::grow_arrays(atom->nmax);
// initialize shear history as if particle is not touching region
@ -355,8 +355,7 @@ void FixWallGranRegion::grow_arrays(int nmax)
if (use_history) {
memory->grow(ncontact,nmax,"fix_wall_gran:ncontact");
memory->grow(walls,nmax,tmax,"fix_wall_gran:walls");
memory->grow(history_many,nmax,tmax,size_history,
"fix_wall_gran:history_many");
memory->grow(history_many,nmax,tmax,size_history,"fix_wall_gran:history_many");
}
if (peratom_flag)
memory->grow(array_atom,nmax,size_peratom_cols,"fix_wall_gran:array_atom");

1
src/GRANULAR/pair_granular.cpp
View File

@ -110,6 +110,7 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp)
PairGranular::~PairGranular()
{
delete[] svector;
delete[] history_transfer_factors;
if (!fix_history) modify->delete_fix("NEIGH_HISTORY_GRANULAR_DUMMY");
else modify->delete_fix("NEIGH_HISTORY_GRANULAR");

2
src/MANYBODY/pair_tersoff_table.cpp
View File

@ -91,7 +91,7 @@ PairTersoffTable::~PairTersoffTable()
memory->destroy(cutsq);
}
deallocateGrids();
deallocatePreLoops();
PairTersoffTable::deallocatePreLoops();
}
/* ---------------------------------------------------------------------- */

2
src/OPENMP/pair_tersoff_table_omp.cpp
View File

@ -53,7 +53,7 @@ PairTersoffTableOMP::PairTersoffTableOMP(LAMMPS *lmp) :
PairTersoffTableOMP::~PairTersoffTableOMP()
{
if (allocated) {
deallocatePreLoops();
PairTersoffTableOMP::deallocatePreLoops();
}
}

2
src/REAXFF/fix_qeq_reaxff.cpp
View File

@ -165,7 +165,7 @@ FixQEqReaxFF::~FixQEqReaxFF()
memory->destroy(s_hist);
memory->destroy(t_hist);
deallocate_storage();
FixQEqReaxFF::deallocate_storage();
deallocate_matrix();
memory->destroy(shld);

338
src/REPLICA/fix_pimd.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
@ -24,18 +23,18 @@
#include "fix_pimd.h"
#include <cmath>
#include <cstring>
#include "universe.h"
#include "comm.h"
#include "force.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "update.h"
#include "error.h"
#include "force.h"
#include "math_const.h"
#include "memory.h"
#include "error.h"
#include "universe.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
@ -47,46 +46,6 @@ enum{PIMD,NMPIMD,CMD};
FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
{
method = PIMD;
fmass = 1.0;
nhc_temp = 298.15;
nhc_nchain = 2;
sp = 1.0;
for (int i=3; i<narg-1; i+=2)
{
if (strcmp(arg[i],"method")==0)
{
if (strcmp(arg[i+1],"pimd")==0) method=PIMD;
else if (strcmp(arg[i+1],"nmpimd")==0) method=NMPIMD;
else if (strcmp(arg[i+1],"cmd")==0) method=CMD;
else error->universe_all(FLERR,"Unknown method parameter for fix pimd");
}
else if (strcmp(arg[i],"fmass")==0)
{
fmass = atof(arg[i+1]);
if (fmass<0.0 || fmass>1.0) error->universe_all(FLERR,"Invalid fmass value for fix pimd");
}
else if (strcmp(arg[i],"sp")==0)
{
sp = atof(arg[i+1]);
if (fmass<0.0) error->universe_all(FLERR,"Invalid sp value for fix pimd");
}
else if (strcmp(arg[i],"temp")==0)
{
nhc_temp = atof(arg[i+1]);
if (nhc_temp<0.0) error->universe_all(FLERR,"Invalid temp value for fix pimd");
}
else if (strcmp(arg[i],"nhc")==0)
{
nhc_nchain = atoi(arg[i+1]);
if (nhc_nchain<2) error->universe_all(FLERR,"Invalid nhc value for fix pimd");
}
else error->universe_all(arg[i],i+1,"Unknown keyword for fix pimd");
}
/* Initiation */
max_nsend = 0;
tag_send = nullptr;
buf_send = nullptr;
@ -110,6 +69,41 @@ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
nhc_eta_dotdot = nullptr;
nhc_eta_mass = nullptr;
method = PIMD;
fmass = 1.0;
nhc_temp = 298.15;
nhc_nchain = 2;
sp = 1.0;
for (int i = 3; i < narg - 1; i += 2) {
if (strcmp(arg[i], "method") == 0) {
if (strcmp(arg[i + 1], "pimd") == 0)
method = PIMD;
else if (strcmp(arg[i + 1], "nmpimd") == 0)
method = NMPIMD;
else if (strcmp(arg[i + 1], "cmd") == 0)
method = CMD;
else
error->universe_all(FLERR, "Unknown method parameter for fix pimd");
} else if (strcmp(arg[i], "fmass") == 0) {
fmass = utils::numeric(FLERR, arg[i + 1], false, lmp);
if (fmass < 0.0 || fmass > 1.0)
error->universe_all(FLERR, "Invalid fmass value for fix pimd");
} else if (strcmp(arg[i], "sp") == 0) {
sp = utils::numeric(FLERR, arg[i + 1], false, lmp);
if (fmass < 0.0) error->universe_all(FLERR, "Invalid sp value for fix pimd");
} else if (strcmp(arg[i], "temp") == 0) {
nhc_temp = utils::numeric(FLERR, arg[i + 1], false, lmp);
if (nhc_temp < 0.0) error->universe_all(FLERR, "Invalid temp value for fix pimd");
} else if (strcmp(arg[i], "nhc") == 0) {
nhc_nchain = utils::inumeric(FLERR, arg[i + 1], false, lmp);
if (nhc_nchain < 2) error->universe_all(FLERR, "Invalid nhc value for fix pimd");
} else
error->universe_all(FLERR, fmt::format("Unknown keyword {} for fix pimd", arg[i]));
}
/* Initiation */
size_peratom_cols = 12 * nhc_nchain + 3;
nhc_offset_one_1 = 3 * nhc_nchain;
@ -138,7 +132,37 @@ FixPIMD::FixPIMD(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
}
/* ---------------------------------------------------------------------- */
FixPIMD::~FixPIMD()
{
delete[] mass;
atom->delete_callback(id, Atom::GROW);
atom->delete_callback(id, Atom::RESTART);
memory->destroy(M_x2xp);
memory->destroy(M_xp2x);
memory->destroy(M_f2fp);
memory->destroy(M_fp2f);
memory->sfree(lam);
if (buf_beads)
for (int i = 0; i < np; i++) memory->sfree(buf_beads[i]);
delete[] buf_beads;
delete[] plan_send;
delete[] plan_recv;
delete[] mode_index;
memory->sfree(tag_send);
memory->sfree(buf_send);
memory->sfree(buf_recv);
memory->destroy(array_atom);
memory->destroy(nhc_eta);
memory->destroy(nhc_eta_dot);
memory->destroy(nhc_eta_dotdot);
memory->destroy(nhc_eta_mass);
}
/* ---------------------------------------------------------------------- */
int FixPIMD::setmask()
{
int mask = 0;
@ -155,7 +179,8 @@ void FixPIMD::init()
if (atom->map_style == Atom::MAP_NONE)
error->all(FLERR, "Fix pimd requires an atom map, see atom_modify");
if (universe->me==0 && screen) fprintf(screen,"Fix pimd initializing Path-Integral ...\n");
if (universe->me == 0 && universe->uscreen)
fprintf(universe->uscreen, "Fix pimd initializing Path-Integral ...\n");
// prepare the constants
@ -200,8 +225,10 @@ void FixPIMD::init()
mass = new double[atom->ntypes + 1];
if (method==CMD || method==NMPIMD) nmpimd_init();
else for (int i=1; i<=atom->ntypes; i++) mass[i] = atom->mass[i] / np * fmass;
if (method == CMD || method == NMPIMD)
nmpimd_init();
else
for (int i = 1; i <= atom->ntypes; i++) mass[i] = atom->mass[i] / np * fmass;
if (!nhc_ready) nhc_init();
}
@ -210,7 +237,8 @@ void FixPIMD::init()
void FixPIMD::setup(int vflag)
{
if (universe->me==0 && screen) fprintf(screen,"Setting up Path-Integral ...\n");
if (universe->me == 0 && universe->uscreen)
fprintf(universe->uscreen, "Setting up Path-Integral ...\n");
post_force(vflag);
}
@ -234,13 +262,13 @@ void FixPIMD::final_integrate()
void FixPIMD::post_force(int /*flag*/)
{
for (int i=0; i<atom->nlocal; i++) for(int j=0; j<3; j++) atom->f[i][j] /= np;
for (int i = 0; i < atom->nlocal; i++)
for (int j = 0; j < 3; j++) atom->f[i][j] /= np;
comm_exec(atom->x);
spring_force();
if (method==CMD || method==NMPIMD)
{
if (method == CMD || method == NMPIMD) {
/* forward comm for the force on ghost atoms */
nmpimd_fill(atom->f);
@ -267,28 +295,32 @@ void FixPIMD::nhc_init()
double mass0 = KT * tau * tau;
int max = 3 * atom->nlocal;
for (int i=0; i<max; i++)
{
for (int ichain=0; ichain<nhc_nchain; ichain++)
{
for (int i = 0; i < max; i++) {
for (int ichain = 0; ichain < nhc_nchain; ichain++) {
nhc_eta[i][ichain] = 0.0;
nhc_eta_dot[i][ichain] = 0.0;
nhc_eta_dot[i][ichain] = 0.0;
nhc_eta_dotdot[i][ichain] = 0.0;
nhc_eta_mass[i][ichain] = mass0;
if ((method==CMD || method==NMPIMD) && universe->iworld==0) ; else nhc_eta_mass[i][ichain] *= fmass;
if ((method == CMD || method == NMPIMD) && universe->iworld == 0)
;
else
nhc_eta_mass[i][ichain] *= fmass;
}
nhc_eta_dot[i][nhc_nchain] = 0.0;
for (int ichain = 1; ichain < nhc_nchain; ichain++)
nhc_eta_dotdot[i][ichain] = (nhc_eta_mass[i][ichain-1] * nhc_eta_dot[i][ichain-1]
* nhc_eta_dot[i][ichain-1] * force->mvv2e - KT) / nhc_eta_mass[i][ichain];
nhc_eta_dotdot[i][ichain] = (nhc_eta_mass[i][ichain - 1] * nhc_eta_dot[i][ichain - 1] *
nhc_eta_dot[i][ichain - 1] * force->mvv2e -
KT) /
nhc_eta_mass[i][ichain];
}
// Zero NH acceleration for CMD
if (method==CMD && universe->iworld==0) for (int i=0; i<max; i++)
if (method == CMD && universe->iworld == 0)
for (int i = 0; i < max; i++)
for (int ichain = 0; ichain < nhc_nchain; ichain++) nhc_eta_dotdot[i][ichain] = 0.0;
nhc_ready = true;
@ -302,8 +334,7 @@ void FixPIMD::nhc_update_x()
double **x = atom->x;
double **v = atom->v;
if (method==CMD || method==NMPIMD)
{
if (method == CMD || method == NMPIMD) {
nmpimd_fill(atom->v);
comm_exec(atom->v);
@ -313,8 +344,7 @@ void FixPIMD::nhc_update_x()
nmpimd_transform(buf_beads, v, M_xp2x[universe->iworld]);
}
for (int i=0; i<n; i++)
{
for (int i = 0; i < n; i++) {
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
@ -330,8 +360,7 @@ void FixPIMD::nhc_update_v()
double **v = atom->v;
double **f = atom->f;
for (int i=0; i<n; i++)
{
for (int i = 0; i < n; i++) {
double dtfm = dtf / mass[type[i]];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
@ -350,8 +379,7 @@ void FixPIMD::nhc_update_v()
double dt4 = 0.25 * update->dt;
double dt8 = 0.125 * update->dt;
for (int i=0; i<nmax; i++)
{
for (int i = 0; i < nmax; i++) {
int iatm = i / 3;
int idim = i % 3;
@ -366,8 +394,7 @@ void FixPIMD::nhc_update_v()
eta_dotdot[0] = (kecurrent - KT) / nhc_eta_mass[i][0];
for (int ichain=nhc_nchain-1; ichain>0; ichain--)
{
for (int ichain = nhc_nchain - 1; ichain > 0; ichain--) {
expfac = exp(-dt8 * eta_dot[ichain + 1]);
eta_dot[ichain] *= expfac;
eta_dot[ichain] += eta_dotdot[ichain] * dt4;
@ -388,19 +415,18 @@ void FixPIMD::nhc_update_v()
kecurrent = force->boltz * t_current;
eta_dotdot[0] = (kecurrent - KT) / nhc_eta_mass[i][0];
for (int ichain=0; ichain<nhc_nchain; ichain++)
eta[ichain] += dthalf * eta_dot[ichain];
for (int ichain = 0; ichain < nhc_nchain; ichain++) eta[ichain] += dthalf * eta_dot[ichain];
eta_dot[0] *= expfac;
eta_dot[0] += eta_dotdot[0] * dt4;
eta_dot[0] *= expfac;
for (int ichain=1; ichain<nhc_nchain; ichain++)
{
for (int ichain = 1; ichain < nhc_nchain; ichain++) {
expfac = exp(-dt8 * eta_dot[ichain + 1]);
eta_dot[ichain] *= expfac;
eta_dotdot[ichain] = (nhc_eta_mass[i][ichain-1] * eta_dot[ichain-1] * eta_dot[ichain-1]
- KT) / nhc_eta_mass[i][ichain];
eta_dotdot[ichain] =
(nhc_eta_mass[i][ichain - 1] * eta_dot[ichain - 1] * eta_dot[ichain - 1] - KT) /
nhc_eta_mass[i][ichain];
eta_dot[ichain] += eta_dotdot[ichain] * dt4;
eta_dot[ichain] *= expfac;
}
@ -429,23 +455,21 @@ void FixPIMD::nmpimd_init()
lam[0] = 0.0;
if (np % 2 == 0) lam[np - 1] = 4.0 * np;
for (int i=2; i<=np/2; i++)
{
for (int i = 2; i <= np / 2; i++) {
lam[2 * i - 3] = lam[2 * i - 2] = 2.0 * np * (1.0 - 1.0 * cos(2.0 * MY_PI * (i - 1) / np));
}
// Set up eigenvectors for non-degenerated modes
for (int i=0; i<np; i++)
{
for (int i = 0; i < np; i++) {
M_x2xp[0][i] = 1.0 / np;
if (np % 2 == 0) M_x2xp[np - 1][i] = 1.0 / np * pow(-1.0, i);
}
// Set up eigenvectors for degenerated modes
for (int i=0; i<(np-1)/2; i++) for(int j=0; j<np; j++)
{
for (int i = 0; i < (np - 1) / 2; i++)
for (int j = 0; j < np; j++) {
M_x2xp[2 * i + 1][j] = sqrt(2.0) * cos(2.0 * MY_PI * (i + 1) * j / np) / np;
M_x2xp[2 * i + 2][j] = -sqrt(2.0) * sin(2.0 * MY_PI * (i + 1) * j / np) / np;
}
@ -453,8 +477,7 @@ void FixPIMD::nmpimd_init()
// Set up Ut
for (int i = 0; i < np; i++)
for (int j=0; j<np; j++)
{
for (int j = 0; j < np; j++) {
M_xp2x[i][j] = M_x2xp[j][i] * np;
M_f2fp[i][j] = M_x2xp[i][j] * np;
M_fp2f[i][j] = M_xp2x[i][j];
@ -464,12 +487,10 @@ void FixPIMD::nmpimd_init()
int iworld = universe->iworld;
for (int i=1; i<=atom->ntypes; i++)
{
for (int i = 1; i <= atom->ntypes; i++) {
mass[i] = atom->mass[i];
if (iworld)
{
if (iworld) {
mass[i] *= lam[iworld];
mass[i] *= fmass;
}
@ -491,8 +512,8 @@ void FixPIMD::nmpimd_transform(double** src, double** des, double *vector)
int n = atom->nlocal;
int m = 0;
for (int i=0; i<n; i++) for(int d=0; d<3; d++)
{
for (int i = 0; i < n; i++)
for (int d = 0; d < 3; d++) {
des[i][d] = 0.0;
for (int j = 0; j < np; j++) { des[i][d] += (src[j][m] * vector[j]); }
m++;
@ -514,8 +535,7 @@ void FixPIMD::spring_force()
double *xlast = buf_beads[x_last];
double *xnext = buf_beads[x_next];
for (int i=0; i<nlocal; i++)
{
for (int i = 0; i < nlocal; i++) {
double delx1 = xlast[0] - x[i][0];
double dely1 = xlast[1] - x[i][1];
double delz1 = xlast[2] - x[i][2];
@ -548,14 +568,12 @@ void FixPIMD::spring_force()
void FixPIMD::comm_init()
{
if (size_plan)
{
if (size_plan) {
delete[] plan_send;
delete[] plan_recv;
}
if (method == PIMD)
{
if (method == PIMD) {
size_plan = 2;
plan_send = new int[2];
plan_recv = new int[2];
@ -566,21 +584,22 @@ void FixPIMD::comm_init()
if (rank_last < 0) rank_last += universe->nprocs;
if (rank_next >= universe->nprocs) rank_next -= universe->nprocs;
plan_send[0] = rank_next; plan_send[1] = rank_last;
plan_recv[0] = rank_last; plan_recv[1] = rank_next;
plan_send[0] = rank_next;
plan_send[1] = rank_last;
plan_recv[0] = rank_last;
plan_recv[1] = rank_next;
mode_index[0] = 0; mode_index[1] = 1;
x_last = 1; x_next = 0;
}
else
{
mode_index[0] = 0;
mode_index[1] = 1;
x_last = 1;
x_next = 0;
} else {
size_plan = np - 1;
plan_send = new int[size_plan];
plan_recv = new int[size_plan];
mode_index = new int[size_plan];
for (int i=0; i<size_plan; i++)
{
for (int i = 0; i < size_plan; i++) {
plan_send[i] = universe->me + comm->nprocs * (i + 1);
if (plan_send[i] >= universe->nprocs) plan_send[i] -= universe->nprocs;
@ -594,9 +613,9 @@ void FixPIMD::comm_init()
x_last = (universe->iworld - 1 + universe->nworlds) % (universe->nworlds);
}
if (buf_beads)
{
for (int i=0; i<np; i++) if (buf_beads[i]) delete [] buf_beads[i];
if (buf_beads) {
for (int i = 0; i < np; i++)
if (buf_beads[i]) delete[] buf_beads[i];
delete[] buf_beads;
}
@ -610,8 +629,7 @@ void FixPIMD::comm_exec(double **ptr)
{
int nlocal = atom->nlocal;
if (nlocal > max_nlocal)
{
if (nlocal > max_nlocal) {
max_nlocal = nlocal + 200;
int size = sizeof(double) * max_nlocal * 3;
buf_recv = (double *) memory->srealloc(buf_recv, size, "FixPIMD:x_recv");
@ -626,47 +644,46 @@ void FixPIMD::comm_exec(double **ptr)
// go over comm plans
for (int iplan = 0; iplan<size_plan; iplan++)
{
for (int iplan = 0; iplan < size_plan; iplan++) {
// sendrecv nlocal
int nsend;
MPI_Sendrecv( &(nlocal), 1, MPI_INT, plan_send[iplan], 0,
&(nsend), 1, MPI_INT, plan_recv[iplan], 0, universe->uworld, MPI_STATUS_IGNORE);
MPI_Sendrecv(&(nlocal), 1, MPI_INT, plan_send[iplan], 0, &(nsend), 1, MPI_INT, plan_recv[iplan],
0, universe->uworld, MPI_STATUS_IGNORE);
// allocate arrays
if (nsend > max_nsend)
{
if (nsend > max_nsend) {
max_nsend = nsend + 200;
tag_send = (tagint*) memory->srealloc(tag_send, sizeof(tagint)*max_nsend, "FixPIMD:tag_send");
buf_send = (double*) memory->srealloc(buf_send, sizeof(double)*max_nsend*3, "FixPIMD:x_send");
tag_send =
(tagint *) memory->srealloc(tag_send, sizeof(tagint) * max_nsend, "FixPIMD:tag_send");
buf_send =
(double *) memory->srealloc(buf_send, sizeof(double) * max_nsend * 3, "FixPIMD:x_send");
}
// send tags
MPI_Sendrecv( atom->tag, nlocal, MPI_LMP_TAGINT, plan_send[iplan], 0,
tag_send, nsend, MPI_LMP_TAGINT, plan_recv[iplan], 0, universe->uworld, MPI_STATUS_IGNORE);
MPI_Sendrecv(atom->tag, nlocal, MPI_LMP_TAGINT, plan_send[iplan], 0, tag_send, nsend,
MPI_LMP_TAGINT, plan_recv[iplan], 0, universe->uworld, MPI_STATUS_IGNORE);
// wrap positions
double *wrap_ptr = buf_send;
int ncpy = sizeof(double) * 3;
for (int i=0; i<nsend; i++)
{
for (int i = 0; i < nsend; i++) {
int index = atom->map(tag_send[i]);
if (index<0)
{
if (index < 0) {
char error_line[256];
sprintf(error_line, "Atom " TAGINT_FORMAT " is missing at world [%d] "
"rank [%d] required by rank [%d] (" TAGINT_FORMAT ", "
TAGINT_FORMAT ", " TAGINT_FORMAT ").\n", tag_send[i],
universe->iworld, comm->me, plan_recv[iplan],
atom->tag[0], atom->tag[1], atom->tag[2]);
sprintf(error_line,
"Atom " TAGINT_FORMAT " is missing at world [%d] "
"rank [%d] required by rank [%d] (" TAGINT_FORMAT ", " TAGINT_FORMAT
", " TAGINT_FORMAT ").\n",
tag_send[i], universe->iworld, comm->me, plan_recv[iplan], atom->tag[0],
atom->tag[1], atom->tag[2]);
error->universe_one(FLERR, error_line);
}
@ -677,8 +694,8 @@ void FixPIMD::comm_exec(double **ptr)
// sendrecv x
MPI_Sendrecv( buf_send, nsend*3, MPI_DOUBLE, plan_recv[iplan], 0,
buf_recv, nlocal*3, MPI_DOUBLE, plan_send[iplan], 0, universe->uworld, MPI_STATUS_IGNORE);
MPI_Sendrecv(buf_send, nsend * 3, MPI_DOUBLE, plan_recv[iplan], 0, buf_recv, nlocal * 3,
MPI_DOUBLE, plan_send[iplan], 0, universe->uworld, MPI_STATUS_IGNORE);
// copy x
@ -688,8 +705,7 @@ void FixPIMD::comm_exec(double **ptr)
/* ---------------------------------------------------------------------- */
int FixPIMD::pack_forward_comm(int n, int *list, double *buf,
int /*pbc_flag*/, int * /*pbc*/)
int FixPIMD::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, int * /*pbc*/)
{
int i, j, m;
@ -763,10 +779,14 @@ int FixPIMD::pack_exchange(int i, double *buf)
int offset = 0;
int pos = i * 3;
memcpy(buf+offset, nhc_eta[pos], nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(buf+offset, nhc_eta_dot[pos], nhc_size_one_2); offset += nhc_offset_one_2;
memcpy(buf+offset, nhc_eta_dotdot[pos], nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(buf+offset, nhc_eta_mass[pos], nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(buf + offset, nhc_eta[pos], nhc_size_one_1);
offset += nhc_offset_one_1;
memcpy(buf + offset, nhc_eta_dot[pos], nhc_size_one_2);
offset += nhc_offset_one_2;
memcpy(buf + offset, nhc_eta_dotdot[pos], nhc_size_one_1);
offset += nhc_offset_one_1;
memcpy(buf + offset, nhc_eta_mass[pos], nhc_size_one_1);
offset += nhc_offset_one_1;
return size_peratom_cols;
}
@ -778,10 +798,14 @@ int FixPIMD::unpack_exchange(int nlocal, double *buf)
int offset = 0;
int pos = nlocal * 3;
memcpy(nhc_eta[pos], buf+offset, nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(nhc_eta_dot[pos], buf+offset, nhc_size_one_2); offset += nhc_offset_one_2;
memcpy(nhc_eta_dotdot[pos], buf+offset, nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(nhc_eta_mass[pos], buf+offset, nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(nhc_eta[pos], buf + offset, nhc_size_one_1);
offset += nhc_offset_one_1;
memcpy(nhc_eta_dot[pos], buf + offset, nhc_size_one_2);
offset += nhc_offset_one_2;
memcpy(nhc_eta_dotdot[pos], buf + offset, nhc_size_one_1);
offset += nhc_offset_one_1;
memcpy(nhc_eta_mass[pos], buf + offset, nhc_size_one_1);
offset += nhc_offset_one_1;
return size_peratom_cols;
}
@ -795,10 +819,14 @@ int FixPIMD::pack_restart(int i, double *buf)
// pack buf[0] this way because other fixes unpack it
buf[offset++] = size_peratom_cols + 1;
memcpy(buf+offset, nhc_eta[pos], nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(buf+offset, nhc_eta_dot[pos], nhc_size_one_2); offset += nhc_offset_one_2;
memcpy(buf+offset, nhc_eta_dotdot[pos], nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(buf+offset, nhc_eta_mass[pos], nhc_size_one_1); offset += nhc_offset_one_1;
memcpy(buf + offset, nhc_eta[pos], nhc_size_one_1);
offset += nhc_offset_one_1;
memcpy(buf + offset, nhc_eta_dot[pos], nhc_size_one_2);
offset += nhc_offset_one_2;
memcpy(buf + offset, nhc_eta_dotdot[pos], nhc_size_one_1);
offset += nhc_offset_one_1;
memcpy(buf + offset, nhc_eta_mass[pos], nhc_size_one_1);
offset += nhc_offset_one_1;
return size_peratom_cols + 1;
}
@ -818,10 +846,14 @@ void FixPIMD::unpack_restart(int nlocal, int nth)
int pos = nlocal * 3;
memcpy(nhc_eta[pos], extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1;
memcpy(nhc_eta_dot[pos], extra[nlocal]+m, nhc_size_one_2); m += nhc_offset_one_2;
memcpy(nhc_eta_dotdot[pos], extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1;
memcpy(nhc_eta_mass[pos], extra[nlocal]+m, nhc_size_one_1); m += nhc_offset_one_1;
memcpy(nhc_eta[pos], extra[nlocal] + m, nhc_size_one_1);
m += nhc_offset_one_1;
memcpy(nhc_eta_dot[pos], extra[nlocal] + m, nhc_size_one_2);
m += nhc_offset_one_2;
memcpy(nhc_eta_dotdot[pos], extra[nlocal] + m, nhc_size_one_1);
m += nhc_offset_one_1;
memcpy(nhc_eta_mass[pos], extra[nlocal] + m, nhc_size_one_1);
m += nhc_offset_one_1;
nhc_ready = true;
}

1
src/REPLICA/fix_pimd.h
View File

@ -27,6 +27,7 @@ namespace LAMMPS_NS {
class FixPIMD : public Fix {
public:
FixPIMD(class LAMMPS *, int, char **);
virtual ~FixPIMD();
int setmask();

6
src/REPLICA/neb.cpp
View File

@ -95,6 +95,7 @@ NEB::~NEB()
MPI_Comm_free(&roots);
memory->destroy(all);
delete[] rdist;
if (fp) fclose(fp);
}
/* ----------------------------------------------------------------------
@ -379,8 +380,8 @@ void NEB::readfile(char *file, int flag)
char line[MAXLINE];
double xx,yy,zz,delx,dely,delz;
if (me_universe == 0 && screen)
fprintf(screen,"Reading NEB coordinate file(s) ...\n");
if (me_universe == 0 && universe->uscreen)
fprintf(universe->uscreen,"Reading NEB coordinate file(s) ...\n");
// flag = 0, universe root reads header of file, bcast to universe
// flag = 1, each replica's root reads header of file, bcast to world
@ -534,6 +535,7 @@ void NEB::readfile(char *file, int flag)
else fclose(fp);
}
}
fp = nullptr;
}
/* ----------------------------------------------------------------------

6
src/SPIN/neb_spin.cpp
View File

@ -79,6 +79,7 @@ NEBSpin::~NEBSpin()
MPI_Comm_free(&roots);
memory->destroy(all);
delete[] rdist;
if (fp) fclose(fp);
}
/* ----------------------------------------------------------------------
@ -373,8 +374,8 @@ void NEBSpin::readfile(char *file, int flag)
double xx,yy,zz;
double musp,spx,spy,spz;
if (me_universe == 0 && screen)
fprintf(screen,"Reading NEBSpin coordinate file(s) ...\n");
if (me_universe == 0 && universe->uscreen)
fprintf(universe->uscreen,"Reading NEBSpin coordinate file(s) ...\n");
// flag = 0, universe root reads header of file, bcast to universe
// flag = 1, each replica's root reads header of file, bcast to world
@ -560,6 +561,7 @@ void NEBSpin::readfile(char *file, int flag)
else fclose(fp);
}
}
fp = nullptr;
}
/* ----------------------------------------------------------------------

10
src/arg_info.cpp
View File

@ -30,7 +30,8 @@ using namespace LAMMPS_NS;
ArgInfo::ArgInfo(const std::string &arg, int allowed) : type(NONE), dim(0), index1(-1), index2(-1)
{
if ((arg.size() > 2) && (arg[1] == '_')) {
if (((arg.size() > 3) && (arg[1] == '2') && (arg[2] == '_'))
|| ((arg.size() > 2) && (arg[1] == '_'))) {
if ((arg[0] == 'c') && (allowed & COMPUTE))
type = COMPUTE;
else if ((arg[0] == 'f') && (allowed & FIX))
@ -46,10 +47,11 @@ ArgInfo::ArgInfo(const std::string &arg, int allowed) : type(NONE), dim(0), inde
name = arg;
return;
}
const int offset = (arg[1] == '_') ? 2 : 3;
std::size_t has_idx1 = arg.find('[', 2);
std::size_t has_idx1 = arg.find('[', offset);
if (has_idx1 != std::string::npos) {
name = arg.substr(2, has_idx1 - 2);
name = arg.substr(offset, has_idx1 - offset);
dim = 1;
std::size_t has_idx2 = arg.find('[', has_idx1 + 1);
@ -79,7 +81,7 @@ ArgInfo::ArgInfo(const std::string &arg, int allowed) : type(NONE), dim(0), inde
}
} else {
index1 = 0;
name = arg.substr(2);
name = arg.substr(offset);
}
} else {
index1 = 0;

5
src/compute_dipole.cpp
View File

@ -29,7 +29,7 @@ enum { MASSCENTER, GEOMCENTER };
ComputeDipole::ComputeDipole(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR, "Illegal compute com command");
if ((narg < 3) || (narg > 4)) error->all(FLERR, "Illegal compute dipole command");
scalar_flag = 1;
vector_flag = 1;
@ -39,11 +39,8 @@ ComputeDipole::ComputeDipole(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, n
vector = new double[size_vector];
vector[0] = vector[1] = vector[2] = 0.0;
usecenter = MASSCENTER;
if ((narg != 3) && (narg != 4)) error->all(FLERR, "Illegal compute dipole command");
if (narg == 4) {
if (utils::strmatch(arg[3], "^geom"))
usecenter = GEOMCENTER;

13
src/read_restart.cpp
View File

@ -44,7 +44,7 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ReadRestart::ReadRestart(LAMMPS *lmp) : Command(lmp) {}
ReadRestart::ReadRestart(LAMMPS *lmp) : Command(lmp), mpiio(nullptr) {}
/* ---------------------------------------------------------------------- */
@ -110,8 +110,7 @@ void ReadRestart::command(int narg, char **arg)
}
fp = fopen(hfile.c_str(),"rb");
if (fp == nullptr)
error->one(FLERR,"Cannot open restart file {}: {}",
hfile, utils::getsyserror());
error->one(FLERR,"Cannot open restart file {}: {}", hfile, utils::getsyserror());
}
// read magic string, endian flag, format revision
@ -336,8 +335,7 @@ void ReadRestart::command(int narg, char **arg)
procfile.replace(procfile.find("%"),1,fmt::format("{}",icluster));
fp = fopen(procfile.c_str(),"rb");
if (fp == nullptr)
error->one(FLERR,"Cannot open restart file {}: {}",
procfile, utils::getsyserror());
error->one(FLERR,"Cannot open restart file {}: {}", procfile, utils::getsyserror());
}
int procsperfile;
@ -525,8 +523,9 @@ void ReadRestart::command(int narg, char **arg)
MPI_Barrier(world);
if (comm->me == 0)
utils::logmesg(lmp," read_restart CPU = {:.3f} seconds\n",
MPI_Wtime()-time1);
utils::logmesg(lmp," read_restart CPU = {:.3f} seconds\n",MPI_Wtime()-time1);
delete mpiio;
}
/* ----------------------------------------------------------------------

31
src/utils.cpp
View File

@ -177,20 +177,23 @@ char *utils::fgets_trunc(char *buf, int size, FILE *fp)
char dummy[MAXDUMMY];
char *ptr = fgets(buf, size, fp);
// EOF
// EOF?
if (!ptr) return nullptr;
int n = strlen(buf);
// line is shorter than buffer, append newline if needed,
if (n < size - 2) {
// check the string being read in:
// - if string is shorter than the buffer make sure it has a final newline and return
// - if string is exactly the size of the buffer and has a final newline return
// - otherwise truncate with final newline and read into dummy buffer until EOF or newline is found
if (n < size - 1) {
if (buf[n - 1] != '\n') {
buf[n] = '\n';
buf[n + 1] = '\0';
}
return buf;
// line fits exactly. overwrite last but one character.
} else if (buf[n - 1] == '\n') {
return buf;
} else
buf[size - 2] = '\n';
@ -203,15 +206,15 @@ char *utils::fgets_trunc(char *buf, int size, FILE *fp)
n = 0;
} while (n == MAXDUMMY - 1 && ptr[MAXDUMMY - 1] != '\n');
// return first chunk
// return truncated chunk
return buf;
}
#define MAXPATHLENBUF 1024
/* like fgets() but aborts with an error or EOF is encountered */
void utils::sfgets(const char *srcname, int srcline, char *s, int size, FILE *fp,
const char *filename, Error *error)
{
constexpr int MAXPATHLENBUF=1024;
char *rv = fgets(s, size, fp);
if (rv == nullptr) { // something went wrong
char buf[MAXPATHLENBUF];
@ -240,6 +243,7 @@ void utils::sfgets(const char *srcname, int srcline, char *s, int size, FILE *fp
void utils::sfread(const char *srcname, int srcline, void *s, size_t size, size_t num, FILE *fp,
const char *filename, Error *error)
{
constexpr int MAXPATHLENBUF=1024;
size_t rv = fread(s, size, num, fp);
if (rv != num) { // something went wrong
char buf[MAXPATHLENBUF];
@ -1093,11 +1097,6 @@ bool utils::file_is_readable(const std::string &path)
search current directory and the LAMMPS_POTENTIALS directory if
specified
------------------------------------------------------------------------- */
#if defined(_WIN32)
#define OS_PATH_VAR_SEP ";"
#else
#define OS_PATH_VAR_SEP ":"
#endif
std::string utils::get_potential_file_path(const std::string &path)
{
@ -1111,8 +1110,11 @@ std::string utils::get_potential_file_path(const std::string &path)
const char *var = getenv("LAMMPS_POTENTIALS");
if (var != nullptr) {
Tokenizer dirs(var, OS_PATH_VAR_SEP);
#if defined(_WIN32)
Tokenizer dirs(var, ";");
#else
Tokenizer dirs(var, ":");
#endif
while (dirs.has_next()) {
auto pot = utils::path_basename(filepath);
auto dir = dirs.next();
@ -1124,7 +1126,6 @@ std::string utils::get_potential_file_path(const std::string &path)
}
return "";
}
#undef OS_PATH_VAR_SEP
/* ----------------------------------------------------------------------
read first line of potential file

9
unittest/formats/test_file_operations.cpp
View File

@ -74,7 +74,8 @@ protected:
}
};
#define MAX_BUF_SIZE 128
static constexpr int MAX_BUF_SIZE=128;
TEST_F(FileOperationsTest, safe_fgets)
{
char buf[MAX_BUF_SIZE];
@ -110,7 +111,6 @@ TEST_F(FileOperationsTest, safe_fgets)
fclose(fp);
}
#define MAX_BUF_SIZE 128
TEST_F(FileOperationsTest, fgets_trunc)
{
char buf[MAX_BUF_SIZE];
@ -119,13 +119,15 @@ TEST_F(FileOperationsTest, fgets_trunc)
FILE *fp = fopen("safe_file_read_test.txt", "rb");
ASSERT_NE(fp, nullptr);
// read line shorter than buffer
memset(buf, 0, MAX_BUF_SIZE);
ptr = utils::fgets_trunc(buf, MAX_BUF_SIZE, fp);
ASSERT_THAT(buf, StrEq("one line\n"));
ASSERT_NE(ptr,nullptr);
// read line of exactly the buffer length
memset(buf, 0, MAX_BUF_SIZE);
ptr = utils::fgets_trunc(buf, MAX_BUF_SIZE, fp);
ptr = utils::fgets_trunc(buf, sizeof("two_lines\n"), fp);
ASSERT_THAT(buf, StrEq("two_lines\n"));
ASSERT_NE(ptr,nullptr);
@ -172,7 +174,6 @@ TEST_F(FileOperationsTest, fgets_trunc)
fclose(fp);
}
#define MAX_BUF_SIZE 128
TEST_F(FileOperationsTest, safe_fread)
{
char buf[MAX_BUF_SIZE];

70
unittest/utils/test_argutils.cpp
View File

@ -169,6 +169,16 @@ TEST(ArgInfo, variable2)
ASSERT_THAT(arg.get_name(), StrEq("x"));
}
TEST(ArgInfo, variable3)
{
ArgInfo arg("v_x[11][5]");
ASSERT_EQ(arg.get_dim(), 2);
ASSERT_EQ(arg.get_type(), ArgInfo::VARIABLE);
ASSERT_EQ(arg.get_index1(), 11);
ASSERT_EQ(arg.get_index2(), 5);
ASSERT_THAT(arg.get_name(), StrEq("x"));
}
TEST(ArgInfo, dname0)
{
ArgInfo arg("d_text", ArgInfo::DNAME);
@ -179,6 +189,36 @@ TEST(ArgInfo, dname0)
ASSERT_THAT(arg.get_name(), StrEq("text"));
}
TEST(ArgInfo, dname1)
{
ArgInfo arg("d2_text", ArgInfo::DNAME | ArgInfo::INAME);
ASSERT_EQ(arg.get_dim(), 0);
ASSERT_EQ(arg.get_type(), ArgInfo::DNAME);
ASSERT_EQ(arg.get_index1(), 0);
ASSERT_EQ(arg.get_index2(), -1);
ASSERT_THAT(arg.get_name(), StrEq("text"));
}
TEST(ArgInfo, dname2)
{
ArgInfo arg("d2_text[11]", ArgInfo::DNAME);
ASSERT_EQ(arg.get_dim(), 1);
ASSERT_EQ(arg.get_type(), ArgInfo::DNAME);
ASSERT_EQ(arg.get_index1(), 11);
ASSERT_EQ(arg.get_index2(), -1);
ASSERT_THAT(arg.get_name(), StrEq("text"));
}
TEST(ArgInfo, dname3)
{
ArgInfo arg("d2_text[24][11]", ArgInfo::DNAME);
ASSERT_EQ(arg.get_dim(), 2);
ASSERT_EQ(arg.get_type(), ArgInfo::DNAME);
ASSERT_EQ(arg.get_index1(), 24);
ASSERT_EQ(arg.get_index2(), 11);
ASSERT_THAT(arg.get_name(), StrEq("text"));
}
TEST(ArgInfo, iname0)
{
ArgInfo arg("i_text", ArgInfo::INAME);
@ -189,6 +229,36 @@ TEST(ArgInfo, iname0)
ASSERT_THAT(arg.get_name(), StrEq("text"));
}
TEST(ArgInfo, iname1)
{
ArgInfo arg("i2_text", ArgInfo::INAME);
ASSERT_EQ(arg.get_dim(), 0);
ASSERT_EQ(arg.get_type(), ArgInfo::INAME);
ASSERT_EQ(arg.get_index1(), 0);
ASSERT_EQ(arg.get_index2(), -1);
ASSERT_THAT(arg.get_name(), StrEq("text"));
}
TEST(ArgInfo, iname2)
{
ArgInfo arg("i2_text[2]", ArgInfo::INAME | ArgInfo::DNAME);
ASSERT_EQ(arg.get_dim(), 1);
ASSERT_EQ(arg.get_type(), ArgInfo::INAME);
ASSERT_EQ(arg.get_index1(), 2);
ASSERT_EQ(arg.get_index2(), -1);
ASSERT_THAT(arg.get_name(), StrEq("text"));
}
TEST(ArgInfo, iname3)
{
ArgInfo arg("i2_text[2][100]", ArgInfo::INAME | ArgInfo::DNAME);
ASSERT_EQ(arg.get_dim(), 2);
ASSERT_EQ(arg.get_type(), ArgInfo::INAME);
ASSERT_EQ(arg.get_index1(), 2);
ASSERT_EQ(arg.get_index2(), 100);
ASSERT_THAT(arg.get_name(), StrEq("text"));
}
TEST(ArgInfo, unsupported1)
{
ArgInfo arg("v_text[02][05]", ArgInfo::COMPUTE | ArgInfo::FIX);