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

enabled and apply clang-format

Browse Source
This commit is contained in:
Axel Kohlmeyer
2024-10-19 11:02:53 -04:00
parent 34ab2f862a
commit 760d871b7a
4 changed files with 405 additions and 360 deletions
Download Patch File Download Diff File Expand all files Collapse all files

156
src/compute_angle_local.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -13,38 +12,40 @@
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "compute_angle_local.h" #include "compute_angle_local.h"
#include <cmath>
#include <cstring> #include "angle.h"
#include "atom.h" #include "atom.h"
#include "atom_vec.h" #include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h" #include "domain.h"
#include "error.h"
#include "force.h" #include "force.h"
#include "angle.h"
#include "input.h" #include "input.h"
#include "variable.h"
#include "math_const.h" #include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "molecule.h"
#include "update.h"
#include "variable.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst; using namespace MathConst;
static constexpr int DELTA = 10000; static constexpr int DELTA = 10000;
enum{THETA,ENG,VARIABLE}; enum { THETA, ENG, VARIABLE };
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), tstr(nullptr), vstr(nullptr),
bstyle(nullptr), vvar(nullptr), tstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr) vlocal(nullptr), alocal(nullptr)
{ {
if (narg < 4) error->all(FLERR,"Illegal compute angle/local command"); if (narg < 4) error->all(FLERR, "Illegal compute angle/local command");
if (atom->avec->angles_allow == 0) if (atom->avec->angles_allow == 0)
error->all(FLERR,"Compute angle/local used when angles are not allowed"); error->all(FLERR, "Compute angle/local used when angles are not allowed");
local_flag = 1; local_flag = 1;
@ -52,7 +53,7 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
nvalues = narg - 3; nvalues = narg - 3;
bstyle = new int[nvalues]; bstyle = new int[nvalues];
vstr = new char*[nvalues]; vstr = new char *[nvalues];
vvar = new int[nvalues]; vvar = new int[nvalues];
nvalues = 0; nvalues = 0;
@ -61,16 +62,17 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
int iarg; int iarg;
for (iarg = 3; iarg < narg; iarg++) { for (iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"theta") == 0) { if (strcmp(arg[iarg], "theta") == 0) {
bstyle[nvalues++] = THETA; bstyle[nvalues++] = THETA;
tflag = 1; tflag = 1;
} else if (strcmp(arg[iarg],"eng") == 0) { } else if (strcmp(arg[iarg], "eng") == 0) {
bstyle[nvalues++] = ENG; bstyle[nvalues++] = ENG;
} else if (strncmp(arg[iarg],"v_",2) == 0) { } else if (strncmp(arg[iarg], "v_", 2) == 0) {
bstyle[nvalues++] = VARIABLE; bstyle[nvalues++] = VARIABLE;
vstr[nvar] = utils::strdup(&arg[iarg][2]); vstr[nvar] = utils::strdup(&arg[iarg][2]);
nvar++; nvar++;
} else break; } else
break;
} }
// optional args // optional args
@ -79,45 +81,46 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
tstr = nullptr; tstr = nullptr;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"set") == 0) { if (strcmp(arg[iarg], "set") == 0) {
setflag = 1; setflag = 1;
if (iarg+3 > narg) error->all(FLERR,"Illegal compute angle/local command"); if (iarg + 3 > narg) error->all(FLERR, "Illegal compute angle/local command");
if (strcmp(arg[iarg+1],"theta") == 0) { if (strcmp(arg[iarg + 1], "theta") == 0) {
delete [] tstr; delete[] tstr;
tstr = utils::strdup(arg[iarg+2]); tstr = utils::strdup(arg[iarg + 2]);
tflag = 1; tflag = 1;
} else error->all(FLERR,"Illegal compute angle/local command"); } else
error->all(FLERR, "Illegal compute angle/local command");
iarg += 3; iarg += 3;
} else error->all(FLERR,"Illegal compute angle/local command"); } else
error->all(FLERR, "Illegal compute angle/local command");
} }
// error check // error check
if (nvar) { if (nvar) {
if (!setflag) if (!setflag) error->all(FLERR, "Compute angle/local variable requires a set variable");
error->all(FLERR,"Compute angle/local variable requires a set variable");
for (int i = 0; i < nvar; i++) { for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]); vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0) if (vvar[i] < 0) error->all(FLERR, "Variable name for copute angle/local does not exist");
error->all(FLERR,"Variable name for copute angle/local does not exist");
if (!input->variable->equalstyle(vvar[i])) if (!input->variable->equalstyle(vvar[i]))
error->all(FLERR,"Variable for compute angle/local is invalid style"); error->all(FLERR, "Variable for compute angle/local is invalid style");
} }
if (tstr) { if (tstr) {
tvar = input->variable->find(tstr); tvar = input->variable->find(tstr);
if (tvar < 0) if (tvar < 0) error->all(FLERR, "Variable name for compute angle/local does not exist");
error->all(FLERR,"Variable name for compute angle/local does not exist");
if (!input->variable->internalstyle(tvar)) if (!input->variable->internalstyle(tvar))
error->all(FLERR,"Variable for compute angle/local is invalid style"); error->all(FLERR, "Variable for compute angle/local is invalid style");
} }
} else if (setflag) } else if (setflag)
error->all(FLERR,"Compute angle/local set with no variable"); error->all(FLERR, "Compute angle/local set with no variable");
// initialize output // initialize output
if (nvalues == 1) size_local_cols = 0; if (nvalues == 1)
else size_local_cols = nvalues; size_local_cols = 0;
else
size_local_cols = nvalues;
nmax = 0; nmax = 0;
vlocal = nullptr; vlocal = nullptr;
@ -128,12 +131,12 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeAngleLocal::~ComputeAngleLocal() ComputeAngleLocal::~ComputeAngleLocal()
{ {
delete [] bstyle; delete[] bstyle;
for (int i = 0; i < nvar; i++) delete [] vstr[i]; for (int i = 0; i < nvar; i++) delete[] vstr[i];
delete [] vstr; delete[] vstr;
delete [] vvar; delete[] vvar;
delete [] tstr; delete[] tstr;
memory->destroy(vlocal); memory->destroy(vlocal);
memory->destroy(alocal); memory->destroy(alocal);
@ -144,19 +147,17 @@ ComputeAngleLocal::~ComputeAngleLocal()
void ComputeAngleLocal::init() void ComputeAngleLocal::init()
{ {
if (force->angle == nullptr) if (force->angle == nullptr)
error->all(FLERR,"No angle style is defined for compute angle/local"); error->all(FLERR, "No angle style is defined for compute angle/local");
if (nvar) { if (nvar) {
for (int i = 0; i < nvar; i++) { for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]); vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0) if (vvar[i] < 0) error->all(FLERR, "Variable name for compute angle/local does not exist");
error->all(FLERR,"Variable name for compute angle/local does not exist");
} }
if (tstr) { if (tstr) {
tvar = input->variable->find(tstr); tvar = input->variable->find(tstr);
if (tvar < 0) if (tvar < 0) error->all(FLERR, "Variable name for compute angle/local does not exist");
error->all(FLERR,"Variable name for compute angle/local does not exist");
} }
} }
@ -194,10 +195,10 @@ void ComputeAngleLocal::compute_local()
int ComputeAngleLocal::compute_angles(int flag) int ComputeAngleLocal::compute_angles(int flag)
{ {
int i,m,na,atom1,atom2,atom3,imol,iatom,atype,ivar; int i, m, na, atom1, atom2, atom3, imol, iatom, atype, ivar;
tagint tagprev; tagint tagprev;
double delx1,dely1,delz1,delx2,dely2,delz2; double delx1, dely1, delz1, delx2, dely2, delz2;
double rsq1,rsq2,r1,r2,c,theta; double rsq1, rsq2, r1, r2, c, theta;
double *ptr; double *ptr;
double **x = atom->x; double **x = atom->x;
@ -224,7 +225,8 @@ int ComputeAngleLocal::compute_angles(int flag)
for (atom2 = 0; atom2 < nlocal; atom2++) { for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue; if (!(mask[atom2] & groupbit)) continue;
if (molecular == Atom::MOLECULAR) na = num_angle[atom2]; if (molecular == Atom::MOLECULAR)
na = num_angle[atom2];
else { else {
if (molindex[atom2] < 0) continue; if (molindex[atom2] < 0) continue;
imol = molindex[atom2]; imol = molindex[atom2];
@ -242,8 +244,8 @@ int ComputeAngleLocal::compute_angles(int flag)
if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue; if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue;
atype = onemols[imol]->angle_type[atom2][i]; atype = onemols[imol]->angle_type[atom2][i];
tagprev = tag[atom2] - iatom - 1; tagprev = tag[atom2] - iatom - 1;
atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i]+tagprev); atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i] + tagprev);
atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i]+tagprev); atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i] + tagprev);
} }
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
@ -261,50 +263,54 @@ int ComputeAngleLocal::compute_angles(int flag)
delx1 = x[atom1][0] - x[atom2][0]; delx1 = x[atom1][0] - x[atom2][0];
dely1 = x[atom1][1] - x[atom2][1]; dely1 = x[atom1][1] - x[atom2][1];
delz1 = x[atom1][2] - x[atom2][2]; delz1 = x[atom1][2] - x[atom2][2];
domain->minimum_image(delx1,dely1,delz1); domain->minimum_image(delx1, dely1, delz1);
rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1;
r1 = sqrt(rsq1); r1 = sqrt(rsq1);
delx2 = x[atom3][0] - x[atom2][0]; delx2 = x[atom3][0] - x[atom2][0];
dely2 = x[atom3][1] - x[atom2][1]; dely2 = x[atom3][1] - x[atom2][1];
delz2 = x[atom3][2] - x[atom2][2]; delz2 = x[atom3][2] - x[atom2][2];
domain->minimum_image(delx2,dely2,delz2); domain->minimum_image(delx2, dely2, delz2);
rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2;
r2 = sqrt(rsq2); r2 = sqrt(rsq2);
// c = cosine of angle // c = cosine of angle
// theta = angle in radians // theta = angle in radians
c = delx1*delx2 + dely1*dely2 + delz1*delz2; c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2;
c /= r1*r2; c /= r1 * r2;
if (c > 1.0) c = 1.0; if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
theta = acos(c); theta = acos(c);
} }
if (nvalues == 1) ptr = &vlocal[m]; if (nvalues == 1)
else ptr = alocal[m]; ptr = &vlocal[m];
else
ptr = alocal[m];
if (nvar) { if (nvar) {
ivar = 0; ivar = 0;
if (tstr) input->variable->internal_set(tvar,theta); if (tstr) input->variable->internal_set(tvar, theta);
} }
for (int n = 0; n < nvalues; n++) { for (int n = 0; n < nvalues; n++) {
switch (bstyle[n]) { switch (bstyle[n]) {
case THETA: case THETA:
ptr[n] = 180.0*theta/MY_PI; ptr[n] = 180.0 * theta / MY_PI;
break; break;
case ENG: case ENG:
if (atype > 0) ptr[n] = angle->single(atype,atom1,atom2,atom3); if (atype > 0)
else ptr[n] = 0.0; ptr[n] = angle->single(atype, atom1, atom2, atom3);
break; else
case VARIABLE: ptr[n] = 0.0;
ptr[n] = input->variable->compute_equal(vvar[ivar]); break;
ivar++; case VARIABLE:
break; ptr[n] = input->variable->compute_equal(vvar[ivar]);
ivar++;
break;
} }
} }
@ -325,11 +331,11 @@ void ComputeAngleLocal::reallocate(int n)
if (nvalues == 1) { if (nvalues == 1) {
memory->destroy(vlocal); memory->destroy(vlocal);
memory->create(vlocal,nmax,"angle/local:vector_local"); memory->create(vlocal, nmax, "angle/local:vector_local");
vector_local = vlocal; vector_local = vlocal;
} else { } else {
memory->destroy(alocal); memory->destroy(alocal);
memory->create(alocal,nmax,nvalues,"angle/local:array_local"); memory->create(alocal, nmax, nvalues, "angle/local:array_local");
array_local = alocal; array_local = alocal;
} }
} }
@ -340,6 +346,6 @@ void ComputeAngleLocal::reallocate(int n)
double ComputeAngleLocal::memory_usage() double ComputeAngleLocal::memory_usage()
{ {
double bytes = (double)nmax*nvalues * sizeof(double); double bytes = (double) nmax * nvalues * sizeof(double);
return bytes; return bytes;
} }

344
src/compute_bond_local.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -35,18 +34,35 @@ using namespace LAMMPS_NS;
static constexpr int DELTA = 10000; static constexpr int DELTA = 10000;
enum{DIST,DX,DY,DZ,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,FX,FY,FZ,VARIABLE,BN}; enum {
DIST,
DX,
DY,
DZ,
VELVIB,
OMEGA,
ENGTRANS,
ENGVIB,
ENGROT,
ENGPOT,
FORCE,
FX,
FY,
FZ,
VARIABLE,
BN
};
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), dstr(nullptr), vstr(nullptr),
bstyle(nullptr), vvar(nullptr), dstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr) vlocal(nullptr), alocal(nullptr)
{ {
if (narg < 4) error->all(FLERR,"Illegal compute bond/local command"); if (narg < 4) error->all(FLERR, "Illegal compute bond/local command");
if (atom->avec->bonds_allow == 0) if (atom->avec->bonds_allow == 0)
error->all(FLERR,"Compute bond/local used when bonds are not allowed"); error->all(FLERR, "Compute bond/local used when bonds are not allowed");
local_flag = 1; local_flag = 1;
comm_forward = 3; comm_forward = 3;
@ -56,7 +72,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
nvalues = narg - 3; nvalues = narg - 3;
bstyle = new int[nvalues]; bstyle = new int[nvalues];
bindex = new int[nvalues]; bindex = new int[nvalues];
vstr = new char*[nvalues]; vstr = new char *[nvalues];
vvar = new int[nvalues]; vvar = new int[nvalues];
nvalues = 0; nvalues = 0;
@ -64,30 +80,45 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
int iarg; int iarg;
for (iarg = 3; iarg < narg; iarg++) { for (iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST; if (strcmp(arg[iarg], "dist") == 0)
else if (strcmp(arg[iarg],"dx") == 0) bstyle[nvalues++] = DX; bstyle[nvalues++] = DIST;
else if (strcmp(arg[iarg],"dy") == 0) bstyle[nvalues++] = DY; else if (strcmp(arg[iarg], "dx") == 0)
else if (strcmp(arg[iarg],"dz") == 0) bstyle[nvalues++] = DZ; bstyle[nvalues++] = DX;
else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT; else if (strcmp(arg[iarg], "dy") == 0)
else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE; bstyle[nvalues++] = DY;
else if (strcmp(arg[iarg],"fx") == 0) bstyle[nvalues++] = FX; else if (strcmp(arg[iarg], "dz") == 0)
else if (strcmp(arg[iarg],"fy") == 0) bstyle[nvalues++] = FY; bstyle[nvalues++] = DZ;
else if (strcmp(arg[iarg],"fz") == 0) bstyle[nvalues++] = FZ; else if (strcmp(arg[iarg], "engpot") == 0)
else if (strcmp(arg[iarg],"engvib") == 0) bstyle[nvalues++] = ENGVIB; bstyle[nvalues++] = ENGPOT;
else if (strcmp(arg[iarg],"engrot") == 0) bstyle[nvalues++] = ENGROT; else if (strcmp(arg[iarg], "force") == 0)
else if (strcmp(arg[iarg],"engtrans") == 0) bstyle[nvalues++] = ENGTRANS; bstyle[nvalues++] = FORCE;
else if (strcmp(arg[iarg],"omega") == 0) bstyle[nvalues++] = OMEGA; else if (strcmp(arg[iarg], "fx") == 0)
else if (strcmp(arg[iarg],"velvib") == 0) bstyle[nvalues++] = VELVIB; bstyle[nvalues++] = FX;
else if (strncmp(arg[iarg],"v_",2) == 0) { else if (strcmp(arg[iarg], "fy") == 0)
bstyle[nvalues++] = FY;
else if (strcmp(arg[iarg], "fz") == 0)
bstyle[nvalues++] = FZ;
else if (strcmp(arg[iarg], "engvib") == 0)
bstyle[nvalues++] = ENGVIB;
else if (strcmp(arg[iarg], "engrot") == 0)
bstyle[nvalues++] = ENGROT;
else if (strcmp(arg[iarg], "engtrans") == 0)
bstyle[nvalues++] = ENGTRANS;
else if (strcmp(arg[iarg], "omega") == 0)
bstyle[nvalues++] = OMEGA;
else if (strcmp(arg[iarg], "velvib") == 0)
bstyle[nvalues++] = VELVIB;
else if (strncmp(arg[iarg], "v_", 2) == 0) {
bstyle[nvalues++] = VARIABLE; bstyle[nvalues++] = VARIABLE;
vstr[nvar] = utils::strdup(&arg[iarg][2]); vstr[nvar] = utils::strdup(&arg[iarg][2]);
nvar++; nvar++;
} else if (utils::strmatch(arg[iarg], "^b\\d+$")) { // b1, b2, b3, ... bN } else if (utils::strmatch(arg[iarg], "^b\\d+$")) { // b1, b2, b3, ... bN
int n = std::stoi(&arg[iarg][1]); int n = std::stoi(&arg[iarg][1]);
if (n <= 0) error->all(FLERR, "Invalid keyword {} in compute bond/local command", arg[iarg]); if (n <= 0) error->all(FLERR, "Invalid keyword {} in compute bond/local command", arg[iarg]);
bstyle[nvalues] = BN; bstyle[nvalues] = BN;
bindex[nvalues++] = n - 1; bindex[nvalues++] = n - 1;
} else break; } else
break;
} }
// optional args // optional args
@ -96,40 +127,39 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
dstr = nullptr; dstr = nullptr;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"set") == 0) { if (strcmp(arg[iarg], "set") == 0) {
setflag = 1; setflag = 1;
if (iarg+3 > narg) utils::missing_cmd_args(FLERR,"compute bond/local set", error); if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "compute bond/local set", error);
if (strcmp(arg[iarg+1],"dist") == 0) { if (strcmp(arg[iarg + 1], "dist") == 0) {
delete [] dstr; delete[] dstr;
dstr = utils::strdup(arg[iarg+2]); dstr = utils::strdup(arg[iarg + 2]);
} else error->all(FLERR,"Unknown compute bond/local set keyword: {}", arg[iarg+2]); } else
error->all(FLERR, "Unknown compute bond/local set keyword: {}", arg[iarg + 2]);
iarg += 3; iarg += 3;
} else error->all(FLERR,"Unknown compute bond/local keyword: {}", arg[iarg]); } else
error->all(FLERR, "Unknown compute bond/local keyword: {}", arg[iarg]);
} }
// error check // error check
if (nvar) { if (nvar) {
if (!setflag) if (!setflag) error->all(FLERR, "Compute bond/local variable requires a set variable");
error->all(FLERR,"Compute bond/local variable requires a set variable");
for (int i = 0; i < nvar; i++) { for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]); vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0) if (vvar[i] < 0)
error->all(FLERR,"Variable name {} for copute bond/local does not exist", vstr[i]); error->all(FLERR, "Variable name {} for copute bond/local does not exist", vstr[i]);
if (!input->variable->equalstyle(vvar[i])) if (!input->variable->equalstyle(vvar[i]))
error->all(FLERR,"Variable {} for compute bond/local is invalid style", vstr[i]); error->all(FLERR, "Variable {} for compute bond/local is invalid style", vstr[i]);
} }
if (dstr) { if (dstr) {
dvar = input->variable->find(dstr); dvar = input->variable->find(dstr);
if (dvar < 0) if (dvar < 0) error->all(FLERR, "Variable name for compute bond/local does not exist");
error->all(FLERR,"Variable name for compute bond/local does not exist");
if (!input->variable->internalstyle(dvar)) if (!input->variable->internalstyle(dvar))
error->all(FLERR,"Variable for compute bond/local is invalid style"); error->all(FLERR, "Variable for compute bond/local is invalid style");
} }
} else if (setflag) } else if (setflag)
error->all(FLERR,"Compute bond/local set used with without a variable"); error->all(FLERR, "Compute bond/local set used with without a variable");
// set singleflag if need to call bond->single() // set singleflag if need to call bond->single()
// set velflag if compute any quantities based on velocities // set velflag if compute any quantities based on velocities
@ -137,16 +167,20 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
singleflag = 0; singleflag = 0;
velflag = 0; velflag = 0;
for (int i = 0; i < nvalues; i++) { for (int i = 0; i < nvalues; i++) {
if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX || if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX || bstyle[i] == FY ||
bstyle[i] == FY || bstyle[i] == FZ || bstyle[i] == BN) singleflag = 1; bstyle[i] == FZ || bstyle[i] == BN)
if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS || singleflag = 1;
bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1; if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS || bstyle[i] == ENGVIB ||
bstyle[i] == ENGROT)
velflag = 1;
} }
// initialize output // initialize output
if (nvalues == 1) size_local_cols = 0; if (nvalues == 1)
else size_local_cols = nvalues; size_local_cols = 0;
else
size_local_cols = nvalues;
nmax = 0; nmax = 0;
vlocal = nullptr; vlocal = nullptr;
@ -157,13 +191,13 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeBondLocal::~ComputeBondLocal() ComputeBondLocal::~ComputeBondLocal()
{ {
delete [] bstyle; delete[] bstyle;
delete [] bindex; delete[] bindex;
for (int i = 0; i < nvar; i++) delete [] vstr[i]; for (int i = 0; i < nvar; i++) delete[] vstr[i];
delete [] vstr; delete[] vstr;
delete [] vvar; delete[] vvar;
delete [] dstr; delete[] dstr;
memory->destroy(vlocal); memory->destroy(vlocal);
memory->destroy(alocal); memory->destroy(alocal);
@ -173,8 +207,7 @@ ComputeBondLocal::~ComputeBondLocal()
void ComputeBondLocal::init() void ComputeBondLocal::init()
{ {
if (force->bond == nullptr) if (force->bond == nullptr) error->all(FLERR, "No bond style is defined for compute bond/local");
error->all(FLERR,"No bond style is defined for compute bond/local");
for (int i = 0; i < nvalues; i++) for (int i = 0; i < nvalues; i++)
if (bstyle[i] == BN && bindex[i] >= force->bond->single_extra) if (bstyle[i] == BN && bindex[i] >= force->bond->single_extra)
@ -183,21 +216,21 @@ void ComputeBondLocal::init()
if (nvar) { if (nvar) {
for (int i = 0; i < nvar; i++) { for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]); vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0) if (vvar[i] < 0) error->all(FLERR, "Variable name for compute bond/local does not exist");
error->all(FLERR,"Variable name for compute bond/local does not exist");
} }
if (dstr) { if (dstr) {
dvar = input->variable->find(dstr); dvar = input->variable->find(dstr);
if (dvar < 0) if (dvar < 0) error->all(FLERR, "Variable name for compute bond/local does not exist");
error->all(FLERR,"Variable name for compute bond/local does not exist");
} }
} }
// set ghostvelflag if need to acquire ghost atom velocities // set ghostvelflag if need to acquire ghost atom velocities
if (velflag && !comm->ghost_velocity) ghostvelflag = 1; if (velflag && !comm->ghost_velocity)
else ghostvelflag = 0; ghostvelflag = 1;
else
ghostvelflag = 0;
// do initial memory allocation so that memory_usage() is correct // do initial memory allocation so that memory_usage() is correct
@ -236,17 +269,17 @@ void ComputeBondLocal::compute_local()
int ComputeBondLocal::compute_bonds(int flag) int ComputeBondLocal::compute_bonds(int flag)
{ {
int i,m,nb,atom1,atom2,imol,iatom,btype,ivar; int i, m, nb, atom1, atom2, imol, iatom, btype, ivar;
tagint tagprev; tagint tagprev;
double dx,dy,dz,rsq; double dx, dy, dz, rsq;
double mass1,mass2,masstotal,invmasstotal; double mass1, mass2, masstotal, invmasstotal;
double xcm[3],vcm[3]; double xcm[3], vcm[3];
double delr1[3],delr2[3],delv1[3],delv2[3]; double delr1[3], delr2[3], delv1[3], delv2[3];
double r12[3],vpar1,vpar2; double r12[3], vpar1, vpar2;
double vvib,vrotsq; double vvib, vrotsq;
double inertia,omegasq; double inertia, omegasq;
double mvv2e; double mvv2e;
double engpot,engtrans,engvib,engrot,fbond; double engpot, engtrans, engvib, engrot, fbond;
double *ptr; double *ptr;
double **x = atom->x; double **x = atom->x;
@ -280,7 +313,8 @@ int ComputeBondLocal::compute_bonds(int flag)
for (atom1 = 0; atom1 < nlocal; atom1++) { for (atom1 = 0; atom1 < nlocal; atom1++) {
if (!(mask[atom1] & groupbit)) continue; if (!(mask[atom1] & groupbit)) continue;
if (molecular == Atom::MOLECULAR) nb = num_bond[atom1]; if (molecular == Atom::MOLECULAR)
nb = num_bond[atom1];
else { else {
if (molindex[atom1] < 0) continue; if (molindex[atom1] < 0) continue;
imol = molindex[atom1]; imol = molindex[atom1];
@ -295,7 +329,7 @@ int ComputeBondLocal::compute_bonds(int flag)
} else { } else {
tagprev = tag[atom1] - iatom - 1; tagprev = tag[atom1] - iatom - 1;
btype = onemols[imol]->bond_type[iatom][i]; btype = onemols[imol]->bond_type[iatom][i];
atom2 = atom->map(onemols[imol]->bond_atom[iatom][i]+tagprev); atom2 = atom->map(onemols[imol]->bond_atom[iatom][i] + tagprev);
} }
if (atom2 < 0 || !(mask[atom2] & groupbit)) continue; if (atom2 < 0 || !(mask[atom2] & groupbit)) continue;
@ -310,55 +344,56 @@ int ComputeBondLocal::compute_bonds(int flag)
dx = x[atom1][0] - x[atom2][0]; dx = x[atom1][0] - x[atom2][0];
dy = x[atom1][1] - x[atom2][1]; dy = x[atom1][1] - x[atom2][1];
dz = x[atom1][2] - x[atom2][2]; dz = x[atom1][2] - x[atom2][2];
domain->minimum_image(dx,dy,dz); domain->minimum_image(dx, dy, dz);
rsq = dx*dx + dy*dy + dz*dz; rsq = dx * dx + dy * dy + dz * dz;
if (btype == 0) { if (btype == 0) {
fbond = 0.0; fbond = 0.0;
} else { } else {
if (singleflag) engpot = bond->single(btype,rsq,atom1,atom2,fbond); if (singleflag)
else fbond = engpot = 0.0; engpot = bond->single(btype, rsq, atom1, atom2, fbond);
else
fbond = engpot = 0.0;
engvib = engrot = engtrans = omegasq = vvib = 0.0; engvib = engrot = engtrans = omegasq = vvib = 0.0;
if (velflag) { if (velflag) {
if (rmass) { if (rmass) {
mass1 = rmass[atom1]; mass1 = rmass[atom1];
mass2 = rmass[atom2]; mass2 = rmass[atom2];
} } else {
else {
mass1 = mass[type[atom1]]; mass1 = mass[type[atom1]];
mass2 = mass[type[atom2]]; mass2 = mass[type[atom2]];
} }
masstotal = mass1+mass2; masstotal = mass1 + mass2;
invmasstotal = 1.0 / (masstotal); invmasstotal = 1.0 / (masstotal);
xcm[0] = (mass1*x[atom1][0] + mass2*x[atom2][0]) * invmasstotal; xcm[0] = (mass1 * x[atom1][0] + mass2 * x[atom2][0]) * invmasstotal;
xcm[1] = (mass1*x[atom1][1] + mass2*x[atom2][1]) * invmasstotal; xcm[1] = (mass1 * x[atom1][1] + mass2 * x[atom2][1]) * invmasstotal;
xcm[2] = (mass1*x[atom1][2] + mass2*x[atom2][2]) * invmasstotal; xcm[2] = (mass1 * x[atom1][2] + mass2 * x[atom2][2]) * invmasstotal;
vcm[0] = (mass1*v[atom1][0] + mass2*v[atom2][0]) * invmasstotal; vcm[0] = (mass1 * v[atom1][0] + mass2 * v[atom2][0]) * invmasstotal;
vcm[1] = (mass1*v[atom1][1] + mass2*v[atom2][1]) * invmasstotal; vcm[1] = (mass1 * v[atom1][1] + mass2 * v[atom2][1]) * invmasstotal;
vcm[2] = (mass1*v[atom1][2] + mass2*v[atom2][2]) * invmasstotal; vcm[2] = (mass1 * v[atom1][2] + mass2 * v[atom2][2]) * invmasstotal;
engtrans = 0.5 * masstotal * MathExtra::lensq3(vcm); engtrans = 0.5 * masstotal * MathExtra::lensq3(vcm);
// r12 = unit bond vector from atom1 to atom2 // r12 = unit bond vector from atom1 to atom2
MathExtra::sub3(x[atom2],x[atom1],r12); MathExtra::sub3(x[atom2], x[atom1], r12);
MathExtra::norm3(r12); MathExtra::norm3(r12);
// delr = vector from COM to each atom // delr = vector from COM to each atom
// delv = velocity of each atom relative to COM // delv = velocity of each atom relative to COM
MathExtra::sub3(x[atom1],xcm,delr1); MathExtra::sub3(x[atom1], xcm, delr1);
MathExtra::sub3(x[atom2],xcm,delr2); MathExtra::sub3(x[atom2], xcm, delr2);
MathExtra::sub3(v[atom1],vcm,delv1); MathExtra::sub3(v[atom1], vcm, delv1);
MathExtra::sub3(v[atom2],vcm,delv2); MathExtra::sub3(v[atom2], vcm, delv2);
// vpar = component of delv parallel to bond vector // vpar = component of delv parallel to bond vector
vpar1 = MathExtra::dot3(delv1,r12); vpar1 = MathExtra::dot3(delv1, r12);
vpar2 = MathExtra::dot3(delv2,r12); vpar2 = MathExtra::dot3(delv2, r12);
engvib = 0.5 * (mass1*vpar1*vpar1 + mass2*vpar2*vpar2); engvib = 0.5 * (mass1 * vpar1 * vpar1 + mass2 * vpar2 * vpar2);
// vvib = relative velocity of 2 atoms along bond direction // vvib = relative velocity of 2 atoms along bond direction
// vvib < 0 for 2 atoms moving towards each other // vvib < 0 for 2 atoms moving towards each other
@ -369,9 +404,8 @@ int ComputeBondLocal::compute_bonds(int flag)
// vrotsq = tangential speed squared of atom1 only // vrotsq = tangential speed squared of atom1 only
// omegasq = omega squared, and is the same for atom1 and atom2 // omegasq = omega squared, and is the same for atom1 and atom2
inertia = mass1*MathExtra::lensq3(delr1) + inertia = mass1 * MathExtra::lensq3(delr1) + mass2 * MathExtra::lensq3(delr2);
mass2*MathExtra::lensq3(delr2); vrotsq = MathExtra::lensq3(delv1) - vpar1 * vpar1;
vrotsq = MathExtra::lensq3(delv1) - vpar1*vpar1;
omegasq = vrotsq / MathExtra::lensq3(delr1); omegasq = vrotsq / MathExtra::lensq3(delr1);
engrot = 0.5 * inertia * omegasq; engrot = 0.5 * inertia * omegasq;
@ -384,12 +418,14 @@ int ComputeBondLocal::compute_bonds(int flag)
engrot *= mvv2e; engrot *= mvv2e;
} }
if (nvalues == 1) ptr = &vlocal[m]; if (nvalues == 1)
else ptr = alocal[m]; ptr = &vlocal[m];
else
ptr = alocal[m];
if (nvar) { if (nvar) {
ivar = 0; ivar = 0;
if (dstr) input->variable->internal_set(dvar,sqrt(rsq)); if (dstr) input->variable->internal_set(dvar, sqrt(rsq));
} }
// to make sure dx, dy and dz are always from the lower to the higher id // to make sure dx, dy and dz are always from the lower to the higher id
@ -397,55 +433,55 @@ int ComputeBondLocal::compute_bonds(int flag)
for (int n = 0; n < nvalues; n++) { for (int n = 0; n < nvalues; n++) {
switch (bstyle[n]) { switch (bstyle[n]) {
case DIST: case DIST:
ptr[n] = sqrt(rsq); ptr[n] = sqrt(rsq);
break; break;
case DX: case DX:
ptr[n] = dx*directionCorrection; ptr[n] = dx * directionCorrection;
break; break;
case DY: case DY:
ptr[n] = dy*directionCorrection; ptr[n] = dy * directionCorrection;
break; break;
case DZ: case DZ:
ptr[n] = dz*directionCorrection; ptr[n] = dz * directionCorrection;
break; break;
case ENGPOT: case ENGPOT:
ptr[n] = engpot; ptr[n] = engpot;
break; break;
case FORCE: case FORCE:
ptr[n] = sqrt(rsq)*fbond; ptr[n] = sqrt(rsq) * fbond;
break; break;
case FX: case FX:
ptr[n] = dx*fbond; ptr[n] = dx * fbond;
break; break;
case FY: case FY:
ptr[n] = dy*fbond; ptr[n] = dy * fbond;
break; break;
case FZ: case FZ:
ptr[n] = dz*fbond; ptr[n] = dz * fbond;
break; break;
case ENGVIB: case ENGVIB:
ptr[n] = engvib; ptr[n] = engvib;
break; break;
case ENGROT: case ENGROT:
ptr[n] = engrot; ptr[n] = engrot;
break; break;
case ENGTRANS: case ENGTRANS:
ptr[n] = engtrans; ptr[n] = engtrans;
break; break;
case OMEGA: case OMEGA:
ptr[n] = sqrt(omegasq); ptr[n] = sqrt(omegasq);
break; break;
case VELVIB: case VELVIB:
ptr[n] = vvib; ptr[n] = vvib;
break; break;
case VARIABLE: case VARIABLE:
ptr[n] = input->variable->compute_equal(vvar[ivar]); ptr[n] = input->variable->compute_equal(vvar[ivar]);
ivar++; ivar++;
break; break;
case BN: case BN:
ptr[n] = bond->svector[bindex[n]]; ptr[n] = bond->svector[bindex[n]];
break; break;
} }
} }
} }
@ -459,10 +495,10 @@ int ComputeBondLocal::compute_bonds(int flag)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf, int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/,
int /*pbc_flag*/, int * /*pbc*/) int * /*pbc*/)
{ {
int i,j,m; int i, j, m;
double **v = atom->v; double **v = atom->v;
@ -481,7 +517,7 @@ int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf,
void ComputeBondLocal::unpack_forward_comm(int n, int first, double *buf) void ComputeBondLocal::unpack_forward_comm(int n, int first, double *buf)
{ {
int i,m,last; int i, m, last;
double **v = atom->v; double **v = atom->v;
@ -504,11 +540,11 @@ void ComputeBondLocal::reallocate(int n)
if (nvalues == 1) { if (nvalues == 1) {
memory->destroy(vlocal); memory->destroy(vlocal);
memory->create(vlocal,nmax,"bond/local:vector_local"); memory->create(vlocal, nmax, "bond/local:vector_local");
vector_local = vlocal; vector_local = vlocal;
} else { } else {
memory->destroy(alocal); memory->destroy(alocal);
memory->create(alocal,nmax,nvalues,"bond/local:array_local"); memory->create(alocal, nmax, nvalues, "bond/local:array_local");
array_local = alocal; array_local = alocal;
} }
} }
@ -519,6 +555,6 @@ void ComputeBondLocal::reallocate(int n)
double ComputeBondLocal::memory_usage() double ComputeBondLocal::memory_usage()
{ {
double bytes = (double)nmax*nvalues * sizeof(double); double bytes = (double) nmax * nvalues * sizeof(double);
return bytes; return bytes;
} }

176
src/compute_dihedral_local.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -13,19 +12,21 @@
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "compute_dihedral_local.h" #include "compute_dihedral_local.h"
#include <cmath>
#include <cstring>
#include "atom.h" #include "atom.h"
#include "atom_vec.h" #include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h" #include "domain.h"
#include "error.h"
#include "force.h" #include "force.h"
#include "input.h" #include "input.h"
#include "variable.h"
#include "math_const.h" #include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "molecule.h"
#include "update.h"
#include "variable.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst; using namespace MathConst;
@ -37,14 +38,13 @@ enum { PHI, VARIABLE };
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), pstr(nullptr), vstr(nullptr),
bstyle(nullptr), vvar(nullptr), pstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr) vlocal(nullptr), alocal(nullptr)
{ {
if (narg < 4) error->all(FLERR,"Illegal compute dihedral/local command"); if (narg < 4) error->all(FLERR, "Illegal compute dihedral/local command");
if (atom->avec->dihedrals_allow == 0) if (atom->avec->dihedrals_allow == 0)
error->all(FLERR, error->all(FLERR, "Compute dihedral/local used when dihedrals are not allowed");
"Compute dihedral/local used when dihedrals are not allowed");
local_flag = 1; local_flag = 1;
@ -52,7 +52,7 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
nvalues = narg - 3; nvalues = narg - 3;
bstyle = new int[nvalues]; bstyle = new int[nvalues];
vstr = new char*[nvalues]; vstr = new char *[nvalues];
vvar = new int[nvalues]; vvar = new int[nvalues];
nvalues = 0; nvalues = 0;
@ -60,13 +60,14 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
int iarg; int iarg;
for (iarg = 3; iarg < narg; iarg++) { for (iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"phi") == 0) { if (strcmp(arg[iarg], "phi") == 0) {
bstyle[nvalues++] = PHI; bstyle[nvalues++] = PHI;
} else if (strncmp(arg[iarg],"v_",2) == 0) { } else if (strncmp(arg[iarg], "v_", 2) == 0) {
bstyle[nvalues++] = VARIABLE; bstyle[nvalues++] = VARIABLE;
vstr[nvar] = utils::strdup(&arg[iarg][2]); vstr[nvar] = utils::strdup(&arg[iarg][2]);
nvar++; nvar++;
} else break; } else
break;
} }
// optional args // optional args
@ -75,47 +76,45 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
pstr = nullptr; pstr = nullptr;
while (iarg < narg) { while (iarg < narg) {
if (strcmp(arg[iarg],"set") == 0) { if (strcmp(arg[iarg], "set") == 0) {
setflag = 1; setflag = 1;
if (iarg+3 > narg) if (iarg + 3 > narg) error->all(FLERR, "Illegal compute dihedral/local command");
error->all(FLERR,"Illegal compute dihedral/local command"); if (strcmp(arg[iarg + 1], "phi") == 0) {
if (strcmp(arg[iarg+1],"phi") == 0) { delete[] pstr;
delete [] pstr; pstr = utils::strdup(arg[iarg + 2]);
pstr = utils::strdup(arg[iarg+2]); } else
} else error->all(FLERR,"Illegal compute dihedral/local command"); error->all(FLERR, "Illegal compute dihedral/local command");
iarg += 3; iarg += 3;
} else error->all(FLERR,"Illegal compute dihedral/local command"); } else
error->all(FLERR, "Illegal compute dihedral/local command");
} }
// error check // error check
if (nvar) { if (nvar) {
if (!setflag) if (!setflag) error->all(FLERR, "Compute dihedral/local variable requires a set variable");
error->all(FLERR,"Compute dihedral/local variable requires a set variable");
for (int i = 0; i < nvar; i++) { for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]); vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0) if (vvar[i] < 0) error->all(FLERR, "Variable name for copute dihedral/local does not exist");
error->all(FLERR,
"Variable name for copute dihedral/local does not exist");
if (!input->variable->equalstyle(vvar[i])) if (!input->variable->equalstyle(vvar[i]))
error->all(FLERR,"Variable for compute dihedral/local is invalid style"); error->all(FLERR, "Variable for compute dihedral/local is invalid style");
} }
if (pstr) { if (pstr) {
pvar = input->variable->find(pstr); pvar = input->variable->find(pstr);
if (pvar < 0) if (pvar < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist");
error->all(FLERR,
"Variable name for compute dihedral/local does not exist");
if (!input->variable->internalstyle(pvar)) if (!input->variable->internalstyle(pvar))
error->all(FLERR,"Variable for compute dihedral/local is invalid style"); error->all(FLERR, "Variable for compute dihedral/local is invalid style");
} }
} else if (setflag) } else if (setflag)
error->all(FLERR,"Compute dihedral/local set with no variable"); error->all(FLERR, "Compute dihedral/local set with no variable");
// initialize output // initialize output
if (nvalues == 1) size_local_cols = 0; if (nvalues == 1)
else size_local_cols = nvalues; size_local_cols = 0;
else
size_local_cols = nvalues;
nmax = 0; nmax = 0;
vlocal = nullptr; vlocal = nullptr;
@ -126,12 +125,12 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) :
ComputeDihedralLocal::~ComputeDihedralLocal() ComputeDihedralLocal::~ComputeDihedralLocal()
{ {
delete [] bstyle; delete[] bstyle;
for (int i = 0; i < nvar; i++) delete [] vstr[i]; for (int i = 0; i < nvar; i++) delete[] vstr[i];
delete [] vstr; delete[] vstr;
delete [] vvar; delete[] vvar;
delete [] pstr; delete[] pstr;
memory->destroy(vlocal); memory->destroy(vlocal);
memory->destroy(alocal); memory->destroy(alocal);
@ -142,21 +141,17 @@ ComputeDihedralLocal::~ComputeDihedralLocal()
void ComputeDihedralLocal::init() void ComputeDihedralLocal::init()
{ {
if (force->dihedral == nullptr) if (force->dihedral == nullptr)
error->all(FLERR,"No dihedral style is defined for compute dihedral/local"); error->all(FLERR, "No dihedral style is defined for compute dihedral/local");
if (nvar) { if (nvar) {
for (int i = 0; i < nvar; i++) { for (int i = 0; i < nvar; i++) {
vvar[i] = input->variable->find(vstr[i]); vvar[i] = input->variable->find(vstr[i]);
if (vvar[i] < 0) if (vvar[i] < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist");
error->all(FLERR,
"Variable name for compute dihedral/local does not exist");
} }
if (pstr) { if (pstr) {
pvar = input->variable->find(pstr); pvar = input->variable->find(pstr);
if (pvar < 0) if (pvar < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist");
error->all(FLERR,
"Variable name for compute dihedral/local does not exist");
} }
} }
@ -191,11 +186,11 @@ void ComputeDihedralLocal::compute_local()
int ComputeDihedralLocal::compute_dihedrals(int flag) int ComputeDihedralLocal::compute_dihedrals(int flag)
{ {
int i,m,nd,atom1,atom2,atom3,atom4,imol,iatom,ivar; int i, m, nd, atom1, atom2, atom3, atom4, imol, iatom, ivar;
tagint tagprev; tagint tagprev;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm;
double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,ra2inv,rb2inv,rabinv; double ax, ay, az, bx, by, bz, rasq, rbsq, rgsq, rg, ra2inv, rb2inv, rabinv;
double s,c,phi; double s, c, phi;
double *ptr; double *ptr;
double **x = atom->x; double **x = atom->x;
@ -220,7 +215,8 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
for (atom2 = 0; atom2 < nlocal; atom2++) { for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue; if (!(mask[atom2] & groupbit)) continue;
if (molecular == Atom::MOLECULAR) nd = num_dihedral[atom2]; if (molecular == Atom::MOLECULAR)
nd = num_dihedral[atom2];
else { else {
if (molindex[atom2] < 0) continue; if (molindex[atom2] < 0) continue;
imol = molindex[atom2]; imol = molindex[atom2];
@ -237,9 +233,9 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
} else { } else {
if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue;
tagprev = tag[atom2] - iatom - 1; tagprev = tag[atom2] - iatom - 1;
atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i]+tagprev); atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i] + tagprev);
atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i]+tagprev); atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i] + tagprev);
atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i]+tagprev); atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i] + tagprev);
} }
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
@ -256,64 +252,66 @@ int ComputeDihedralLocal::compute_dihedrals(int flag)
vb1x = x[atom1][0] - x[atom2][0]; vb1x = x[atom1][0] - x[atom2][0];
vb1y = x[atom1][1] - x[atom2][1]; vb1y = x[atom1][1] - x[atom2][1];
vb1z = x[atom1][2] - x[atom2][2]; vb1z = x[atom1][2] - x[atom2][2];
domain->minimum_image(vb1x,vb1y,vb1z); domain->minimum_image(vb1x, vb1y, vb1z);
vb2x = x[atom3][0] - x[atom2][0]; vb2x = x[atom3][0] - x[atom2][0];
vb2y = x[atom3][1] - x[atom2][1]; vb2y = x[atom3][1] - x[atom2][1];
vb2z = x[atom3][2] - x[atom2][2]; vb2z = x[atom3][2] - x[atom2][2];
domain->minimum_image(vb2x,vb2y,vb2z); domain->minimum_image(vb2x, vb2y, vb2z);
vb2xm = -vb2x; vb2xm = -vb2x;
vb2ym = -vb2y; vb2ym = -vb2y;
vb2zm = -vb2z; vb2zm = -vb2z;
domain->minimum_image(vb2xm,vb2ym,vb2zm); domain->minimum_image(vb2xm, vb2ym, vb2zm);
vb3x = x[atom4][0] - x[atom3][0]; vb3x = x[atom4][0] - x[atom3][0];
vb3y = x[atom4][1] - x[atom3][1]; vb3y = x[atom4][1] - x[atom3][1];
vb3z = x[atom4][2] - x[atom3][2]; vb3z = x[atom4][2] - x[atom3][2];
domain->minimum_image(vb3x,vb3y,vb3z); domain->minimum_image(vb3x, vb3y, vb3z);
ax = vb1y*vb2zm - vb1z*vb2ym; ax = vb1y * vb2zm - vb1z * vb2ym;
ay = vb1z*vb2xm - vb1x*vb2zm; ay = vb1z * vb2xm - vb1x * vb2zm;
az = vb1x*vb2ym - vb1y*vb2xm; az = vb1x * vb2ym - vb1y * vb2xm;
bx = vb3y*vb2zm - vb3z*vb2ym; bx = vb3y * vb2zm - vb3z * vb2ym;
by = vb3z*vb2xm - vb3x*vb2zm; by = vb3z * vb2xm - vb3x * vb2zm;
bz = vb3x*vb2ym - vb3y*vb2xm; bz = vb3x * vb2ym - vb3y * vb2xm;
rasq = ax*ax + ay*ay + az*az; rasq = ax * ax + ay * ay + az * az;
rbsq = bx*bx + by*by + bz*bz; rbsq = bx * bx + by * by + bz * bz;
rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm;
rg = sqrt(rgsq); rg = sqrt(rgsq);
ra2inv = rb2inv = 0.0; ra2inv = rb2inv = 0.0;
if (rasq > 0) ra2inv = 1.0/rasq; if (rasq > 0) ra2inv = 1.0 / rasq;
if (rbsq > 0) rb2inv = 1.0/rbsq; if (rbsq > 0) rb2inv = 1.0 / rbsq;
rabinv = sqrt(ra2inv*rb2inv); rabinv = sqrt(ra2inv * rb2inv);
c = (ax*bx + ay*by + az*bz)*rabinv; c = (ax * bx + ay * by + az * bz) * rabinv;
s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); s = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z);
if (c > 1.0) c = 1.0; if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
phi = atan2(s,c); phi = atan2(s, c);
if (nvalues == 1) ptr = &vlocal[m]; if (nvalues == 1)
else ptr = alocal[m]; ptr = &vlocal[m];
else
ptr = alocal[m];
if (nvar) { if (nvar) {
ivar = 0; ivar = 0;
if (pstr) input->variable->internal_set(pvar,phi); if (pstr) input->variable->internal_set(pvar, phi);
} }
for (int n = 0; n < nvalues; n++) { for (int n = 0; n < nvalues; n++) {
switch (bstyle[n]) { switch (bstyle[n]) {
case PHI: case PHI:
ptr[n] = 180.0*phi/MY_PI; ptr[n] = 180.0 * phi / MY_PI;
break; break;
case VARIABLE: case VARIABLE:
ptr[n] = input->variable->compute_equal(vvar[ivar]); ptr[n] = input->variable->compute_equal(vvar[ivar]);
ivar++; ivar++;
break; break;
} }
} }
@ -334,11 +332,11 @@ void ComputeDihedralLocal::reallocate(int n)
if (nvalues == 1) { if (nvalues == 1) {
memory->destroy(vlocal); memory->destroy(vlocal);
memory->create(vlocal,nmax,"dihedral/local:vector_local"); memory->create(vlocal, nmax, "dihedral/local:vector_local");
vector_local = vlocal; vector_local = vlocal;
} else { } else {
memory->destroy(alocal); memory->destroy(alocal);
memory->create(alocal,nmax,nvalues,"dihedral/local:array_local"); memory->create(alocal, nmax, nvalues, "dihedral/local:array_local");
array_local = alocal; array_local = alocal;
} }
} }
@ -349,6 +347,6 @@ void ComputeDihedralLocal::reallocate(int n)
double ComputeDihedralLocal::memory_usage() double ComputeDihedralLocal::memory_usage()
{ {
double bytes = (double)nmax*nvalues * sizeof(double); double bytes = (double) nmax * nvalues * sizeof(double);
return bytes; return bytes;
} }

89
src/compute_improper_local.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ---------------------------------------------------------------------- /* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories https://www.lammps.org/, Sandia National Laboratories
@ -13,36 +12,35 @@
------------------------------------------------------------------------- */ ------------------------------------------------------------------------- */
#include "compute_improper_local.h" #include "compute_improper_local.h"
#include <cmath>
#include <cstring>
#include "atom.h" #include "atom.h"
#include "atom_vec.h" #include "atom_vec.h"
#include "molecule.h"
#include "update.h"
#include "domain.h" #include "domain.h"
#include "error.h"
#include "force.h" #include "force.h"
#include "math_const.h" #include "math_const.h"
#include "memory.h" #include "memory.h"
#include "error.h" #include "molecule.h"
#include "update.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS; using namespace LAMMPS_NS;
using namespace MathConst; using namespace MathConst;
static constexpr int DELTA = 10000; static constexpr int DELTA = 10000;
static constexpr double SMALL = 0.001;
static constexpr double SMALL = 0.001;
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) : ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), Compute(lmp, narg, arg), vlocal(nullptr), alocal(nullptr)
vlocal(nullptr), alocal(nullptr)
{ {
if (narg < 4) error->all(FLERR,"Illegal compute improper/local command"); if (narg < 4) error->all(FLERR, "Illegal compute improper/local command");
if (atom->avec->impropers_allow == 0) if (atom->avec->impropers_allow == 0)
error->all(FLERR, error->all(FLERR, "Compute improper/local used when impropers are not allowed");
"Compute improper/local used when impropers are not allowed");
local_flag = 1; local_flag = 1;
nvalues = narg - 3; nvalues = narg - 3;
@ -50,12 +48,16 @@ ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) :
nvalues = 0; nvalues = 0;
for (int iarg = 3; iarg < narg; iarg++) { for (int iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"chi") == 0) cflag = nvalues++; if (strcmp(arg[iarg], "chi") == 0)
else error->all(FLERR,"Invalid keyword in compute improper/local command"); cflag = nvalues++;
else
error->all(FLERR, "Invalid keyword in compute improper/local command");
} }
if (nvalues == 1) size_local_cols = 0; if (nvalues == 1)
else size_local_cols = nvalues; size_local_cols = 0;
else
size_local_cols = nvalues;
nmax = 0; nmax = 0;
vlocal = nullptr; vlocal = nullptr;
@ -75,7 +77,7 @@ ComputeImproperLocal::~ComputeImproperLocal()
void ComputeImproperLocal::init() void ComputeImproperLocal::init()
{ {
if (force->improper == nullptr) if (force->improper == nullptr)
error->all(FLERR,"No improper style is defined for compute improper/local"); error->all(FLERR, "No improper style is defined for compute improper/local");
// do initial memory allocation so that memory_usage() is correct // do initial memory allocation so that memory_usage() is correct
@ -108,11 +110,11 @@ void ComputeImproperLocal::compute_local()
int ComputeImproperLocal::compute_impropers(int flag) int ComputeImproperLocal::compute_impropers(int flag)
{ {
int i,m,n,ni,atom1,atom2,atom3,atom4,imol,iatom; int i, m, n, ni, atom1, atom2, atom3, atom4, imol, iatom;
tagint tagprev; tagint tagprev;
double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z; double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z;
double ss1,ss2,ss3,r1,r2,r3,c0,c1,c2,s1,s2; double ss1, ss2, ss3, r1, r2, r3, c0, c1, c2, s1, s2;
double s12,c; double s12, c;
double *cbuf; double *cbuf;
double **x = atom->x; double **x = atom->x;
@ -135,8 +137,10 @@ int ComputeImproperLocal::compute_impropers(int flag)
if (nvalues == 1) { if (nvalues == 1) {
if (cflag >= 0) cbuf = vlocal; if (cflag >= 0) cbuf = vlocal;
} else { } else {
if (cflag >= 0 && alocal) cbuf = &alocal[0][cflag]; if (cflag >= 0 && alocal)
else cbuf = nullptr; cbuf = &alocal[0][cflag];
else
cbuf = nullptr;
} }
} }
@ -144,7 +148,8 @@ int ComputeImproperLocal::compute_impropers(int flag)
for (atom2 = 0; atom2 < nlocal; atom2++) { for (atom2 = 0; atom2 < nlocal; atom2++) {
if (!(mask[atom2] & groupbit)) continue; if (!(mask[atom2] & groupbit)) continue;
if (molecular == Atom::MOLECULAR) ni = num_improper[atom2]; if (molecular == Atom::MOLECULAR)
ni = num_improper[atom2];
else { else {
if (molindex[atom2] < 0) continue; if (molindex[atom2] < 0) continue;
imol = molindex[atom2]; imol = molindex[atom2];
@ -161,9 +166,9 @@ int ComputeImproperLocal::compute_impropers(int flag)
} else { } else {
if (tag[atom2] != onemols[imol]->improper_atom2[atom2][i]) continue; if (tag[atom2] != onemols[imol]->improper_atom2[atom2][i]) continue;
tagprev = tag[atom2] - iatom - 1; tagprev = tag[atom2] - iatom - 1;
atom1 = atom->map(onemols[imol]->improper_atom1[atom2][i]+tagprev); atom1 = atom->map(onemols[imol]->improper_atom1[atom2][i] + tagprev);
atom3 = atom->map(onemols[imol]->improper_atom3[atom2][i]+tagprev); atom3 = atom->map(onemols[imol]->improper_atom3[atom2][i] + tagprev);
atom4 = atom->map(onemols[imol]->improper_atom4[atom2][i]+tagprev); atom4 = atom->map(onemols[imol]->improper_atom4[atom2][i] + tagprev);
} }
if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; if (atom1 < 0 || !(mask[atom1] & groupbit)) continue;
@ -178,21 +183,21 @@ int ComputeImproperLocal::compute_impropers(int flag)
vb1x = x[atom1][0] - x[atom2][0]; vb1x = x[atom1][0] - x[atom2][0];
vb1y = x[atom1][1] - x[atom2][1]; vb1y = x[atom1][1] - x[atom2][1];
vb1z = x[atom1][2] - x[atom2][2]; vb1z = x[atom1][2] - x[atom2][2];
domain->minimum_image(vb1x,vb1y,vb1z); domain->minimum_image(vb1x, vb1y, vb1z);
vb2x = x[atom3][0] - x[atom2][0]; vb2x = x[atom3][0] - x[atom2][0];
vb2y = x[atom3][1] - x[atom2][1]; vb2y = x[atom3][1] - x[atom2][1];
vb2z = x[atom3][2] - x[atom2][2]; vb2z = x[atom3][2] - x[atom2][2];
domain->minimum_image(vb2x,vb2y,vb2z); domain->minimum_image(vb2x, vb2y, vb2z);
vb3x = x[atom4][0] - x[atom3][0]; vb3x = x[atom4][0] - x[atom3][0];
vb3y = x[atom4][1] - x[atom3][1]; vb3y = x[atom4][1] - x[atom3][1];
vb3z = x[atom4][2] - x[atom3][2]; vb3z = x[atom4][2] - x[atom3][2];
domain->minimum_image(vb3x,vb3y,vb3z); domain->minimum_image(vb3x, vb3y, vb3z);
ss1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); ss1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z);
ss2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); ss2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z);
ss3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); ss3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z);
r1 = sqrt(ss1); r1 = sqrt(ss1);
r2 = sqrt(ss2); r2 = sqrt(ss2);
@ -202,20 +207,20 @@ int ComputeImproperLocal::compute_impropers(int flag)
c1 = (vb1x * vb2x + vb1y * vb2y + vb1z * vb2z) * r1 * r2; c1 = (vb1x * vb2x + vb1y * vb2y + vb1z * vb2z) * r1 * r2;
c2 = -(vb3x * vb2x + vb3y * vb2y + vb3z * vb2z) * r3 * r2; c2 = -(vb3x * vb2x + vb3y * vb2y + vb3z * vb2z) * r3 * r2;
s1 = 1.0 - c1*c1; s1 = 1.0 - c1 * c1;
if (s1 < SMALL) s1 = SMALL; if (s1 < SMALL) s1 = SMALL;
s1 = 1.0 / s1; s1 = 1.0 / s1;
s2 = 1.0 - c2*c2; s2 = 1.0 - c2 * c2;
if (s2 < SMALL) s2 = SMALL; if (s2 < SMALL) s2 = SMALL;
s2 = 1.0 / s2; s2 = 1.0 / s2;
s12 = sqrt(s1*s2); s12 = sqrt(s1 * s2);
c = (c1*c2 + c0) * s12; c = (c1 * c2 + c0) * s12;
if (c > 1.0) c = 1.0; if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0; if (c < -1.0) c = -1.0;
cbuf[n] = 180.0*acos(c)/MY_PI; cbuf[n] = 180.0 * acos(c) / MY_PI;
} }
n += nvalues; n += nvalues;
} }
@ -237,11 +242,11 @@ void ComputeImproperLocal::reallocate(int n)
if (nvalues == 1) { if (nvalues == 1) {
memory->destroy(vlocal); memory->destroy(vlocal);
memory->create(vlocal,nmax,"improper/local:vector_local"); memory->create(vlocal, nmax, "improper/local:vector_local");
vector_local = vlocal; vector_local = vlocal;
} else { } else {
memory->destroy(alocal); memory->destroy(alocal);
memory->create(alocal,nmax,nvalues,"improper/local:array_local"); memory->create(alocal, nmax, nvalues, "improper/local:array_local");
array_local = alocal; array_local = alocal;
} }
} }
@ -252,6 +257,6 @@ void ComputeImproperLocal::reallocate(int n)
double ComputeImproperLocal::memory_usage() double ComputeImproperLocal::memory_usage()
{ {
double bytes = (double)nmax*nvalues * sizeof(double); double bytes = (double) nmax * nvalues * sizeof(double);
return bytes; return bytes;
} }