diff --git a/src/USER-OMP/npair_full_multi_omp.cpp b/src/USER-OMP/npair_full_multi_omp.cpp index cf897f9dea..11b8bae196 100755 --- a/src/USER-OMP/npair_full_multi_omp.cpp +++ b/src/USER-OMP/npair_full_multi_omp.cpp @@ -30,7 +30,7 @@ NPairFullMultiOmp::NPairFullMultiOmp(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- binned neighbor list construction for all neighbors - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent every neighbor pair appears in list of both atoms i and j ------------------------------------------------------------------------- */ @@ -46,7 +46,7 @@ void NPairFullMultiOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,jbin,which,ns,imol,iatom; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,which,ns,imol,iatom; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; @@ -80,6 +80,7 @@ void NPairFullMultiOmp::build(NeighList *list) neighptr = ipage.vget(); itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -89,27 +90,28 @@ void NPairFullMultiOmp::build(NeighList *list) tagprev = tag[i] - iatom - 1; } - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); // loop over all atoms in surrounding bins in stencil including self // skip i = j - // use full stencil for all type combinations + // use full stencil for all group combinations - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { if (i == j) continue; + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; diff --git a/src/USER-OMP/npair_half_multi_newtoff_omp.cpp b/src/USER-OMP/npair_half_multi_newtoff_omp.cpp index bf7644f01f..e3e03ee8ed 100755 --- a/src/USER-OMP/npair_half_multi_newtoff_omp.cpp +++ b/src/USER-OMP/npair_half_multi_newtoff_omp.cpp @@ -30,7 +30,7 @@ NPairHalfMultiNewtoffOmp::NPairHalfMultiNewtoffOmp(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- binned neighbor list construction with partial Newton's 3rd law - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks own bin and other bins in stencil pair stored once if i,j are both owned and i < j pair stored by me if j is ghost (also stored by proc owning j) @@ -48,7 +48,7 @@ void NPairHalfMultiNewtoffOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,jbin,which,ns,imol,iatom; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,which,ns,imol,iatom; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; @@ -82,6 +82,7 @@ void NPairHalfMultiNewtoffOmp::build(NeighList *list) neighptr = ipage.vget(); itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -91,29 +92,30 @@ void NPairHalfMultiNewtoffOmp::build(NeighList *list) tagprev = tag[i] - iatom - 1; } - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); // loop over all atoms in other bins in stencil including self // only store pair if i < j // stores own/own pairs only once // stores own/ghost pairs on both procs - // use full stencil for all type combinations + // use full stencil for all group combinations - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >=0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >=0; j = bins_multi[jgroup][j]) { if (j <= i) continue; + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; diff --git a/src/USER-OMP/npair_half_multi_newton_omp.cpp b/src/USER-OMP/npair_half_multi_newton_omp.cpp index c68a3a150f..28db1b3c0d 100755 --- a/src/USER-OMP/npair_half_multi_newton_omp.cpp +++ b/src/USER-OMP/npair_half_multi_newton_omp.cpp @@ -30,7 +30,7 @@ NPairHalfMultiNewtonOmp::NPairHalfMultiNewtonOmp(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- binned neighbor list construction with full Newton's 3rd law - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks its own bin and other bins in Newton stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ @@ -47,7 +47,7 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,jbin,which,ns,imol,iatom; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,which,ns,imol,iatom; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; @@ -81,6 +81,7 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list) neighptr = ipage.vget(); itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -90,28 +91,28 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list) tagprev = tag[i] - iatom - 1; } - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); - if(cutneighsq[itype][itype] == cutneighsq[jtype][jtype]){ + if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){ // if same size: use half stencil - if(itype == jtype){ + if(igroup == jgroup){ - // if same type, implement with: + // if same group, implement with: // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i - js = bins_multi[itype][i]; + js = bins_multi[igroup][i]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { if (j >= nlocal) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { @@ -120,6 +121,7 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list) } } + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; @@ -145,14 +147,14 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list) } } else { - // if different types, implement with: - // loop over all atoms in jtype bin + // if different groups, implement with: + // loop over all atoms in jgroup bin // if j is owned atom, store it if j > i // if j is ghost, only store if j coords are "above and to the right" of i - js = binhead_multi[jtype][jbin]; + js = binhead_multi[jgroup][jbin]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { if(j < i) continue; if (j >= nlocal) { @@ -163,8 +165,9 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list) } } + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; @@ -189,20 +192,21 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list) } } - // for all types, loop over all atoms in other bins in stencil, store every pair + // for all groups, loop over all atoms in other bins in stencil, store every pair // stencil is empty if i larger than j // stencil is half if i same size as j // stencil is full if i smaller than j - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { - - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; diff --git a/src/USER-OMP/npair_half_multi_newton_tri_omp.cpp b/src/USER-OMP/npair_half_multi_newton_tri_omp.cpp index 3f3ac37b6d..0289967d30 100755 --- a/src/USER-OMP/npair_half_multi_newton_tri_omp.cpp +++ b/src/USER-OMP/npair_half_multi_newton_tri_omp.cpp @@ -31,7 +31,7 @@ NPairHalfMultiNewtonTriOmp::NPairHalfMultiNewtonTriOmp(LAMMPS *lmp) : /* ---------------------------------------------------------------------- binned neighbor list construction with Newton's 3rd law for triclinic - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks its own bin and other bins in triclinic stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ @@ -48,7 +48,7 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,jbin,which,ns,imol,iatom; + int i,j,k,n,itype,jtype,ibin,jbin,igroup,jgroup,which,ns,imol,iatom; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; @@ -82,6 +82,7 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list) neighptr = ipage.vget(); itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -91,14 +92,14 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list) tagprev = tag[i] - iatom - 1; } - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); // loop over all atoms in bins in stencil // stencil is empty if i larger than j @@ -109,15 +110,15 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list) // (equal zyx and j <= i) // latter excludes self-self interaction but allows superposed atoms - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { - // if same size (e.g. same type), use half stencil - if(cutneighsq[itype][itype] == cutneighsq[jtype][jtype]){ + // if same size (same group), use half stencil + if(cutmultisq[igroup][igroup] == cutmutlisq[jgroup][jgroup]){ if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; @@ -128,6 +129,7 @@ void NPairHalfMultiNewtonTriOmp::build(NeighList *list) } } + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; diff --git a/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp b/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp index 8201ba8cb6..35cfaad5e6 100755 --- a/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp +++ b/src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp @@ -29,7 +29,7 @@ NPairHalfSizeMultiNewtoffOmp::NPairHalfSizeMultiNewtoffOmp(LAMMPS *lmp) : NPair( /* ---------------------------------------------------------------------- size particles binned neighbor list construction with partial Newton's 3rd law - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks own bin and other bins in stencil pair stored once if i,j are both owned and i < j pair stored by me if j is ghost (also stored by proc owning j) @@ -47,7 +47,7 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,jbin,ns; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,ns; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -75,34 +75,36 @@ void NPairHalfSizeMultiNewtoffOmp::build(NeighList *list) neighptr = ipage.vget(); itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); // loop over all atoms in other bins in stencil including self // only store pair if i < j // stores own/own pairs only once // stores own/ghost pairs on both procs - // use full stencil for all type combinations + // use full stencil for all group combinations - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >=0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >=0; j = bins_multi[jgroup][j]) { if (j <= i) continue; + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; diff --git a/src/USER-OMP/npair_half_size_multi_newton_omp.cpp b/src/USER-OMP/npair_half_size_multi_newton_omp.cpp index 197b41c12e..974b8c5809 100755 --- a/src/USER-OMP/npair_half_size_multi_newton_omp.cpp +++ b/src/USER-OMP/npair_half_size_multi_newton_omp.cpp @@ -29,7 +29,7 @@ NPairHalfSizeMultiNewtonOmp::NPairHalfSizeMultiNewtonOmp(LAMMPS *lmp) : NPair(lm /* ---------------------------------------------------------------------- size particles binned neighbor list construction with full Newton's 3rd law - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks its own bin and other bins in Newton stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ @@ -46,7 +46,7 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,jbin,ns; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,ns; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -74,33 +74,34 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) neighptr = ipage.vget(); itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_group; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); - if(cutneighsq[itype][itype] == cutneighsq[jtype][jtype]){ + if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){ // if same size: use half stencil - if(itype == jtype){ + if(igroup == jgroup){ - // if same type, implement with: + // if same group, implement with: // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i - js = bins_multi[itype][i]; + js = bins_multi[igroup][i]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { if (j >= nlocal) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { @@ -109,6 +110,7 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) } } + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; @@ -127,14 +129,14 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) } } else { - // if different types, implement with: - // loop over all atoms in jtype bin + // if different groups, implement with: + // loop over all atoms in jgroup bin // if j is owned atom, store it if j > i // if j is ghost, only store if j coords are "above and to the right" of i - js = binhead_multi[jtype][jbin]; + js = binhead_multi[jgroup][jbin]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { if(j < i) continue; if (j >= nlocal) { @@ -145,6 +147,7 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) } } + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; @@ -164,18 +167,19 @@ void NPairHalfSizeMultiNewtonOmp::build(NeighList *list) } } - // for all types, loop over all atoms in other bins in stencil, store every pair + // for all groups, loop over all atoms in other bins in stencil, store every pair // stencil is empty if i larger than j // stencil is half if i same size as j // stencil is full if i smaller than j - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; diff --git a/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp index cd123202b1..d1c5a97498 100755 --- a/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp +++ b/src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp @@ -30,7 +30,7 @@ NPairHalfSizeMultiNewtonTriOmp::NPairHalfSizeMultiNewtonTriOmp(LAMMPS *lmp) : /* ---------------------------------------------------------------------- size particles binned neighbor list construction with Newton's 3rd law for triclinic - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks its own bin and other bins in triclinic stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ @@ -47,7 +47,7 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) #endif NPAIR_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,jbin,ns; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,ns; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -75,19 +75,20 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) neighptr = ipage.vget(); itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); // loop over all atoms in bins in stencil @@ -99,15 +100,15 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) // (equal zyx and j <= i) // latter excludes self-self interaction but allows superposed atoms - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { - // if same size (e.g. same type), use half stencil - if(cutneighsq[itype][itype] == cutneighsq[jtype][jtype]){ + // if same size (same group), use half stencil + if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){ if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; @@ -118,6 +119,7 @@ void NPairHalfSizeMultiNewtonTriOmp::build(NeighList *list) } } + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; diff --git a/src/comm_brick.cpp b/src/comm_brick.cpp index 528456ab7c..74865476f5 100644 --- a/src/comm_brick.cpp +++ b/src/comm_brick.cpp @@ -252,9 +252,9 @@ void CommBrick::setup() jgroup = map_type_multi[j]; if(cutmultisq[jgroup][jgroup] > cutmultisq[igroup][igroup]) continue; - cutghostmulti[i][0] = MAX(cutghostmulti[i][0],sqrt(cutmultisq[igroup][jgroup])); - cutghostmulti[i][1] = MAX(cutghostmulti[i][1],sqrt(cutmultisq[igroup][jgroup])); - cutghostmulti[i][2] = MAX(cutghostmulti[i][2],sqrt(cutmultisq[igroup][jgroup])); + cutghostmulti[i][0] = length0 * MAX(cutghostmulti[i][0],sqrt(cutmultisq[igroup][jgroup])); + cutghostmulti[i][1] = length1 * MAX(cutghostmulti[i][1],sqrt(cutmultisq[igroup][jgroup])); + cutghostmulti[i][2] = length2 * MAX(cutghostmulti[i][2],sqrt(cutmultisq[igroup][jgroup])); } } } else { diff --git a/src/nbin.cpp b/src/nbin.cpp index 084d3bfc4a..cb10014b36 100644 --- a/src/nbin.cpp +++ b/src/nbin.cpp @@ -43,8 +43,11 @@ NBin::NBin(LAMMPS *lmp) : Pointers(lmp) bins_multi = nullptr; atom2bin_multi = nullptr; maxbins_multi = nullptr; + + map_type_multi = nullptr; + cutmultisq = nullptr; - maxtypes = 0; + maxgroups = 0; neighbor->last_setup_bins = -1; @@ -84,7 +87,7 @@ NBin::~NBin() memory->destroy(bininvy_multi); memory->destroy(bininvz_multi); - for (int n = 1; n <= maxtypes; n++) { + for (int n = 0; n < maxgroups; n++) { memory->destroy(binhead_multi[n]); memory->destroy(bins_multi[n]); memory->destroy(atom2bin_multi[n]); @@ -117,6 +120,10 @@ void NBin::copy_neighbor_info() binsize_user = neighbor->binsize_user; bboxlo = neighbor->bboxlo; bboxhi = neighbor->bboxhi; + + n_multi_groups = neighbor->n_multi_groups; + map_type_multi = neighbor->map_type_multi; + cutmultisq = neighbor->cutmultisq; // overwrite Neighbor cutoff with custom value set by requestor // only works for style = BIN (checked by Neighbor class) @@ -174,10 +181,10 @@ int NBin::coord2bin(double *x) /* ---------------------------------------------------------------------- - convert atom coords into local bin # for a particular type + convert atom coords into local bin # for a particular grouping ------------------------------------------------------------------------- */ -int NBin::coord2bin_multi(double *x, int it) +int NBin::coord2bin_multi(double *x, int ig) { int ix,iy,iz; int ibin; @@ -186,33 +193,33 @@ int NBin::coord2bin_multi(double *x, int it) error->one(FLERR,"Non-numeric positions - simulation unstable"); if (x[0] >= bboxhi[0]) - ix = static_cast ((x[0]-bboxhi[0])*bininvx_multi[it]) + nbinx_multi[it]; + ix = static_cast ((x[0]-bboxhi[0])*bininvx_multi[ig]) + nbinx_multi[ig]; else if (x[0] >= bboxlo[0]) { - ix = static_cast ((x[0]-bboxlo[0])*bininvx_multi[it]); - ix = MIN(ix,nbinx_multi[it]-1); + ix = static_cast ((x[0]-bboxlo[0])*bininvx_multi[ig]); + ix = MIN(ix,nbinx_multi[ig]-1); } else - ix = static_cast ((x[0]-bboxlo[0])*bininvx_multi[it]) - 1; + ix = static_cast ((x[0]-bboxlo[0])*bininvx_multi[ig]) - 1; if (x[1] >= bboxhi[1]) - iy = static_cast ((x[1]-bboxhi[1])*bininvy_multi[it]) + nbiny_multi[it]; + iy = static_cast ((x[1]-bboxhi[1])*bininvy_multi[ig]) + nbiny_multi[ig]; else if (x[1] >= bboxlo[1]) { - iy = static_cast ((x[1]-bboxlo[1])*bininvy_multi[it]); - iy = MIN(iy,nbiny_multi[it]-1); + iy = static_cast ((x[1]-bboxlo[1])*bininvy_multi[ig]); + iy = MIN(iy,nbiny_multi[ig]-1); } else - iy = static_cast ((x[1]-bboxlo[1])*bininvy_multi[it]) - 1; + iy = static_cast ((x[1]-bboxlo[1])*bininvy_multi[ig]) - 1; if (x[2] >= bboxhi[2]) - iz = static_cast ((x[2]-bboxhi[2])*bininvz_multi[it]) + nbinz_multi[it]; + iz = static_cast ((x[2]-bboxhi[2])*bininvz_multi[ig]) + nbinz_multi[ig]; else if (x[2] >= bboxlo[2]) { - iz = static_cast ((x[2]-bboxlo[2])*bininvz_multi[it]); - iz = MIN(iz,nbinz_multi[it]-1); + iz = static_cast ((x[2]-bboxlo[2])*bininvz_multi[ig]); + iz = MIN(iz,nbinz_multi[ig]-1); } else - iz = static_cast ((x[2]-bboxlo[2])*bininvz_multi[it]) - 1; + iz = static_cast ((x[2]-bboxlo[2])*bininvz_multi[ig]) - 1; - ibin = (iz-mbinzlo_multi[it])*mbiny_multi[it]*mbinx_multi[it] - + (iy-mbinylo_multi[it])*mbinx_multi[it] - + (ix-mbinxlo_multi[it]); + ibin = (iz-mbinzlo_multi[ig])*mbiny_multi[ig]*mbinx_multi[ig] + + (iy-mbinylo_multi[ig])*mbinx_multi[ig] + + (ix-mbinxlo_multi[ig]); return ibin; } diff --git a/src/nbin.h b/src/nbin.h index 8f95ad987e..19035d619e 100644 --- a/src/nbin.h +++ b/src/nbin.h @@ -75,6 +75,9 @@ class NBin : protected Pointers { int binsizeflag; double binsize_user; double *bboxlo,*bboxhi; + int n_multi_groups; + int *map_type_multi; + double **cutmultisq; // data common to all NBin variants @@ -84,15 +87,12 @@ class NBin : protected Pointers { // data for standard NBin int maxatom; // size of bins array - - // data for standard NBin - int maxbin; // size of binhead array - // data for multi/multi NBin + // data for multi NBin - int maxtypes; // size of multi arrays - int * maxbins_multi; // size of 2nd dimension of binhead_multi array + int maxgroups; // size of multi arrays + int * maxbins_multi; // size of 2nd dimension of binhead_multi array // methods diff --git a/src/nbin_multi.cpp b/src/nbin_multi.cpp index 302fe30d88..b39489c3fe 100644 --- a/src/nbin_multi.cpp +++ b/src/nbin_multi.cpp @@ -38,7 +38,7 @@ void NBinMulti::bin_atoms_setup(int nall) { // binhead_multi[n] = per-bin vector mbins in length mbins_multi[n] - for (int n = 1; n <= maxtypes; n++) { + for (int n = 0; n < maxgroups; n++) { if (mbins_multi[n] > maxbins_multi[n]) { maxbins_multi[n] = mbins_multi[n]; memory->destroy(binhead_multi[n]); @@ -48,10 +48,12 @@ void NBinMulti::bin_atoms_setup(int nall) // bins_multi[n] and atom2bin_multi[n] = per-atom vectors // for both local and ghost atoms + + // TODO: maxatom should be maxatom per group if (nall > maxatom) { maxatom = nall; - for (int n = 1; n <= maxtypes; n++) { + for (int n = 0; n < maxgroups; n++) { memory->destroy(bins_multi[n]); memory->destroy(atom2bin_multi[n]); memory->create(bins_multi[n], maxatom, "neigh:bin_multi"); @@ -61,23 +63,16 @@ void NBinMulti::bin_atoms_setup(int nall) } /* --------------------------------------------------------------------- - Identify index of type with smallest cutoff + Identify index of group with smallest cutoff ------------------------------------------------------------------------ */ -int NBinMulti::itype_min() { +int NBinMulti::igroup_min() +{ + int imin = 0; + for (int n = 0; n < maxgroups; n++) + if (cutmultisq[n][n] < cutmultisq[imin][imin]) imin = n; - int itypemin = 1; - double ** cutneighsq; - - cutneighsq = neighbor->cutneighsq; - - for (int n = 1; n <= atom->ntypes; n++) { - if (cutneighsq[n][n] < cutneighsq[itypemin][itypemin]) { - itypemin = n; - } - } - - return itypemin; + return imin; } /* ---------------------------------------------------------------------- @@ -87,15 +82,15 @@ int NBinMulti::itype_min() { void NBinMulti::setup_bins(int style) { int n; - int itypemin; + int igroupmin; // Initialize arrays - if (atom->ntypes > maxtypes) { + if (n_multi_groups > maxgroups) { // Clear any/all memory for existing types - for (n = 1; n <= maxtypes; n++) { + for (n = 0; n < maxgroups; n++) { memory->destroy(atom2bin_multi[n]); memory->destroy(binhead_multi[n]); memory->destroy(bins_multi[n]); @@ -106,59 +101,60 @@ void NBinMulti::setup_bins(int style) // Realloacte at updated maxtypes - maxtypes = atom->ntypes; + maxgroups = n_multi_groups; - atom2bin_multi = new int*[maxtypes+1](); - binhead_multi = new int*[maxtypes+1](); - bins_multi = new int*[maxtypes+1](); + atom2bin_multi = new int*[maxgroups](); + binhead_multi = new int*[maxgroups](); + bins_multi = new int*[maxgroups](); memory->destroy(nbinx_multi); memory->destroy(nbiny_multi); memory->destroy(nbinz_multi); - memory->create(nbinx_multi, maxtypes+1, "neigh:nbinx_multi"); - memory->create(nbiny_multi, maxtypes+1, "neigh:nbiny_multi"); - memory->create(nbinz_multi, maxtypes+1, "neigh:nbinz_multi"); + memory->create(nbinx_multi, maxgroups, "neigh:nbinx_multi"); + memory->create(nbiny_multi, maxgroups, "neigh:nbiny_multi"); + memory->create(nbinz_multi, maxgroups, "neigh:nbinz_multi"); memory->destroy(mbins_multi); memory->destroy(mbinx_multi); memory->destroy(mbiny_multi); memory->destroy(mbinz_multi); - memory->create(mbins_multi, maxtypes+1, "neigh:mbins_multi"); - memory->create(mbinx_multi, maxtypes+1, "neigh:mbinx_multi"); - memory->create(mbiny_multi, maxtypes+1, "neigh:mbiny_multi"); - memory->create(mbinz_multi, maxtypes+1, "neigh:mbinz_multi"); + memory->create(mbins_multi, maxgroups, "neigh:mbins_multi"); + memory->create(mbinx_multi, maxgroups, "neigh:mbinx_multi"); + memory->create(mbiny_multi, maxgroups, "neigh:mbiny_multi"); + memory->create(mbinz_multi, maxgroups, "neigh:mbinz_multi"); memory->destroy(mbinxlo_multi); memory->destroy(mbinylo_multi); memory->destroy(mbinzlo_multi); - memory->create(mbinxlo_multi, maxtypes+1, "neigh:mbinxlo_multi"); - memory->create(mbinylo_multi, maxtypes+1, "neigh:mbinylo_multi"); - memory->create(mbinzlo_multi, maxtypes+1, "neigh:mbinzlo_multi"); + memory->create(mbinxlo_multi, maxgroups, "neigh:mbinxlo_multi"); + memory->create(mbinylo_multi, maxgroups, "neigh:mbinylo_multi"); + memory->create(mbinzlo_multi, maxgroups, "neigh:mbinzlo_multi"); memory->destroy(binsizex_multi); memory->destroy(binsizey_multi); memory->destroy(binsizez_multi); - memory->create(binsizex_multi, maxtypes+1, "neigh:binsizex_multi"); - memory->create(binsizey_multi, maxtypes+1, "neigh:binsizey_multi"); - memory->create(binsizez_multi, maxtypes+1, "neigh:binsizez_multi"); + memory->create(binsizex_multi, maxgroups, "neigh:binsizex_multi"); + memory->create(binsizey_multi, maxgroups, "neigh:binsizey_multi"); + memory->create(binsizez_multi, maxgroups, "neigh:binsizez_multi"); memory->destroy(bininvx_multi); memory->destroy(bininvy_multi); memory->destroy(bininvz_multi); - memory->create(bininvx_multi, maxtypes+1, "neigh:bininvx_multi"); - memory->create(bininvy_multi, maxtypes+1, "neigh:bininvy_multi"); - memory->create(bininvz_multi, maxtypes+1, "neigh:bininvz_multi"); + memory->create(bininvx_multi, maxgroups, "neigh:bininvx_multi"); + memory->create(bininvy_multi, maxgroups, "neigh:bininvy_multi"); + memory->create(bininvz_multi, maxgroups, "neigh:bininvz_multi"); memory->destroy(maxbins_multi); - memory->create(maxbins_multi, maxtypes+1, "neigh:maxbins_multi"); + memory->create(maxbins_multi, maxgroups, "neigh:maxbins_multi"); + // make sure reallocation occurs in bin_atoms_setup() - for (n = 1; n <= maxtypes; n++) { + for (n = 0; n < maxgroups; n++) { maxbins_multi[n] = 0; } maxatom = 0; } - itypemin = itype_min(); + igroupmin = igroup_min(); // bbox = size of bbox of entire domain // bsubbox lo/hi = bounding box of my subdomain extended by comm->cutghost @@ -193,16 +189,16 @@ void NBinMulti::setup_bins(int style) bbox[1] = bboxhi[1] - bboxlo[1]; bbox[2] = bboxhi[2] - bboxlo[2]; - // For each type... + // For each grouping... - for (n = 1; n <= atom->ntypes; n++) { + for (n = 0; n < maxgroups; n++) { // binsize_user only relates to smallest type // optimal bin size is roughly 1/2 the type-type cutoff // special case of all cutoffs = 0.0, binsize = box size double binsize_optimal; - if (n == itypemin && binsizeflag) binsize_optimal = binsize_user; - else binsize_optimal = 0.5*sqrt(neighbor->cutneighsq[n][n]); + if (n == igroupmin && binsizeflag) binsize_optimal = binsize_user; + else binsize_optimal = 0.5*sqrt(cutmultisq[n][n]); if (binsize_optimal == 0.0) binsize_optimal = bbox[0]; double binsizeinv = 1.0/binsize_optimal; @@ -306,7 +302,7 @@ void NBinMulti::bin_atoms() int i,ibin,n; last_bin = update->ntimestep; - for (n = 1; n <= maxtypes; n++) { + for (n = 0; n < maxgroups; n++) { for (i = 0; i < mbins_multi[n]; i++) binhead_multi[n][i] = -1; } @@ -323,7 +319,7 @@ void NBinMulti::bin_atoms() int bitmask = group->bitmask[includegroup]; for (i = nall-1; i >= nlocal; i--) { if (mask[i] & bitmask) { - n = type[i]; + n = map_type_multi[type[i]]; ibin = coord2bin_multi(x[i], n); atom2bin_multi[n][i] = ibin; bins_multi[n][i] = binhead_multi[n][ibin]; @@ -331,7 +327,7 @@ void NBinMulti::bin_atoms() } } for (i = atom->nfirst-1; i >= 0; i--) { - n = type[i]; + n = map_type_multi[type[i]]; ibin = coord2bin_multi(x[i], n); atom2bin_multi[n][i] = ibin; bins_multi[n][i] = binhead_multi[n][ibin]; @@ -339,7 +335,7 @@ void NBinMulti::bin_atoms() } } else { for (i = nall-1; i >= 0; i--) { - n = type[i]; + n = map_type_multi[type[i]]; ibin = coord2bin_multi(x[i], n); atom2bin_multi[n][i] = ibin; bins_multi[n][i] = binhead_multi[n][ibin]; @@ -354,7 +350,7 @@ double NBinMulti::memory_usage() { double bytes = 0; - for (int m = 1; m < maxtypes; m++) { + for (int m = 0; m < maxgroups; m++) { bytes += maxbins_multi[m]*sizeof(int); bytes += 2*maxatom*sizeof(int); } diff --git a/src/nbin_multi.h b/src/nbin_multi.h index 856f3fa722..84ec0348bc 100644 --- a/src/nbin_multi.h +++ b/src/nbin_multi.h @@ -38,7 +38,7 @@ class NBinMulti : public NBin { private: - int itype_min(); + int igroup_min(); }; } diff --git a/src/neighbor.cpp b/src/neighbor.cpp index b08de6cb02..94add1f324 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -350,8 +350,8 @@ void Neighbor::init() int igroup, jgroup; // If not defined, create default mapping if(not map_type_multi) { - memory->create(map_type_multi,n+1,"neigh:map_type_multi"); n_multi_groups = n; + memory->create(map_type_multi,n+1,"neigh:map_type_multi"); for(i = 1; i <= n; i++) map_type_multi[i] = i-1; } diff --git a/src/npair.cpp b/src/npair.cpp index 07e476e1dd..b2ecc953ea 100644 --- a/src/npair.cpp +++ b/src/npair.cpp @@ -95,6 +95,12 @@ void NPair::copy_neighbor_info() special_flag = neighbor->special_flag; + // multi info + + n_multi_groups = neighbor->n_multi_groups; + map_type_multi = neighbor->map_type_multi; + cutmultisq = neighbor->cutmultisq; + // overwrite per-type Neighbor cutoffs with custom value set by requestor // only works for style = BIN (checked by Neighbor class) @@ -265,10 +271,10 @@ int NPair::coord2bin(double *x, int &ix, int &iy, int &iz) /* ---------------------------------------------------------------------- - same as coord2bin in NbinMulti2 + multi version of coord2bin for a given grouping ------------------------------------------------------------------------- */ -int NPair::coord2bin(double *x, int it) +int NPair::coord2bin(double *x, int ig) { int ix,iy,iz; int ibin; @@ -277,32 +283,32 @@ int NPair::coord2bin(double *x, int it) error->one(FLERR,"Non-numeric positions - simulation unstable"); if (x[0] >= bboxhi[0]) - ix = static_cast ((x[0]-bboxhi[0])*bininvx_multi[it]) + nbinx_multi[it]; + ix = static_cast ((x[0]-bboxhi[0])*bininvx_multi[ig]) + nbinx_multi[ig]; else if (x[0] >= bboxlo[0]) { - ix = static_cast ((x[0]-bboxlo[0])*bininvx_multi[it]); - ix = MIN(ix,nbinx_multi[it]-1); + ix = static_cast ((x[0]-bboxlo[0])*bininvx_multi[ig]); + ix = MIN(ix,nbinx_multi[ig]-1); } else - ix = static_cast ((x[0]-bboxlo[0])*bininvx_multi[it]) - 1; + ix = static_cast ((x[0]-bboxlo[0])*bininvx_multi[ig]) - 1; if (x[1] >= bboxhi[1]) - iy = static_cast ((x[1]-bboxhi[1])*bininvy_multi[it]) + nbiny_multi[it]; + iy = static_cast ((x[1]-bboxhi[1])*bininvy_multi[ig]) + nbiny_multi[ig]; else if (x[1] >= bboxlo[1]) { - iy = static_cast ((x[1]-bboxlo[1])*bininvy_multi[it]); - iy = MIN(iy,nbiny_multi[it]-1); + iy = static_cast ((x[1]-bboxlo[1])*bininvy_multi[ig]); + iy = MIN(iy,nbiny_multi[ig]-1); } else - iy = static_cast ((x[1]-bboxlo[1])*bininvy_multi[it]) - 1; + iy = static_cast ((x[1]-bboxlo[1])*bininvy_multi[ig]) - 1; if (x[2] >= bboxhi[2]) - iz = static_cast ((x[2]-bboxhi[2])*bininvz_multi[it]) + nbinz_multi[it]; + iz = static_cast ((x[2]-bboxhi[2])*bininvz_multi[ig]) + nbinz_multi[ig]; else if (x[2] >= bboxlo[2]) { - iz = static_cast ((x[2]-bboxlo[2])*bininvz_multi[it]); - iz = MIN(iz,nbinz_multi[it]-1); + iz = static_cast ((x[2]-bboxlo[2])*bininvz_multi[ig]); + iz = MIN(iz,nbinz_multi[ig]-1); } else - iz = static_cast ((x[2]-bboxlo[2])*bininvz_multi[it]) - 1; + iz = static_cast ((x[2]-bboxlo[2])*bininvz_multi[ig]) - 1; - ibin = (iz-mbinzlo_multi[it])*mbiny_multi[it]*mbinx_multi[it] - + (iy-mbinylo_multi[it])*mbinx_multi[it] - + (ix-mbinxlo_multi[it]); + ibin = (iz-mbinzlo_multi[ig])*mbiny_multi[ig]*mbinx_multi[ig] + + (iy-mbinylo_multi[ig])*mbinx_multi[ig] + + (ix-mbinxlo_multi[ig]); return ibin; } \ No newline at end of file diff --git a/src/npair.h b/src/npair.h index c56c6bdb20..4d99e3a8e6 100644 --- a/src/npair.h +++ b/src/npair.h @@ -48,7 +48,10 @@ class NPair : protected Pointers { double cut_middle_sq; double cut_middle_inside_sq; double *bboxlo,*bboxhi; - + int n_multi_groups; + int *map_type_multi; + double **cutmultisq; + // exclusion data from Neighbor class int nex_type; // # of entries in type exclusion list diff --git a/src/npair_full_multi.cpp b/src/npair_full_multi.cpp index 4463964d5f..f29f04107a 100644 --- a/src/npair_full_multi.cpp +++ b/src/npair_full_multi.cpp @@ -28,13 +28,13 @@ NPairFullMulti::NPairFullMulti(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- binned neighbor list construction for all neighbors - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent every neighbor pair appears in list of both atoms i and j ------------------------------------------------------------------------- */ void NPairFullMulti::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,jbin,which,ns,imol,iatom,moltemplate; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,which,ns,imol,iatom,moltemplate; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; @@ -67,8 +67,8 @@ void NPairFullMulti::build(NeighList *list) for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -78,29 +78,30 @@ void NPairFullMulti::build(NeighList *list) tagprev = tag[i] - iatom - 1; } - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); // loop over all atoms in surrounding bins in stencil including self // skip i = j - // use full stencil for all type combinations + // use full stencil for all group combinations - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { if (i == j) continue; - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; diff --git a/src/npair_half_multi_newtoff.cpp b/src/npair_half_multi_newtoff.cpp index c0dcaf0c98..85832b11b1 100755 --- a/src/npair_half_multi_newtoff.cpp +++ b/src/npair_half_multi_newtoff.cpp @@ -28,7 +28,7 @@ NPairHalfMultiNewtoff::NPairHalfMultiNewtoff(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- binned neighbor list construction with partial Newton's 3rd law - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks own bin and other bins in stencil pair stored once if i,j are both owned and i < j pair stored by me if j is ghost (also stored by proc owning j) @@ -36,7 +36,7 @@ NPairHalfMultiNewtoff::NPairHalfMultiNewtoff(LAMMPS *lmp) : NPair(lmp) {} void NPairHalfMultiNewtoff::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,jbin,which,ns,imol,iatom,moltemplate; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,which,ns,imol,iatom,moltemplate; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; @@ -69,8 +69,8 @@ void NPairHalfMultiNewtoff::build(NeighList *list) for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -80,30 +80,31 @@ void NPairHalfMultiNewtoff::build(NeighList *list) tagprev = tag[i] - iatom - 1; } - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); // loop over all atoms in other bins in stencil including self // only store pair if i < j // stores own/own pairs only once // stores own/ghost pairs on both procs - // use full stencil for all type combinations + // use full stencil for all group combinations - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >=0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >=0; j = bins_multi[jgroup][j]) { if (j <= i) continue; - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/npair_half_multi_newton.cpp b/src/npair_half_multi_newton.cpp index e097052fe5..d18b009460 100755 --- a/src/npair_half_multi_newton.cpp +++ b/src/npair_half_multi_newton.cpp @@ -28,14 +28,14 @@ NPairHalfMultiNewton::NPairHalfMultiNewton(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- binned neighbor list construction with full Newton's 3rd law - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks its own bin and other bins in Newton stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void NPairHalfMultiNewton::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,jbin,which,ns,imol,iatom,moltemplate; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,which,ns,imol,iatom,moltemplate; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; @@ -68,8 +68,8 @@ void NPairHalfMultiNewton::build(NeighList *list) for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -79,28 +79,28 @@ void NPairHalfMultiNewton::build(NeighList *list) tagprev = tag[i] - iatom - 1; } - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); - if(cutneighsq[itype][itype] == cutneighsq[jtype][jtype]){ + if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){ // if same size: use half stencil - if(itype == jtype){ + if(igroup == jgroup){ - // if same type, implement with: + // if same group, implement with: // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i - js = bins_multi[itype][i]; + js = bins_multi[igroup][i]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { if (j >= nlocal) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { @@ -109,7 +109,8 @@ void NPairHalfMultiNewton::build(NeighList *list) } } - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -134,14 +135,14 @@ void NPairHalfMultiNewton::build(NeighList *list) } } else { - // if different types, implement with: - // loop over all atoms in jtype bin + // if different groups, implement with: + // loop over all atoms in jgroup bin // if j is owned atom, store it if j > i // if j is ghost, only store if j coords are "above and to the right" of i - js = binhead_multi[jtype][jbin]; + js = binhead_multi[jgroup][jbin]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { if(j < i) continue; if (j >= nlocal) { @@ -152,7 +153,8 @@ void NPairHalfMultiNewton::build(NeighList *list) } } - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -178,20 +180,21 @@ void NPairHalfMultiNewton::build(NeighList *list) } } - // for all types, loop over all atoms in other bins in stencil, store every pair + // for all groups, loop over all atoms in other bins in stencil, store every pair // stencil is empty if i larger than j // stencil is half if i same size as j // stencil is full if i smaller than j - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { - - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; diff --git a/src/npair_half_multi_newton_tri.cpp b/src/npair_half_multi_newton_tri.cpp index fabaf5f202..ebf871ad99 100755 --- a/src/npair_half_multi_newton_tri.cpp +++ b/src/npair_half_multi_newton_tri.cpp @@ -28,14 +28,14 @@ NPairHalfMultiNewtonTri::NPairHalfMultiNewtonTri(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- binned neighbor list construction with Newton's 3rd law for triclinic - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks its own bin and other bins in triclinic stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void NPairHalfMultiNewtonTri::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,jbin,which,ns,imol,iatom,moltemplate; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,which,ns,imol,iatom,moltemplate; tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; @@ -68,8 +68,8 @@ void NPairHalfMultiNewtonTri::build(NeighList *list) for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; @@ -79,14 +79,14 @@ void NPairHalfMultiNewtonTri::build(NeighList *list) tagprev = tag[i] - iatom - 1; } - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); // loop over all atoms in bins in stencil // stencil is empty if i larger than j @@ -97,15 +97,15 @@ void NPairHalfMultiNewtonTri::build(NeighList *list) // (equal zyx and j <= i) // latter excludes self-self interaction but allows superposed atoms - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { - // if same size (e.g. same type), use half stencil - if(cutneighsq[itype][itype] == cutneighsq[jtype][jtype]){ + // if same size (same group), use half stencil + if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){ if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; @@ -116,8 +116,9 @@ void NPairHalfMultiNewtonTri::build(NeighList *list) } } - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; diff --git a/src/npair_half_size_multi_newtoff.cpp b/src/npair_half_size_multi_newtoff.cpp index 65b6b68fb2..ec4fbad9aa 100644 --- a/src/npair_half_size_multi_newtoff.cpp +++ b/src/npair_half_size_multi_newtoff.cpp @@ -28,7 +28,7 @@ NPairHalfSizeMultiNewtoff::NPairHalfSizeMultiNewtoff(LAMMPS *lmp) : NPair(lmp) { /* ---------------------------------------------------------------------- size particles binned neighbor list construction with partial Newton's 3rd law - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks own bin and other bins in stencil pair stored once if i,j are both owned and i < j pair stored by me if j is ghost (also stored by proc owning j) @@ -36,7 +36,7 @@ NPairHalfSizeMultiNewtoff::NPairHalfSizeMultiNewtoff(LAMMPS *lmp) : NPair(lmp) { void NPairHalfSizeMultiNewtoff::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,jbin,ns; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,ns; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -64,37 +64,38 @@ void NPairHalfSizeMultiNewtoff::build(NeighList *list) for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); // loop over all atoms in other bins in stencil including self // only store pair if i < j // stores own/own pairs only once // stores own/ghost pairs on both procs - // use full stencil for all type combinations + // use full stencil for all group combinations - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >=0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >=0; j = bins_multi[jgroup][j]) { if (j <= i) continue; - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/npair_half_size_multi_newton.cpp b/src/npair_half_size_multi_newton.cpp index 8fa3b9bc1b..c3a13f4b81 100755 --- a/src/npair_half_size_multi_newton.cpp +++ b/src/npair_half_size_multi_newton.cpp @@ -27,14 +27,14 @@ NPairHalfSizeMultiNewton::NPairHalfSizeMultiNewton(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- size particles binned neighbor list construction with full Newton's 3rd law - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks its own bin and other bins in Newton stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void NPairHalfSizeMultiNewton::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,jbin,ns,js; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,ns,js; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -61,35 +61,35 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); - if(cutneighsq[itype][itype] == cutneighsq[jtype][jtype]){ + if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){ // if same size: use half stencil - if(itype == jtype){ + if(igroup == jgroup){ - // if same type, implement with: + // if same group, implement with: // loop over rest of atoms in i's bin, ghosts are at end of linked list // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above and to the right" of i - js = bins_multi[itype][i]; + js = bins_multi[igroup][i]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { if (j >= nlocal) { if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { @@ -98,7 +98,8 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) } } - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -116,14 +117,14 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) } } else { - // if different types, implement with: - // loop over all atoms in jtype bin + // if different groups, implement with: + // loop over all atoms in jgroup bin // if j is owned atom, store it if j > i // if j is ghost, only store if j coords are "above and to the right" of i - js = binhead_multi[jtype][jbin]; + js = binhead_multi[jgroup][jbin]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { if(j < i) continue; if (j >= nlocal) { @@ -134,8 +135,9 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) } } - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; @@ -153,20 +155,21 @@ void NPairHalfSizeMultiNewton::build(NeighList *list) } } - // for all types, loop over all atoms in other bins in stencil, store every pair + // for all groups, loop over all atoms in other bins in stencil, store every pair // stencil is empty if i larger than j // stencil is half if i same size as j // stencil is full if i smaller than j - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { - if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; - + jtype = type[j]; + if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; diff --git a/src/npair_half_size_multi_newton_tri.cpp b/src/npair_half_size_multi_newton_tri.cpp index 780fa95f0b..7c792c6a19 100755 --- a/src/npair_half_size_multi_newton_tri.cpp +++ b/src/npair_half_size_multi_newton_tri.cpp @@ -27,14 +27,14 @@ NPairHalfSizeMultiNewtonTri::NPairHalfSizeMultiNewtonTri(LAMMPS *lmp) : NPair(lm /* ---------------------------------------------------------------------- size particles binned neighbor list construction with Newton's 3rd law for triclinic - multi-type stencil is itype-jtype dependent + multi stencil is igroup-jgroup dependent each owned atom i checks its own bin and other bins in triclinic stencil every pair stored exactly once by some processor ------------------------------------------------------------------------- */ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) { - int i,j,k,n,itype,jtype,ibin,jbin,ns,js; + int i,j,k,n,itype,jtype,igroup,jgroup,ibin,jbin,ns,js; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; double radi,radsum,cutdistsq; int *neighptr,*s; @@ -61,21 +61,21 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) for (i = 0; i < nlocal; i++) { n = 0; neighptr = ipage->vget(); - itype = type[i]; + igroup = map_type_multi[itype]; xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; radi = radius[i]; - ibin = atom2bin_multi[itype][i]; + ibin = atom2bin_multi[igroup][i]; - // loop through stencils for all types - for (jtype = 1; jtype <= atom->ntypes; jtype++) { + // loop through stencils for all groups + for (jgroup = 0; jgroup < n_multi_groups; jgroup++) { - // if same type use own bin - if(itype == jtype) jbin = ibin; - else jbin = coord2bin(x[i], jtype); + // if same group use own bin + if(igroup == jgroup) jbin = ibin; + else jbin = coord2bin(x[i], jgroup); // loop over all atoms in bins in stencil @@ -87,15 +87,15 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) // (equal zyx and j <= i) // latter excludes self-self interaction but allows superposed atoms - s = stencil_multi[itype][jtype]; - ns = nstencil_multi[itype][jtype]; + s = stencil_multi[igroup][jgroup]; + ns = nstencil_multi[igroup][jgroup]; for (k = 0; k < ns; k++) { - js = binhead_multi[jtype][jbin + s[k]]; - for (j = js; j >= 0; j = bins_multi[jtype][j]) { + js = binhead_multi[jgroup][jbin + s[k]]; + for (j = js; j >= 0; j = bins_multi[jgroup][j]) { - // if same size (e.g. same type), use half stencil - if(cutneighsq[itype][itype] == cutneighsq[jtype][jtype]){ + // if same size (same group), use half stencil + if(cutmultisq[igroup][igroup] == cutmultisq[jgroup][jgroup]){ if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; @@ -106,6 +106,7 @@ void NPairHalfSizeMultiNewtonTri::build(NeighList *list) } } + jtype = type[j]; if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue; delx = xtmp - x[j][0]; diff --git a/src/nstencil.cpp b/src/nstencil.cpp index f6281a1964..acd7272245 100644 --- a/src/nstencil.cpp +++ b/src/nstencil.cpp @@ -52,14 +52,14 @@ using namespace LAMMPS_NS; cutoff is not cutneighmaxsq, but max cutoff for that atom type no versions that allow ghost on (any need for it?) for multi: - create one stencil for each itype-jtype pairing - the full/half stencil label refers to the same-type stencil - a half list with newton on has a half same-type stencil - a full list or half list with newton of has a full same-type stencil - cross type stencils are always full to allow small-to-large lookups + create one stencil for each igroup-jgroup pairing + the full/half stencil label refers to the same-group stencil + a half list with newton on has a half same-group stencil + a full list or half list with newton of has a full same-group stencil + cross group stencils are always full to allow small-to-large lookups for orthogonal boxes, a half stencil includes bins to the "upper right" of central bin for triclinic, a half stencil includes bins in the z (3D) or y (2D) plane of self and above - cutoff is not cutneighmaxsq, but max cutoff for that atom type + cutoff is not cutneighmaxsq, but max cutoff for that atom group no versions that allow ghost on (any need for it?) ------------------------------------------------------------------------- */ @@ -81,7 +81,7 @@ NStencil::NStencil(LAMMPS *lmp) : Pointers(lmp) flag_half_multi = nullptr; flag_skip_multi = nullptr; - bin_type_multi = nullptr; + bin_group_multi = nullptr; dimension = domain->dimension; } @@ -107,10 +107,10 @@ NStencil::~NStencil() if (stencil_multi) { - int n = atom->ntypes; + int n = n_multi_groups; memory->destroy(nstencil_multi); - for (int i = 1; i <= n; i++) { - for (int j = 0; j <= n; j++) { + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { if (! flag_skip_multi[i][j]) memory->destroy(stencil_multi[i][j]); } @@ -120,7 +120,7 @@ NStencil::~NStencil() memory->destroy(maxstencil_multi); memory->destroy(flag_half_multi); memory->destroy(flag_skip_multi); - memory->destroy(bin_type_multi); + memory->destroy(bin_group_multi); memory->destroy(stencil_sx_multi); memory->destroy(stencil_sy_multi); @@ -155,6 +155,10 @@ void NStencil::copy_neighbor_info() cutneighmaxsq = neighbor->cutneighmaxsq; cuttypesq = neighbor->cuttypesq; cutneighsq = neighbor->cutneighsq; + + n_multi_groups = neighbor->n_multi_groups; + map_type_multi = neighbor->map_type_multi; + cutmultisq = neighbor->cutmultisq; // overwrite Neighbor cutoff with custom value set by requestor // only works for style = BIN (checked by Neighbor class) @@ -265,44 +269,44 @@ void NStencil::create_setup() } } } else { - int i, j, bin_type, smax; + int i, j, bin_group, smax; double stencil_range; - int n = atom->ntypes; + int n = n_multi_groups; if(nb) copy_bin_info_multi(); // Allocate arrays to store stencil information - memory->create(flag_half_multi, n+1, n+1, + memory->create(flag_half_multi, n, n, "neighstencil:flag_half_multi"); - memory->create(flag_skip_multi, n+1, n+1, + memory->create(flag_skip_multi, n, n, "neighstencil:flag_skip_multi"); - memory->create(bin_type_multi, n+1, n+1, - "neighstencil:bin_type_multi"); + memory->create(bin_group_multi, n, n, + "neighstencil:bin_group_multi"); - memory->create(stencil_sx_multi, n+1, n+1, + memory->create(stencil_sx_multi, n, n, "neighstencil:stencil_sx_multi"); - memory->create(stencil_sy_multi, n+1, n+1, + memory->create(stencil_sy_multi, n, n, "neighstencil:stencil_sy_multi"); - memory->create(stencil_sz_multi, n+1, n+1, + memory->create(stencil_sz_multi, n, n, "neighstencil:stencil_sz_multi"); - memory->create(stencil_binsizex_multi, n+1, n+1, + memory->create(stencil_binsizex_multi, n, n, "neighstencil:stencil_binsizex_multi"); - memory->create(stencil_binsizey_multi, n+1, n+1, + memory->create(stencil_binsizey_multi, n, n, "neighstencil:stencil_binsizey_multi"); - memory->create(stencil_binsizez_multi, n+1, n+1, + memory->create(stencil_binsizez_multi, n, n, "neighstencil:stencil_binsizez_multi"); - memory->create(stencil_mbinx_multi, n+1, n+1, + memory->create(stencil_mbinx_multi, n, n, "neighstencil:stencil_mbinx_multi"); - memory->create(stencil_mbiny_multi, n+1, n+1, + memory->create(stencil_mbiny_multi, n, n, "neighstencil:stencil_mbiny_multi"); - memory->create(stencil_mbinz_multi, n+1, n+1, + memory->create(stencil_mbinz_multi, n, n, "neighstencil:stencil_mbinz_multi"); // Skip all stencils by default, initialize smax - for (i = 1; i <= n; i++) { - for (j = 1; j <= n; j++) { + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { flag_skip_multi[i][j] = 1; } } @@ -313,43 +317,43 @@ void NStencil::create_setup() // Allocate arrays to store stencils if (!maxstencil_multi) { - memory->create(maxstencil_multi, n+1, n+1, "neighstencil::stencil_multi"); - memory->create(nstencil_multi, n+1, n+1, "neighstencil::nstencil_multi"); - stencil_multi = new int**[n+1](); - for (i = 1; i <= n; ++i) { - stencil_multi[i] = new int*[n+1](); - for (j = 1; j <= n; ++j) { + memory->create(maxstencil_multi, n, n, "neighstencil::stencil_multi"); + memory->create(nstencil_multi, n, n, "neighstencil::nstencil_multi"); + stencil_multi = new int**[n](); + for (i = 0; i < n; ++i) { + stencil_multi[i] = new int*[n](); + for (j = 0; j < n; ++j) { maxstencil_multi[i][j] = 0; nstencil_multi[i][j] = 0; } } } - for (i = 1; i <= n; i++) { - for (j = 1; j <= n; j++) { + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { // Skip creation of unused stencils if (flag_skip_multi[i][j]) continue; // Copy bin info for this pair of atom types - bin_type = bin_type_multi[i][j]; + bin_group = bin_group_multi[i][j]; - stencil_binsizex_multi[i][j] = binsizex_multi[bin_type]; - stencil_binsizey_multi[i][j] = binsizey_multi[bin_type]; - stencil_binsizez_multi[i][j] = binsizez_multi[bin_type]; + stencil_binsizex_multi[i][j] = binsizex_multi[bin_group]; + stencil_binsizey_multi[i][j] = binsizey_multi[bin_group]; + stencil_binsizez_multi[i][j] = binsizez_multi[bin_group]; - stencil_mbinx_multi[i][j] = mbinx_multi[bin_type]; - stencil_mbiny_multi[i][j] = mbiny_multi[bin_type]; - stencil_mbinz_multi[i][j] = mbinz_multi[bin_type]; + stencil_mbinx_multi[i][j] = mbinx_multi[bin_group]; + stencil_mbiny_multi[i][j] = mbiny_multi[bin_group]; + stencil_mbinz_multi[i][j] = mbinz_multi[bin_group]; - stencil_range = sqrt(cutneighsq[i][j]); + stencil_range = sqrt(cutmultisq[i][j]); - sx = static_cast (stencil_range*bininvx_multi[bin_type]); - if (sx*binsizex_multi[bin_type] < stencil_range) sx++; - sy = static_cast (stencil_range*bininvy_multi[bin_type]); - if (sy*binsizey_multi[bin_type] < stencil_range) sy++; - sz = static_cast (stencil_range*bininvz_multi[bin_type]); - if (sz*binsizez_multi[bin_type] < stencil_range) sz++; + sx = static_cast (stencil_range*bininvx_multi[bin_group]); + if (sx*binsizex_multi[bin_group] < stencil_range) sx++; + sy = static_cast (stencil_range*bininvy_multi[bin_group]); + if (sy*binsizey_multi[bin_group] < stencil_range) sy++; + sz = static_cast (stencil_range*bininvz_multi[bin_group]); + if (sz*binsizez_multi[bin_group] < stencil_range) sz++; stencil_sx_multi[i][j] = sx; stencil_sy_multi[i][j] = sy; @@ -393,24 +397,24 @@ double NStencil::bin_distance(int i, int j, int k) /* ---------------------------------------------------------------------- - compute closest distance for a given atom type + compute closest distance for a given atom grouping ------------------------------------------------------------------------- */ -double NStencil::bin_distance_multi(int i, int j, int k, int type) +double NStencil::bin_distance_multi(int i, int j, int k, int group) { double delx,dely,delz; - if (i > 0) delx = (i-1)*binsizex_multi[type]; + if (i > 0) delx = (i-1)*binsizex_multi[group]; else if (i == 0) delx = 0.0; - else delx = (i+1)*binsizex_multi[type]; + else delx = (i+1)*binsizex_multi[group]; - if (j > 0) dely = (j-1)*binsizey_multi[type]; + if (j > 0) dely = (j-1)*binsizey_multi[group]; else if (j == 0) dely = 0.0; - else dely = (j+1)*binsizey_multi[type]; + else dely = (j+1)*binsizey_multi[group]; - if (k > 0) delz = (k-1)*binsizez_multi[type]; + if (k > 0) delz = (k-1)*binsizez_multi[group]; else if (k == 0) delz = 0.0; - else delz = (k+1)*binsizez_multi[type]; + else delz = (k+1)*binsizez_multi[group]; return (delx*delx + dely*dely + delz*delz); } @@ -427,9 +431,9 @@ double NStencil::memory_usage() bytes += atom->ntypes*maxstencil_multi_old * sizeof(int); bytes += atom->ntypes*maxstencil_multi_old * sizeof(double); } else if (neighstyle == Neighbor::MULTI) { - int n = atom->ntypes; - for (int i = 1; i <= n; i++) { - for (int j = 1; j <= n; j++) { + int n = n_multi_groups; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { bytes += maxstencil_multi[i][j] * sizeof(int); bytes += maxstencil_multi[i][j] * sizeof(int); bytes += maxstencil_multi[i][j] * sizeof(double); diff --git a/src/nstencil.h b/src/nstencil.h index 9ca8bab55d..b06f340d6f 100644 --- a/src/nstencil.h +++ b/src/nstencil.h @@ -44,7 +44,7 @@ class NStencil : protected Pointers { // Arrays to store options for multi itype-jtype stencils bool **flag_half_multi; // flag creation of a half stencil for itype-jtype bool **flag_skip_multi; // skip creation of itype-jtype stencils (for newton on) - int **bin_type_multi; // what type to use for bin information + int **bin_group_multi; // what group to use for bin information NStencil(class LAMMPS *); virtual ~NStencil(); @@ -66,7 +66,10 @@ class NStencil : protected Pointers { double cutneighmaxsq; double *cuttypesq; double **cutneighsq; - + double **cutmultisq; + int n_multi_groups; + int *map_type_multi; + // data from NBin class int mbinx,mbiny,mbinz; diff --git a/src/nstencil_full_multi_2d.cpp b/src/nstencil_full_multi_2d.cpp index 61557f2d00..4c96f41a56 100644 --- a/src/nstencil_full_multi_2d.cpp +++ b/src/nstencil_full_multi_2d.cpp @@ -29,16 +29,16 @@ NStencilFullMulti2d::NStencilFullMulti2d(LAMMPS *lmp) : NStencil(lmp) {} void NStencilFullMulti2d::set_stencil_properties() { - int n = atom->ntypes; + int n = n_multi_groups; int i, j; // Always look up neighbor using full stencil and neighbor's bin - for (i = 1; i <= n; i++) { - for (j = 1; j <= n; j++) { + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { flag_half_multi[i][j] = 0; flag_skip_multi[i][j] = 0; - bin_type_multi[i][j] = j; + bin_group_multi[i][j] = j; } } } @@ -49,33 +49,33 @@ void NStencilFullMulti2d::set_stencil_properties() void NStencilFullMulti2d::create() { - int itype, jtype, bin_type, i, j, k, ns; - int n = atom->ntypes; + int igroup, jgroup, bin_group, i, j, k, ns; + int n = n_multi_groups; double cutsq; - for (itype = 1; itype <= n; itype++) { - for (jtype = 1; jtype <= n; jtype++) { - if (flag_skip_multi[itype][jtype]) continue; + for (igroup = 0; igroup < n; igroup++) { + for (jgroup = 0; jgroup < n; jgroup++) { + if (flag_skip_multi[igroup][jgroup]) continue; ns = 0; - sx = stencil_sx_multi[itype][jtype]; - sy = stencil_sy_multi[itype][jtype]; + sx = stencil_sx_multi[igroup][jgroup]; + sy = stencil_sy_multi[igroup][jgroup]; - mbinx = stencil_mbinx_multi[itype][jtype]; - mbiny = stencil_mbiny_multi[itype][jtype]; + mbinx = stencil_mbinx_multi[igroup][jgroup]; + mbiny = stencil_mbiny_multi[igroup][jgroup]; - bin_type = bin_type_multi[itype][jtype]; + bin_group = bin_group_multi[igroup][jgroup]; - cutsq = cutneighsq[itype][jtype]; + cutsq = cutmultisq[igroup][jgroup]; for (j = -sy; j <= sy; j++) for (i = -sx; i <= sx; i++) - if (bin_distance_multi(i,j,0,bin_type) < cutsq) - stencil_multi[itype][jtype][ns++] = j*mbinx + i; + if (bin_distance_multi(i,j,0,bin_group) < cutsq) + stencil_multi[igroup][jgroup][ns++] = j*mbinx + i; - nstencil_multi[itype][jtype] = ns; + nstencil_multi[igroup][jgroup] = ns; } } } diff --git a/src/nstencil_full_multi_3d.cpp b/src/nstencil_full_multi_3d.cpp index f4935b70e0..bc76433723 100644 --- a/src/nstencil_full_multi_3d.cpp +++ b/src/nstencil_full_multi_3d.cpp @@ -29,17 +29,17 @@ NStencilFullMulti3d::NStencilFullMulti3d(LAMMPS *lmp) : NStencil(lmp) {} void NStencilFullMulti3d::set_stencil_properties() { - int n = atom->ntypes; + int n = n_multi_groups; int i, j; // Always look up neighbor using full stencil and neighbor's bin // Stencil cutoff set by i-j cutoff - for (i = 1; i <= n; i++) { - for (j = 1; j <= n; j++) { + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { flag_half_multi[i][j] = 0; flag_skip_multi[i][j] = 0; - bin_type_multi[i][j] = j; + bin_group_multi[i][j] = j; } } } @@ -50,37 +50,37 @@ void NStencilFullMulti3d::set_stencil_properties() void NStencilFullMulti3d::create() { - int itype, jtype, bin_type, i, j, k, ns; - int n = atom->ntypes; + int igroup, jgroup, bin_group, i, j, k, ns; + int n = n_multi_groups; double cutsq; - for (itype = 1; itype <= n; itype++) { - for (jtype = 1; jtype <= n; jtype++) { - if (flag_skip_multi[itype][jtype]) continue; + for (igroup = 0; igroup < n; igroup++) { + for (jgroup = 0; jgroup < n; jgroup++) { + if (flag_skip_multi[igroup][jgroup]) continue; ns = 0; - sx = stencil_sx_multi[itype][jtype]; - sy = stencil_sy_multi[itype][jtype]; - sz = stencil_sz_multi[itype][jtype]; + sx = stencil_sx_multi[igroup][jgroup]; + sy = stencil_sy_multi[igroup][jgroup]; + sz = stencil_sz_multi[igroup][jgroup]; - mbinx = stencil_mbinx_multi[itype][jtype]; - mbiny = stencil_mbiny_multi[itype][jtype]; - mbinz = stencil_mbinz_multi[itype][jtype]; + mbinx = stencil_mbinx_multi[igroup][jgroup]; + mbiny = stencil_mbiny_multi[igroup][jgroup]; + mbinz = stencil_mbinz_multi[igroup][jgroup]; - bin_type = bin_type_multi[itype][jtype]; + bin_group = bin_group_multi[igroup][jgroup]; - cutsq = cutneighsq[itype][jtype]; + cutsq = cutmultisq[igroup][jgroup]; for (k = -sz; k <= sz; k++) for (j = -sy; j <= sy; j++) for (i = -sx; i <= sx; i++) - if (bin_distance_multi(i,j,k,bin_type) < cutsq) - stencil_multi[itype][jtype][ns++] = + if (bin_distance_multi(i,j,k,bin_group) < cutsq) + stencil_multi[igroup][jgroup][ns++] = k*mbiny*mbinx + j*mbinx + i; - nstencil_multi[itype][jtype] = ns; + nstencil_multi[igroup][jgroup] = ns; } } } diff --git a/src/nstencil_half_multi_2d.cpp b/src/nstencil_half_multi_2d.cpp index b407746602..8be8a1c151 100644 --- a/src/nstencil_half_multi_2d.cpp +++ b/src/nstencil_half_multi_2d.cpp @@ -30,26 +30,26 @@ NStencilHalfMulti2d::NStencilHalfMulti2d(LAMMPS *lmp) : void NStencilHalfMulti2d::set_stencil_properties() { - int n = atom->ntypes; + int n = n_multi_groups; int i, j; - // Cross types: use full stencil, looking one way through hierarchy + // Cross groups: use full stencil, looking one way through hierarchy // smaller -> larger => use full stencil in larger bin // larger -> smaller => no nstencil required // If cut offs are same, use half stencil - for (i = 1; i <= n; i++) { - for (j = 1; j <= n; j++) { - if(cutneighsq[i][i] > cutneighsq[j][j]) continue; + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + if(cutmultisq[i][i] > cutmultisq[j][j]) continue; flag_skip_multi[i][j] = 0; - if(cutneighsq[i][i] == cutneighsq[j][j]){ + if(cutmultisq[i][i] == cutmultisq[j][j]){ flag_half_multi[i][j] = 1; - bin_type_multi[i][j] = i; + bin_group_multi[i][j] = i; } else { flag_half_multi[i][j] = 0; - bin_type_multi[i][j] = j; + bin_group_multi[i][j] = j; } } } @@ -61,42 +61,42 @@ void NStencilHalfMulti2d::set_stencil_properties() void NStencilHalfMulti2d::create() { - int itype, jtype, bin_type, i, j, ns; - int n = atom->ntypes; + int igroup, jgroup, bin_group, i, j, ns; + int n = n_multi_groups; double cutsq; - for (itype = 1; itype <= n; itype++) { - for (jtype = 1; jtype <= n; jtype++) { - if (flag_skip_multi[itype][jtype]) continue; + for (igroup = 0; igroup < n; igroup++) { + for (jgroup = 0; jgroup < n; jgroup++) { + if (flag_skip_multi[igroup][jgroup]) continue; ns = 0; - sx = stencil_sx_multi[itype][jtype]; - sy = stencil_sy_multi[itype][jtype]; + sx = stencil_sx_multi[igroup][jgroup]; + sy = stencil_sy_multi[igroup][jgroup]; - mbinx = stencil_mbinx_multi[itype][jtype]; - mbiny = stencil_mbiny_multi[itype][jtype]; + mbinx = stencil_mbinx_multi[igroup][jgroup]; + mbiny = stencil_mbiny_multi[igroup][jgroup]; - bin_type = bin_type_multi[itype][jtype]; + bin_group = bin_group_multi[igroup][jgroup]; - cutsq = cutneighsq[itype][jtype]; + cutsq = cutmultisq[igroup][jgroup]; - if (flag_half_multi[itype][jtype]) { + if (flag_half_multi[igroup][jgroup]) { for (j = 0; j <= sy; j++) for (i = -sx; i <= sx; i++) if (j > 0 || (j == 0 && i > 0)) { - if (bin_distance_multi(i,j,0,bin_type) < cutsq) - stencil_multi[itype][jtype][ns++] = j*mbinx + i; + if (bin_distance_multi(i,j,0,bin_group) < cutsq) + stencil_multi[igroup][jgroup][ns++] = j*mbinx + i; } } else { for (j = -sy; j <= sy; j++) for (i = -sx; i <= sx; i++) - if (bin_distance_multi(i,j,0,bin_type) < cutsq) - stencil_multi[itype][jtype][ns++] = j*mbinx + i; + if (bin_distance_multi(i,j,0,bin_group) < cutsq) + stencil_multi[igroup][jgroup][ns++] = j*mbinx + i; } - nstencil_multi[itype][jtype] = ns; + nstencil_multi[igroup][jgroup] = ns; } } } diff --git a/src/nstencil_half_multi_2d_tri.cpp b/src/nstencil_half_multi_2d_tri.cpp index cd9f85da60..4a4e74c4a0 100755 --- a/src/nstencil_half_multi_2d_tri.cpp +++ b/src/nstencil_half_multi_2d_tri.cpp @@ -30,26 +30,26 @@ NStencilHalfMulti2dTri::NStencilHalfMulti2dTri(LAMMPS *lmp) : void NStencilHalfMulti2dTri::set_stencil_properties() { - int n = atom->ntypes; + int n = n_multi_groups; int i, j; - // Cross types: use full stencil, looking one way through hierarchy + // Cross groups: use full stencil, looking one way through hierarchy // smaller -> larger => use full stencil in larger bin // larger -> smaller => no nstencil required // If cut offs are same, use half stencil - for (i = 1; i <= n; i++) { - for (j = 1; j <= n; j++) { - if(cutneighsq[i][i] > cutneighsq[j][j]) continue; + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + if(cutmultisq[i][i] > cutmultisq[j][j]) continue; flag_skip_multi[i][j] = 0; - if(cutneighsq[i][i] == cutneighsq[j][j]){ + if(cutmultisq[i][i] == cutmultisq[j][j]){ flag_half_multi[i][j] = 1; - bin_type_multi[i][j] = i; + bin_group_multi[i][j] = i; } else { flag_half_multi[i][j] = 0; - bin_type_multi[i][j] = j; + bin_group_multi[i][j] = j; } } } @@ -61,40 +61,40 @@ void NStencilHalfMulti2dTri::set_stencil_properties() void NStencilHalfMulti2dTri::create() { - int itype, jtype, bin_type, i, j, ns; - int n = atom->ntypes; + int igroup, jgroup, bin_group, i, j, ns; + int n = n_multi_groups; double cutsq; - for (itype = 1; itype <= n; itype++) { - for (jtype = 1; jtype <= n; jtype++) { - if (flag_skip_multi[itype][jtype]) continue; + for (igroup = 0; igroup < n; igroup++) { + for (jgroup = 0; jgroup < n; jgroup++) { + if (flag_skip_multi[igroup][jgroup]) continue; ns = 0; - sx = stencil_sx_multi[itype][jtype]; - sy = stencil_sy_multi[itype][jtype]; + sx = stencil_sx_multi[igroup][jgroup]; + sy = stencil_sy_multi[igroup][jgroup]; - mbinx = stencil_mbinx_multi[itype][jtype]; - mbiny = stencil_mbiny_multi[itype][jtype]; + mbinx = stencil_mbinx_multi[igroup][jgroup]; + mbiny = stencil_mbiny_multi[igroup][jgroup]; - bin_type = bin_type_multi[itype][jtype]; + bin_group = bin_group_multi[igroup][jgroup]; - cutsq = cutneighsq[itype][jtype]; + cutsq = cutmultisq[igroup][jgroup]; - if (flag_half_multi[itype][jtype]) { + if (flag_half_multi[igroup][jgroup]) { for (j = 0; j <= sy; j++) for (i = -sx; i <= sx; i++) - if (bin_distance_multi(i,j,0,bin_type) < cutsq) - stencil_multi[itype][jtype][ns++] = j*mbinx + i; + if (bin_distance_multi(i,j,0,bin_group) < cutsq) + stencil_multi[igroup][jgroup][ns++] = j*mbinx + i; } else { for (j = -sy; j <= sy; j++) for (i = -sx; i <= sx; i++) - if (bin_distance_multi(i,j,0,bin_type) < cutsq) - stencil_multi[itype][jtype][ns++] = j*mbinx + i; + if (bin_distance_multi(i,j,0,bin_group) < cutsq) + stencil_multi[igroup][jgroup][ns++] = j*mbinx + i; } - nstencil_multi[itype][jtype] = ns; + nstencil_multi[igroup][jgroup] = ns; } } } diff --git a/src/nstencil_half_multi_3d.cpp b/src/nstencil_half_multi_3d.cpp index d66cfe98b9..ad4970603c 100644 --- a/src/nstencil_half_multi_3d.cpp +++ b/src/nstencil_half_multi_3d.cpp @@ -30,26 +30,26 @@ NStencilHalfMulti3d::NStencilHalfMulti3d(LAMMPS *lmp) : void NStencilHalfMulti3d::set_stencil_properties() { - int n = atom->ntypes; + int n = n_multi_groups; int i, j; - // Cross types: use full stencil, looking one way through hierarchy + // Cross groups: use full stencil, looking one way through hierarchy // smaller -> larger => use full stencil in larger bin // larger -> smaller => no nstencil required // If cut offs are same, use half stencil - for (i = 1; i <= n; i++) { - for (j = 1; j <= n; j++) { - if(cutneighsq[i][i] > cutneighsq[j][j]) continue; + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + if(cutmultisq[i][i] > cutmultisq[j][j]) continue; flag_skip_multi[i][j] = 0; - if(cutneighsq[i][i] == cutneighsq[j][j]){ + if(cutmultisq[i][i] == cutmultisq[j][j]){ flag_half_multi[i][j] = 1; - bin_type_multi[i][j] = i; + bin_group_multi[i][j] = i; } else { flag_half_multi[i][j] = 0; - bin_type_multi[i][j] = j; + bin_group_multi[i][j] = j; } } } @@ -61,48 +61,48 @@ void NStencilHalfMulti3d::set_stencil_properties() void NStencilHalfMulti3d::create() { - int itype, jtype, bin_type, i, j, k, ns; - int n = atom->ntypes; + int igroup, jgroup, bin_group, i, j, k, ns; + int n = n_multi_groups; double cutsq; - for (itype = 1; itype <= n; itype++) { - for (jtype = 1; jtype <= n; jtype++) { - if (flag_skip_multi[itype][jtype]) continue; + for (igroup = 0; igroup < n; igroup++) { + for (jgroup = 0; jgroup < n; jgroup++) { + if (flag_skip_multi[igroup][jgroup]) continue; ns = 0; - sx = stencil_sx_multi[itype][jtype]; - sy = stencil_sy_multi[itype][jtype]; - sz = stencil_sz_multi[itype][jtype]; + sx = stencil_sx_multi[igroup][jgroup]; + sy = stencil_sy_multi[igroup][jgroup]; + sz = stencil_sz_multi[igroup][jgroup]; - mbinx = stencil_mbinx_multi[itype][jtype]; - mbiny = stencil_mbiny_multi[itype][jtype]; - mbinz = stencil_mbinz_multi[itype][jtype]; + mbinx = stencil_mbinx_multi[igroup][jgroup]; + mbiny = stencil_mbiny_multi[igroup][jgroup]; + mbinz = stencil_mbinz_multi[igroup][jgroup]; - bin_type = bin_type_multi[itype][jtype]; + bin_group = bin_group_multi[igroup][jgroup]; - cutsq = cutneighsq[itype][jtype]; + cutsq = cutmultisq[igroup][jgroup]; - if (flag_half_multi[itype][jtype]) { + if (flag_half_multi[igroup][jgroup]) { for (k = 0; k <= sz; k++) for (j = -sy; j <= sy; j++) for (i = -sx; i <= sx; i++) if (k > 0 || j > 0 || (j == 0 && i > 0)) { - if (bin_distance_multi(i,j,k,bin_type) < cutsq) - stencil_multi[itype][jtype][ns++] = + if (bin_distance_multi(i,j,k,bin_group) < cutsq) + stencil_multi[igroup][jgroup][ns++] = k*mbiny*mbinx + j*mbinx + i; } } else { for (k = -sz; k <= sz; k++) for (j = -sy; j <= sy; j++) for (i = -sx; i <= sx; i++) - if (bin_distance_multi(i,j,k,bin_type) < cutsq) - stencil_multi[itype][jtype][ns++] = + if (bin_distance_multi(i,j,k,bin_group) < cutsq) + stencil_multi[igroup][jgroup][ns++] = k*mbiny*mbinx + j*mbinx + i; } - nstencil_multi[itype][jtype] = ns; + nstencil_multi[igroup][jgroup] = ns; } } } diff --git a/src/nstencil_half_multi_3d_tri.cpp b/src/nstencil_half_multi_3d_tri.cpp index e4a5747ca6..779f0865a5 100755 --- a/src/nstencil_half_multi_3d_tri.cpp +++ b/src/nstencil_half_multi_3d_tri.cpp @@ -30,26 +30,26 @@ NStencilHalfMulti3dTri::NStencilHalfMulti3dTri(LAMMPS *lmp) : void NStencilHalfMulti3dTri::set_stencil_properties() { - int n = atom->ntypes; + int n = n_multi_groups; int i, j; - // Cross types: use full stencil, looking one way through hierarchy + // Cross groups: use full stencil, looking one way through hierarchy // smaller -> larger => use full stencil in larger bin // larger -> smaller => no nstencil required // If cut offs are same, use half stencil - for (i = 1; i <= n; i++) { - for (j = 1; j <= n; j++) { - if(cutneighsq[i][i] > cutneighsq[j][j]) continue; + for (i = 0; i < n; i++) { + for (j = 0; j < n; j++) { + if(cutmultisq[i][i] > cutmultisq[j][j]) continue; flag_skip_multi[i][j] = 0; - if(cutneighsq[i][i] == cutneighsq[j][j]){ + if(cutmultisq[i][i] == cutmultisq[j][j]){ flag_half_multi[i][j] = 1; - bin_type_multi[i][j] = i; + bin_group_multi[i][j] = i; } else { flag_half_multi[i][j] = 0; - bin_type_multi[i][j] = j; + bin_group_multi[i][j] = j; } } } @@ -61,46 +61,46 @@ void NStencilHalfMulti3dTri::set_stencil_properties() void NStencilHalfMulti3dTri::create() { - int itype, jtype, bin_type, i, j, k, ns; - int n = atom->ntypes; + int igroup, jgroup, bin_group, i, j, k, ns; + int n = n_multi_groups; double cutsq; - for (itype = 1; itype <= n; itype++) { - for (jtype = 1; jtype <= n; jtype++) { - if (flag_skip_multi[itype][jtype]) continue; + for (igroup = 0; igroup < n; igroup++) { + for (jgroup = 0; jgroup < n; jgroup++) { + if (flag_skip_multi[igroup][jgroup]) continue; ns = 0; - sx = stencil_sx_multi[itype][jtype]; - sy = stencil_sy_multi[itype][jtype]; - sz = stencil_sz_multi[itype][jtype]; + sx = stencil_sx_multi[igroup][jgroup]; + sy = stencil_sy_multi[igroup][jgroup]; + sz = stencil_sz_multi[igroup][jgroup]; - mbinx = stencil_mbinx_multi[itype][jtype]; - mbiny = stencil_mbiny_multi[itype][jtype]; - mbinz = stencil_mbinz_multi[itype][jtype]; + mbinx = stencil_mbinx_multi[igroup][jgroup]; + mbiny = stencil_mbiny_multi[igroup][jgroup]; + mbinz = stencil_mbinz_multi[igroup][jgroup]; - bin_type = bin_type_multi[itype][jtype]; + bin_group = bin_group_multi[igroup][jgroup]; - cutsq = cutneighsq[itype][jtype]; + cutsq = cutmultisq[igroup][jgroup]; - if (flag_half_multi[itype][jtype]) { + if (flag_half_multi[igroup][jgroup]) { for (k = 0; k <= sz; k++) for (j = -sy; j <= sy; j++) for (i = -sx; i <= sx; i++) - if (bin_distance_multi(i,j,k,bin_type) < cutsq) - stencil_multi[itype][jtype][ns++] = + if (bin_distance_multi(i,j,k,bin_group) < cutsq) + stencil_multi[igroup][jgroup][ns++] = k*mbiny*mbinx + j*mbinx + i; } else { for (k = -sz; k <= sz; k++) for (j = -sy; j <= sy; j++) for (i = -sx; i <= sx; i++) - if (bin_distance_multi(i,j,k,bin_type) < cutsq) - stencil_multi[itype][jtype][ns++] = + if (bin_distance_multi(i,j,k,bin_group) < cutsq) + stencil_multi[igroup][jgroup][ns++] = k*mbiny*mbinx + j*mbinx + i; } - nstencil_multi[itype][jtype] = ns; + nstencil_multi[igroup][jgroup] = ns; } } }