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

change tagint to double storage

Browse Source
This commit is contained in:
Steve Plimpton
2022-04-19 11:49:55 -06:00
parent 44a4f3cff8
commit cc5bfd934d
14 changed files with 92 additions and 93 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/AMOEBA/amoeba_hal.cpp
View File

@ -152,8 +152,8 @@ void PairAmoeba::hal()
// increment the total van der Waals energy and derivatives
// if jv < 0, trigger an error, needed H-bond partner is missing
iv = ired2local[i];
jv = ired2local[j];
iv = red2local[i];
jv = red2local[j];
if (jv < 0)
error->one(FLERR,"AMOEBA hal cannot find H bond partner - "
"ghost comm is too short");

6
src/AMOEBA/amoeba_induce.cpp
View File

@ -772,8 +772,7 @@ void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp,
// owned atoms
pval = atom->dvector[index_pval];
double *pval = atom->dvector[index_pval];
double **x = atom->x;
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
@ -1441,9 +1440,8 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
// owned atoms
pval = atom->dvector[index_pval];
double **x = atom->x;
double *pval = atom->dvector[index_pval];
// neigh list

3
src/AMOEBA/amoeba_multipole.cpp
View File

@ -193,8 +193,7 @@ void PairAmoeba::multipole_real()
// owned atoms
pval = atom->dvector[index_pval];
double *pval = atom->dvector[index_pval];
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;

3
src/AMOEBA/amoeba_polar.cpp
View File

@ -275,8 +275,7 @@ void PairAmoeba::polar_real()
// owned atoms
pval = atom->dvector[index_pval];
double *pval = atom->dvector[index_pval];
double **x = atom->x;
double **f = atom->f;
int nlocal = atom->nlocal;

101
src/AMOEBA/amoeba_utils.cpp
View File

@ -51,11 +51,10 @@ void PairAmoeba::kmpole()
tagint angleneigh[36];
amtype = atom->ivector[index_amtype];
xaxis = atom->ivector[index_xaxis];
yaxis = atom->ivector[index_yaxis];
zaxis = atom->ivector[index_zaxis];
polaxe = atom->ivector[index_polaxe];
int *polaxe = atom->ivector[index_polaxe];
double *xaxis = atom->dvector[index_xaxis];
double *yaxis = atom->dvector[index_yaxis];
double *zaxis = atom->dvector[index_zaxis];
double **pole = fixpole->astore;
int **nspecial = atom->nspecial;
@ -126,9 +125,9 @@ void PairAmoeba::kmpole()
if (ktype == xtype) {
if (ytype == 0 && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = 0;
zaxis[i] = ubuf(jneigh).d;
xaxis[i] = ubuf(kneigh).d;
yaxis[i] = 0.0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
@ -142,9 +141,9 @@ void PairAmoeba::kmpole()
mtype = amtype[m];
if (mtype == ytype && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = mneigh;
zaxis[i] = ubuf(jneigh).d;
xaxis[i] = ubuf(kneigh).d;
yaxis[i] = ubuf(mneigh).d;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
@ -184,9 +183,9 @@ void PairAmoeba::kmpole()
if (ktype == xtype) {
if (ytype == 0 && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = 0;
zaxis[i] = ubuf(jneigh).d;
xaxis[i] = ubuf(kneigh).d;
yaxis[i] = 0.0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
@ -204,9 +203,9 @@ void PairAmoeba::kmpole()
if (!path) continue;
if (mtype == ytype && !flag) {
flag = 1;
zaxis[i] = jneigh;
xaxis[i] = kneigh;
yaxis[i] = mneigh;
zaxis[i] = ubuf(jneigh).d;
xaxis[i] = ubuf(kneigh).d;
yaxis[i] = ubuf(mneigh).d;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
@ -235,8 +234,8 @@ void PairAmoeba::kmpole()
if (jtype == ztype) {
if (xtype == 0 && !flag) {
flag = 1;
zaxis[i] = jtype;
xaxis[i] = yaxis[i] = 0;
zaxis[i] = ubuf(jneigh).d;
xaxis[i] = yaxis[i] = 0.0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
@ -255,7 +254,7 @@ void PairAmoeba::kmpole()
ztype = zpole[itype][iframe];
if (ztype == 0 && !flag) {
flag = 1;
zaxis[i] = xaxis[i] = yaxis[i] = 0;
zaxis[i] = xaxis[i] = yaxis[i] = 0.0;
polaxe[i] = mpaxis[itype][iframe];
for (j = 0; j < 13; j++)
pole[i][j] = fpole[itype][iframe][j];
@ -295,15 +294,18 @@ void PairAmoeba::chkpole(int i)
double xcd,ycd,zcd;
double c1,c2,c3,vol;
polaxe = atom->ivector[index_polaxe];
double **pole = fixpole->astore;
int *polaxe = atom->ivector[index_polaxe];
double *yaxis = atom->dvector[index_yaxis];
tagint yaxisID = (tagint) ubuf(yaxis[i]).i;
// test for chirality inversion
// if not, return
check = true;
if (polaxe[i] != ZTHENX) check = false;
if (yaxis[i] == 0) check = false;
if (yaxisID == 0) check = false;
if (!check) return;
ib = zaxis2local[i];
@ -331,9 +333,8 @@ void PairAmoeba::chkpole(int i)
// invert atomic multipole components involving the y-axis
// flip sign in permanent yaxis, not yaxis2local
int k = yaxis[i];
if ((k < 0 && vol > 0.0) || (k > 0 && vol < 0.0)) {
yaxis[i] = -k;
if ((yaxisID < 0 && vol > 0.0) || (yaxisID > 0 && vol < 0.0)) {
yaxis[i] = -ubuf(yaxisID).d;
pole[i][2] = -pole[i][2];
pole[i][5] = -pole[i][5];
pole[i][7] = -pole[i][7];
@ -359,7 +360,7 @@ void PairAmoeba::rotmat(int i)
double dx2,dy2,dz2;
double dx3,dy3,dz3;
polaxe = atom->ivector[index_polaxe];
int *polaxe = atom->ivector[index_polaxe];
// get coordinates and frame definition for the multipole site
@ -671,32 +672,34 @@ void PairAmoeba::add_onefive_neighbors()
/* ----------------------------------------------------------------------
update local indices of hydrogen neighbors for owned and ghost atoms
ired2local = used for offset of hydrogen positions in Vdwl term
red2local = used for offset of hydrogen positions in Vdwl term
called on reneighboring steps, only for AMOEBA
------------------------------------------------------------------------- */
void PairAmoeba::find_hydrogen_neighbors()
{
int index;
tagint id;
// grab current ptr for ired
// ired[i] = atom ID that atom I is bonded to
// ired2local[i] = local index of that atom
// grab current ptr for redID
// redID[i] = atom ID that atom I is bonded to
// red2local[i] = local index of that atom
// for furthest away ghost atoms, bond partner can be missing
// in that case ired2local = -1, but only an error if accessed in hal()
// in that case red2local = -1, but only an error if accessed in hal()
ired = atom->ivector[index_ired];
double *redID = atom->dvector[index_redID];
int nlocal = atom->nlocal;
int nall = nlocal + atom->nghost;
int nmissing = 0;
for (int i = 0; i < nall; i++) {
if (ired[i] == 0) ired2local[i] = i;
if (redID[i] == 0.0) red2local[i] = i;
else {
index = atom->map(ired[i]);
id = (tagint) ubuf(redID[i]).i;
index = atom->map(id);
if (index >= 0) index = domain->closest_image(i,index);
ired2local[i] = index;
red2local[i] = index;
}
}
}
@ -710,22 +713,27 @@ void PairAmoeba::find_hydrogen_neighbors()
void PairAmoeba::find_multipole_neighbors()
{
int index,ypos;
tagint xaxisID,yaxisID,zaxisID;
// grab current pts for xaxis,yaxis,zaxis
// xyz axis[i] = atom IDs that atom I uses for its multipole orientation
// can be zero if not used, in which case set local index to self
// yaxis can be negative, in which case use absolute value
xaxis = atom->ivector[index_xaxis];
yaxis = atom->ivector[index_yaxis];
zaxis = atom->ivector[index_zaxis];
double *xaxis = atom->dvector[index_xaxis];
double *yaxis = atom->dvector[index_yaxis];
double *zaxis = atom->dvector[index_zaxis];
int nlocal = atom->nlocal;
int nmissing = 0;
for (int i = 0; i < nlocal; i++) {
if (xaxis[i]) {
index = atom->map(xaxis[i]);
xaxisID = (tagint) ubuf(xaxis[i]).i;
yaxisID = (tagint) ubuf(yaxis[i]).i;
zaxisID = (tagint) ubuf(zaxis[i]).i;
if (xaxisID) {
index = atom->map(xaxisID);
if (index == -1) nmissing++;
else {
index = domain->closest_image(i,index);
@ -733,10 +741,9 @@ void PairAmoeba::find_multipole_neighbors()
}
} else xaxis2local[i] = i;
if (yaxis[i]) {
if (yaxis[i] > 0) ypos = yaxis[i];
else ypos = -yaxis[i];
index = atom->map(ypos);
if (yaxisID) {
if (yaxis[i] < 0) yaxisID = -yaxisID;
index = atom->map(yaxisID);
if (index == -1) nmissing++;
else {
index = domain->closest_image(i,index);
@ -744,8 +751,8 @@ void PairAmoeba::find_multipole_neighbors()
}
} else yaxis2local[i] = i;
if (zaxis[i]) {
index = atom->map(zaxis[i]);
if (zaxisID) {
index = atom->map(zaxisID);
if (index == -1) nmissing++;
else {
index = domain->closest_image(i,index);
@ -838,7 +845,7 @@ void PairAmoeba::torque2force(int i, double *trq,
// get the local frame type and the frame-defining atoms
polaxe = atom->ivector[index_polaxe];
int *polaxe = atom->ivector[index_polaxe];
axetyp = polaxe[i];
if (axetyp == NOFRAME) return;

43
src/AMOEBA/pair_amoeba.cpp
View File

@ -84,7 +84,7 @@ PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
rpole = NULL;
tq = NULL;
ired2local = NULL;
red2local = NULL;
xred = NULL;
uind = uinp = udirp = NULL;
@ -166,7 +166,7 @@ PairAmoeba::~PairAmoeba()
memory->destroy(rpole);
memory->destroy(tq);
memory->destroy(ired2local);
memory->destroy(red2local);
memory->destroy(xred);
memory->destroy(uind);
@ -289,7 +289,7 @@ void PairAmoeba::compute(int eflag, int vflag)
kmpole();
if (hippo) {
pval = atom->dvector[index_pval];
double *pval = atom->dvector[index_pval];
double **pole = fixpole->astore;
int nlocal = atom->nlocal;
int itype,iclass;
@ -316,7 +316,7 @@ void PairAmoeba::compute(int eflag, int vflag)
// if reneighboring step:
// augment neighbor list to include 1-5 neighbor flags
// re-create xyz axis2local and ired2local
// re-create xyz axis2local and red2local
// re-create induce neighbor list
if (neighbor->ago == 0) {
@ -341,7 +341,7 @@ void PairAmoeba::compute(int eflag, int vflag)
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
j = ired2local[i];
j = red2local[i];
iclass = amtype2class[amtype[i]];
rdn = kred[iclass];
xred[i][0] = rdn*(x[i][0]-x[j][0]) + x[j][0];
@ -707,17 +707,17 @@ void PairAmoeba::init_style()
if (index_amgroup < 0 || flag || cols)
error->all(FLERR,"Pair amoeba amgroup is not defined");
index_ired = atom->find_custom("ired",flag,cols);
if (index_ired < 0 || flag || cols)
error->all(FLERR,"Pair amoeba ired is not defined");
index_redID = atom->find_custom("redID",flag,cols);
if (index_redID < 0 || !flag || cols)
error->all(FLERR,"Pair amoeba redID is not defined");
index_xaxis = atom->find_custom("xaxis",flag,cols);
if (index_xaxis < 0 || flag || cols)
if (index_xaxis < 0 || !flag || cols)
error->all(FLERR,"Pair amoeba xaxis is not defined");
index_yaxis = atom->find_custom("yaxis",flag,cols);
if (index_yaxis < 0 || flag || cols)
if (index_yaxis < 0 || !flag || cols)
error->all(FLERR,"Pair amoeba yaxis is not defined");
index_zaxis = atom->find_custom("zaxis",flag,cols);
if (index_zaxis < 0 || flag || cols)
if (index_zaxis < 0 || !flag || cols)
error->all(FLERR,"Pair amoeba zaxis is not defined");
index_polaxe = atom->find_custom("polaxe",flag,cols);
@ -879,7 +879,7 @@ void PairAmoeba::init_style()
// pval is not set until first call to compute(), and only for HIPPO
if (first_flag) {
pval = atom->dvector[index_pval];
double *pval = atom->dvector[index_pval];
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) pval[i] = 0.0;
}
@ -918,14 +918,13 @@ void PairAmoeba::init_style()
fixupalt = (FixStore *) modify->fix[ifix];
}
// assign hydrogen neighbors (ired) to each owned atom
// assign hydrogen neighbors (redID) to each owned atom
// only set if kred[i] is non-zero and I is bonded to a single atom
// conceptually: non-zero if I is a hydrogen bonded to another atom
// NOTE: ired needs to be a tagint vector? but atom ivector is not!
// conceptually: non-zero if I is hydrogen bonded to another atom
if (amoeba) {
amtype = atom->ivector[index_amtype];
ired = atom->ivector[index_ired];
double *redID = atom->dvector[index_redID];
int **nspecial = atom->nspecial;
tagint **special = atom->special;
@ -937,13 +936,13 @@ void PairAmoeba::init_style()
itype = amtype[i];
iclass = amtype2class[itype];
if (kred[iclass] == 0.0) {
ired[i] = 0;
redID[i] = 0.0;
}
else if (nspecial[i][0] != 1) {
ired[i] = 0;
redID[i] = 0.0;
}
else {
ired[i] = special[i][0];
redID[i] = ubuf(special[i][0]).d;
}
}
}
@ -1207,6 +1206,7 @@ int PairAmoeba::pack_forward_comm(int n, int *list, double *buf,
}
} else if (cfstyle == PVAL) {
double *pval = atom->dvector[index_pval];
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = pval[j];
@ -1310,6 +1310,7 @@ void PairAmoeba::unpack_forward_comm(int n, int first, double *buf)
}
} else if (cfstyle == PVAL) {
double *pval = atom->dvector[index_pval];
for (i = first; i < last; i++) {
pval[i] = buf[m++];
}
@ -2112,7 +2113,7 @@ void PairAmoeba::grow_local()
memory->destroy(tq);
if (amoeba) {
memory->destroy(ired2local);
memory->destroy(red2local);
memory->destroy(xred);
}
@ -2163,7 +2164,7 @@ void PairAmoeba::grow_local()
memory->create(tq,nmax,3,"amoeba:tq");
if (amoeba) {
memory->create(ired2local,nmax,"amoeba:ired2local");
memory->create(red2local,nmax,"amoeba:red2local");
memory->create(xred,nmax,3,"amoeba:xred");
}

11
src/AMOEBA/pair_amoeba.h
View File

@ -158,16 +158,11 @@ class PairAmoeba : public Pair {
// static per-atom properties, must persist as atoms migrate
int index_amtype,index_amgroup,index_ired;
int index_xaxis,index_yaxis,index_zaxis;
int index_polaxe,index_pval;
int index_amtype,index_amgroup,index_redID;
int index_xaxis,index_yaxis,index_zaxis,index_polaxe,index_pval;
int *amtype; // AMOEBA type, 1 to N_amtype
int *amgroup; // AMOEBA polarization group, 1 to Ngroup
tagint *ired; // ID of atom that H is bonded to
tagint *xaxis,*yaxis,*zaxis; // IDs of nearby atoms for multipole def
int *polaxe;
double *pval;
char *id_pole,*id_udalt,*id_upalt;
class FixStore *fixpole; // stores pole = multipole components
@ -252,7 +247,7 @@ class PairAmoeba : public Pair {
// just for owned atoms
// set to self if not defined
int *ired2local; // local indices of ired IDs, computed for owned and ghost
int *red2local; // local indices of ired IDs, computed for owned and ghost
double **xred; // altered coords for H atoms for Vdwl, comm to ghosts
double **tq; // torque from pairwise multipole, reverse comm from ghosts