diff --git a/source/lmp/fix_dplr.cpp b/source/lmp/fix_dplr.cpp index e66cb643f9..8ee0344d5c 100644 --- a/source/lmp/fix_dplr.cpp +++ b/source/lmp/fix_dplr.cpp @@ -175,11 +175,9 @@ void FixDPLR::get_valid_pairs(vector > &pairs) { int nall = nlocal + nghost; vector dtype(nall); // get type - { - int *type = atom->type; - for (int ii = 0; ii < nall; ++ii) { - dtype[ii] = type[ii] - 1; - } + int *type = atom->type; + for (int ii = 0; ii < nall; ++ii) { + dtype[ii] = type[ii] - 1; } int **bondlist = neighbor->bondlist; @@ -190,21 +188,51 @@ void FixDPLR::get_valid_pairs(vector > &pairs) { if (!binary_search(bond_type.begin(), bond_type.end(), bd_type)) { continue; } - if (binary_search(sel_type.begin(), sel_type.end(), - dtype[bondlist[ii][0]]) && - binary_search(dpl_type.begin(), dpl_type.end(), - dtype[bondlist[ii][1]])) { - idx0 = bondlist[ii][0]; - idx1 = bondlist[ii][1]; - } else if (binary_search(sel_type.begin(), sel_type.end(), - dtype[bondlist[ii][1]]) && - binary_search(dpl_type.begin(), dpl_type.end(), - dtype[bondlist[ii][0]])) { - idx0 = bondlist[ii][1]; - idx1 = bondlist[ii][0]; + std::vector::iterator it = + find(sel_type.begin(), sel_type.end(), dtype[bondlist[ii][0]]); + if (it != sel_type.end()) { + int idx_type = distance(sel_type.begin(), it); + if (dtype[bondlist[ii][1]] == dpl_type[idx_type]) { + idx0 = bondlist[ii][0]; + idx1 = bondlist[ii][1]; + } else { + char str[300]; + sprintf(str, + "Invalid pair: %d %d \n A virtual atom of type %d is " + "expected, but the type of atom %d is " + "%d.\n Please check your data file carefully.\n", + atom->tag[bondlist[ii][0]], atom->tag[bondlist[ii][1]], + dpl_type[idx_type] + 1, atom->tag[bondlist[ii][1]], + type[bondlist[ii][1]]); + error->all(FLERR, str); + } } else { - error->all(FLERR, - "find a bonded pair the types of which are not associated"); + it = find(sel_type.begin(), sel_type.end(), dtype[bondlist[ii][1]]); + if (it != sel_type.end()) { + int idx_type = distance(sel_type.begin(), it); + if (dtype[bondlist[ii][0]] == dpl_type[idx_type]) { + idx0 = bondlist[ii][1]; + idx1 = bondlist[ii][0]; + } else { + char str[300]; + sprintf(str, + "Invalid pair: %d %d \n A virtual atom of type %d is " + "expected, but the type of atom %d is %d.\n Please " + "check your data file carefully.\n", + atom->tag[bondlist[ii][0]], atom->tag[bondlist[ii][1]], + dpl_type[idx_type] + 1, atom->tag[bondlist[ii][0]], + type[bondlist[ii][0]]); + error->all(FLERR, str); + } + } else { + char str[300]; + sprintf(str, + "Invalid pair: %d %d \n They are not expected to have " + "Wannier centroid.\n Please check your data file " + "carefully.\n", + atom->tag[bondlist[ii][0]], atom->tag[bondlist[ii][1]]); + error->all(FLERR, str); + } } if (!(idx0 < nlocal && idx1 < nlocal)) { error->all(FLERR, @@ -306,11 +334,6 @@ void FixDPLR::pre_force(int vflag) { // } // } - // selected type - vector dpl_type; - for (int ii = 0; ii < sel_type.size(); ++ii) { - dpl_type.push_back(type_asso[sel_type[ii]]); - } vector sel_fwd, sel_bwd; int sel_nghost; deepmd_compat::select_by_type(sel_fwd, sel_bwd, sel_nghost, dcoord, dtype, @@ -368,10 +391,6 @@ void FixDPLR::post_force(int vflag) { } PPPMDPLR *pppm_dplr = (PPPMDPLR *)force->kspace_match("pppm/dplr", 1); - if (!pppm_dplr) { - error->all(FLERR, "kspace_style pppm/dplr should be set before this fix\n"); - } - const vector &dfele_(pppm_dplr->get_fele()); int nlocal = atom->nlocal; int nghost = atom->nghost; int nall = nlocal + nghost; @@ -397,10 +416,13 @@ void FixDPLR::post_force(int vflag) { dcoord[ii * 3 + dd] = x[ii][dd] - domain->boxlo[dd]; } } - assert(dfele_.size() == nlocal * 3); // revise force according to efield - for (int ii = 0; ii < nlocal * 3; ++ii) { - dfele[ii] = dfele_[ii]; + if (pppm_dplr) { + const vector &dfele_(pppm_dplr->get_fele()); + assert(dfele_.size() == nlocal * 3); + for (int ii = 0; ii < nlocal * 3; ++ii) { + dfele[ii] += dfele_[ii]; + } } // revise force and virial according to efield double *q = atom->q;