diff --git a/examples/amoeba/in.ubiquitin b/examples/amoeba/in.ubiquitin index 9f3fefba73..ef6f55c6da 100644 --- a/examples/amoeba/in.ubiquitin +++ b/examples/amoeba/in.ubiquitin @@ -20,7 +20,7 @@ fix bit all amoeba/bitorsion bitorsion.ubiquitin.data fix_modify bit energy yes fix extra all property/atom & - i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes + i_amgroup d_redID d_xaxis d_yaxis d_zaxis d_pval ghost yes fix extra2 all property/atom i_polaxe #read_data data.ubiquitin fix amtype NULL "Tinker Types" diff --git a/examples/amoeba/in.water_box.amoeba b/examples/amoeba/in.water_box.amoeba index 3a37784f9d..1e1a72463e 100644 --- a/examples/amoeba/in.water_box.amoeba +++ b/examples/amoeba/in.water_box.amoeba @@ -12,7 +12,7 @@ dihedral_style none fix amtype all property/atom i_amtype ghost yes fix extra all property/atom & - i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes + i_amgroup d_redID d_xaxis d_yaxis d_zaxis d_pval ghost yes fix extra2 all property/atom i_polaxe # read data file diff --git a/examples/amoeba/in.water_box.hippo b/examples/amoeba/in.water_box.hippo index 8e1d9a8e70..49acb58317 100644 --- a/examples/amoeba/in.water_box.hippo +++ b/examples/amoeba/in.water_box.hippo @@ -12,7 +12,7 @@ dihedral_style none fix amtype all property/atom i_amtype ghost yes fix extra all property/atom & - i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes + i_amgroup d_redID d_xaxis d_yaxis d_zaxis d_pval ghost yes fix extra2 all property/atom i_polaxe # read data file diff --git a/examples/amoeba/in.water_dimer.amoeba b/examples/amoeba/in.water_dimer.amoeba index c0e6ee35e4..9e9cae8ac3 100644 --- a/examples/amoeba/in.water_dimer.amoeba +++ b/examples/amoeba/in.water_dimer.amoeba @@ -12,7 +12,7 @@ dihedral_style none fix amtype all property/atom i_amtype ghost yes fix extra all property/atom & - i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes + i_amgroup d_redID d_xaxis d_yaxis d_zaxis d_pval ghost yes fix extra2 all property/atom i_polaxe # read data file diff --git a/examples/amoeba/in.water_dimer.hippo b/examples/amoeba/in.water_dimer.hippo index 55b173b789..1e9755570a 100644 --- a/examples/amoeba/in.water_dimer.hippo +++ b/examples/amoeba/in.water_dimer.hippo @@ -12,7 +12,7 @@ dihedral_style none fix amtype all property/atom i_amtype ghost yes fix extra all property/atom & - i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes + i_amgroup d_redID d_xaxis d_yaxis d_zaxis d_pval ghost yes fix extra2 all property/atom i_polaxe # read data file diff --git a/examples/amoeba/in.water_hexamer.amoeba b/examples/amoeba/in.water_hexamer.amoeba index eb16206f45..0ce596a06c 100644 --- a/examples/amoeba/in.water_hexamer.amoeba +++ b/examples/amoeba/in.water_hexamer.amoeba @@ -12,7 +12,7 @@ dihedral_style none fix amtype all property/atom i_amtype ghost yes fix extra all property/atom & - i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes + i_amgroup d_redID d_xaxis d_yaxis d_zaxis d_pval ghost yes fix extra2 all property/atom i_polaxe # read data file diff --git a/examples/amoeba/in.water_hexamer.hippo b/examples/amoeba/in.water_hexamer.hippo index 1287a1f527..52f6e863a2 100644 --- a/examples/amoeba/in.water_hexamer.hippo +++ b/examples/amoeba/in.water_hexamer.hippo @@ -12,7 +12,7 @@ dihedral_style none fix amtype all property/atom i_amtype ghost yes fix extra all property/atom & - i_amgroup i_ired i_xaxis i_yaxis i_zaxis d_pval ghost yes + i_amgroup d_redID d_xaxis d_yaxis d_zaxis d_pval ghost yes fix extra2 all property/atom i_polaxe # read data file diff --git a/src/AMOEBA/amoeba_hal.cpp b/src/AMOEBA/amoeba_hal.cpp index f402221520..01123251e7 100644 --- a/src/AMOEBA/amoeba_hal.cpp +++ b/src/AMOEBA/amoeba_hal.cpp @@ -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"); diff --git a/src/AMOEBA/amoeba_induce.cpp b/src/AMOEBA/amoeba_induce.cpp index beac1186f7..cbd0a2c480 100644 --- a/src/AMOEBA/amoeba_induce.cpp +++ b/src/AMOEBA/amoeba_induce.cpp @@ -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 diff --git a/src/AMOEBA/amoeba_multipole.cpp b/src/AMOEBA/amoeba_multipole.cpp index d0a59a7a0b..5f66f41deb 100644 --- a/src/AMOEBA/amoeba_multipole.cpp +++ b/src/AMOEBA/amoeba_multipole.cpp @@ -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; diff --git a/src/AMOEBA/amoeba_polar.cpp b/src/AMOEBA/amoeba_polar.cpp index 0e9b5250f3..da84053f0d 100644 --- a/src/AMOEBA/amoeba_polar.cpp +++ b/src/AMOEBA/amoeba_polar.cpp @@ -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; diff --git a/src/AMOEBA/amoeba_utils.cpp b/src/AMOEBA/amoeba_utils.cpp index e7eb99cced..b1be3ef30b 100644 --- a/src/AMOEBA/amoeba_utils.cpp +++ b/src/AMOEBA/amoeba_utils.cpp @@ -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; diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp index 17ea496480..6869aa57e4 100644 --- a/src/AMOEBA/pair_amoeba.cpp +++ b/src/AMOEBA/pair_amoeba.cpp @@ -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"); } diff --git a/src/AMOEBA/pair_amoeba.h b/src/AMOEBA/pair_amoeba.h index 15ba4d8526..0057219b98 100644 --- a/src/AMOEBA/pair_amoeba.h +++ b/src/AMOEBA/pair_amoeba.h @@ -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