From 4c19eab64ca88c9a7edeccd724c7a0cf8f39dc40 Mon Sep 17 00:00:00 2001 From: Mingjian Wen Date: Tue, 23 Apr 2019 14:39:47 -0500 Subject: [PATCH] Bugfix no 3 nearest neighbors for ghost atoms near boundary --- potentials/C.drip | 10 ++++++--- src/USER-MISC/pair_drip.cpp | 44 +++++++++++++++++++++++++++++-------- src/USER-MISC/pair_drip.h | 4 ++-- 3 files changed, 44 insertions(+), 14 deletions(-) diff --git a/potentials/C.drip b/potentials/C.drip index 43d5ca4208..435efe40a9 100644 --- a/potentials/C.drip +++ b/potentials/C.drip @@ -6,10 +6,14 @@ # Cite as M. Wen, S. Carr, S. Fang, E. Kaxiras, and E. B. Tadmor, Phys. Rev. B, 98, 235404 (2018). -# C0 C2 C4 C delta lambda A z0 B eta rho_cut r_cut -C C 1.1598e-02 1.2981e-02 3.2515e-02 7.8151e-03 8.3679e-01 2.7158 2.2216e-02 3.34 7.6799e-03 1.1432 1.562 12.0 +# C0 C2 C4 C delta lambda A z0 B eta rho_cut r_cut normal_cut +C C 1.1598e-02 1.2981e-02 3.2515e-02 7.8151e-03 8.3679e-01 2.7158 2.2216e-02 3.34 7.6799e-03 1.1432 1.562 12.0 3.7 # C0, C2, C4, C, A, and B in [eV] -# delta, z0, eta, rho_cut, and r_cut in [Angstrom] +# delta, z0, eta, rho_cut, r_cut, and normal_cut in [Angstrom] # lambda in [1/Angstrom] +# +# normal_cut is a parameter not present in the Wen paper, but specific to the +# LAMMPS implementation, which is used to find the 3 nearest neighbors of an +# atom to construct the normal. diff --git a/src/USER-MISC/pair_drip.cpp b/src/USER-MISC/pair_drip.cpp index df7f39e8cc..338891c77a 100644 --- a/src/USER-MISC/pair_drip.cpp +++ b/src/USER-MISC/pair_drip.cpp @@ -197,6 +197,14 @@ double PairDRIP::init_one(int i, int j) { if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + int itype = map[i]; + int jtype = map[j]; + int iparam_ij = elem2param[itype][jtype]; + Param& p = params[iparam_ij]; + + // max cutoff is the main cutoff plus the normal cutoff such that + double cutmax = p.rcut + p.ncut; + return cutmax; } @@ -206,7 +214,7 @@ double PairDRIP::init_one(int i, int j) void PairDRIP::read_file(char *filename) { - int params_per_line = 14; + int params_per_line = 15; char **words = new char*[params_per_line+1]; memory->sfree(params); int nparams = 0; @@ -311,13 +319,12 @@ void PairDRIP::read_file(char *filename) params[nparams].eta = atof(words[11]); params[nparams].rhocut = atof(words[12]); params[nparams].rcut = atof(words[13]); + params[nparams].ncut = atof(words[14]); // convenient precomputations params[nparams].rhocutsq = params[nparams].rhocut * params[nparams].rhocut; params[nparams].rcutsq = params[nparams].rcut * params[nparams].rcut; - - // set max cutoff - if(params[nparams].rcut > cutmax) cutmax = params[nparams].rcut; + params[nparams].ncutsq = params[nparams].ncut * params[nparams].ncut; nparams++; } @@ -371,6 +378,9 @@ void PairDRIP::compute(int eflag, int vflag) for (ii = 0; ii < inum; ii++) { i = ilist[ii]; + if (nearest3neigh[i][0] == -1) { + continue; + } itag = tag[i]; xtmp = x[i][0]; ytmp = x[i][1]; @@ -387,6 +397,9 @@ void PairDRIP::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; j &= NEIGHMASK; + if (nearest3neigh[j][0] == -1) { + continue; + } jtype = map[type[j]]; jtag = tag[j]; @@ -604,7 +617,7 @@ double PairDRIP::calc_repulsive(int const i, int const j, Param& p, void PairDRIP::find_nearest3neigh() { - int i, j, ii, jj, n, allnum, jnum, itype, jtype, size; + int i, j, ii, jj, n, allnum, inum, jnum, itype, jtype, size; double xtmp, ytmp, ztmp, delx, dely, delz, rsq; int *ilist, *jlist, *numneigh, **firstneigh; @@ -613,6 +626,7 @@ void PairDRIP::find_nearest3neigh() allnum = list->inum + list->gnum; + inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; @@ -647,6 +661,7 @@ void PairDRIP::find_nearest3neigh() double nb3_rsq = 3.0e10; for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; j &= NEIGHMASK; jtype = map[type[j]]; @@ -656,9 +671,9 @@ void PairDRIP::find_nearest3neigh() rsq = delx * delx + dely * dely + delz * delz; int iparam_ij = elem2param[itype][jtype]; - double rcutsq = params[iparam_ij].rcutsq; + double ncutsq = params[iparam_ij].ncutsq; - if (rsq < rcutsq && atom->molecule[i] == atom->molecule[j]) { + if (rsq < ncutsq && atom->molecule[i] == atom->molecule[j]) { // find the 3 nearest neigh if (rsq < nb1_rsq) { nb3 = nb2; @@ -683,8 +698,19 @@ void PairDRIP::find_nearest3neigh() } // loop over jj // store neighbors to be used later to compute normal - if (nb1_rsq >= 1.0e10 || nb2_rsq >= 1.0e10 || nb3_rsq >= 1.0e10) { - error->one(FLERR, "No enough neighbors to construct normal."); + if (nb3_rsq >= 1.0e10) { + if (ione(FLERR, "No enough neighbors to construct normal."); + } else { + // This only happens for ghost atoms that are near the boundary of the + // domain (i.e. r > r_cut + n_cut). These ghost atoms will not be + // the i j atoms in the compute function, but only neighbors of j atoms. + // It is allowed not to have three neighbors for these atoms, since + // their normals are not needed. + nearest3neigh[i][0] = -1; + nearest3neigh[i][1] = -1; + nearest3neigh[i][2] = -1; + } } else{ nearest3neigh[i][0] = nb1; diff --git a/src/USER-MISC/pair_drip.h b/src/USER-MISC/pair_drip.h index b229f1aced..3035d88ad8 100644 --- a/src/USER-MISC/pair_drip.h +++ b/src/USER-MISC/pair_drip.h @@ -55,8 +55,8 @@ protected: struct Param { int ielement, jelement; - double C0, C2, C4, C, delta, lambda, A, z0, B, eta, rhocut, rcut; - double rhocutsq, rcutsq; + double C0, C2, C4, C, delta, lambda, A, z0, B, eta, rhocut, rcut, ncut; + double rhocutsq, rcutsq, ncutsq; }; Param *params; // parameter set for I-J interactions int **nearest3neigh; // nearest 3 neighbors of atoms