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

Adding custom grouping option

Browse Source
This commit is contained in:
Joel Clemmer
2020-12-26 11:03:29 -07:00
parent b421c3d676
commit 2458eaf4f9
30 changed files with 599 additions and 551 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
src/USER-OMP/npair_full_multi_omp.cpp
View File

@ -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];

28
src/USER-OMP/npair_half_multi_newtoff_omp.cpp
View File

@ -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];

48
src/USER-OMP/npair_half_multi_newton_omp.cpp
View File

@ -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,6 +165,7 @@ void NPairHalfMultiNewtonOmp::build(NeighList *list)
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
@ -189,18 +192,19 @@ 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]) {
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];

30
src/USER-OMP/npair_half_multi_newton_tri_omp.cpp
View File

@ -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];

28
src/USER-OMP/npair_half_size_multi_newtoff_omp.cpp
View File

@ -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];

48
src/USER-OMP/npair_half_size_multi_newton_omp.cpp
View File

@ -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];

30
src/USER-OMP/npair_half_size_multi_newton_tri_omp.cpp
View File

@ -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];

6
src/comm_brick.cpp
View File

@ -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 {

45
src/nbin.cpp
View File

@ -44,7 +44,10 @@ NBin::NBin(LAMMPS *lmp) : Pointers(lmp)
atom2bin_multi = nullptr;
maxbins_multi = nullptr;
maxtypes = 0;
map_type_multi = nullptr;
cutmultisq = nullptr;
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]);
@ -118,6 +121,10 @@ void NBin::copy_neighbor_info()
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<int> ((x[0]-bboxhi[0])*bininvx_multi[it]) + nbinx_multi[it];
ix = static_cast<int> ((x[0]-bboxhi[0])*bininvx_multi[ig]) + nbinx_multi[ig];
else if (x[0] >= bboxlo[0]) {
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx_multi[it]);
ix = MIN(ix,nbinx_multi[it]-1);
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx_multi[ig]);
ix = MIN(ix,nbinx_multi[ig]-1);
} else
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx_multi[it]) - 1;
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx_multi[ig]) - 1;
if (x[1] >= bboxhi[1])
iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy_multi[it]) + nbiny_multi[it];
iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy_multi[ig]) + nbiny_multi[ig];
else if (x[1] >= bboxlo[1]) {
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy_multi[it]);
iy = MIN(iy,nbiny_multi[it]-1);
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy_multi[ig]);
iy = MIN(iy,nbiny_multi[ig]-1);
} else
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy_multi[it]) - 1;
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy_multi[ig]) - 1;
if (x[2] >= bboxhi[2])
iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz_multi[it]) + nbinz_multi[it];
iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz_multi[ig]) + nbinz_multi[ig];
else if (x[2] >= bboxlo[2]) {
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz_multi[it]);
iz = MIN(iz,nbinz_multi[it]-1);
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz_multi[ig]);
iz = MIN(iz,nbinz_multi[ig]-1);
} else
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz_multi[it]) - 1;
iz = static_cast<int> ((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;
}

10
src/nbin.h
View File

@ -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,14 +87,11 @@ 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 maxgroups; // size of multi arrays
int * maxbins_multi; // size of 2nd dimension of binhead_multi array
// methods

98
src/nbin_multi.cpp
View File

@ -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]);
@ -49,9 +49,11 @@ 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);
}

2
src/nbin_multi.h
View File

@ -38,7 +38,7 @@ class NBinMulti : public NBin {
private:
int itype_min();
int igroup_min();
};
}

2
src/neighbor.cpp
View File

@ -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;
}

40
src/npair.cpp
View File

@ -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<int> ((x[0]-bboxhi[0])*bininvx_multi[it]) + nbinx_multi[it];
ix = static_cast<int> ((x[0]-bboxhi[0])*bininvx_multi[ig]) + nbinx_multi[ig];
else if (x[0] >= bboxlo[0]) {
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx_multi[it]);
ix = MIN(ix,nbinx_multi[it]-1);
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx_multi[ig]);
ix = MIN(ix,nbinx_multi[ig]-1);
} else
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx_multi[it]) - 1;
ix = static_cast<int> ((x[0]-bboxlo[0])*bininvx_multi[ig]) - 1;
if (x[1] >= bboxhi[1])
iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy_multi[it]) + nbiny_multi[it];
iy = static_cast<int> ((x[1]-bboxhi[1])*bininvy_multi[ig]) + nbiny_multi[ig];
else if (x[1] >= bboxlo[1]) {
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy_multi[it]);
iy = MIN(iy,nbiny_multi[it]-1);
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy_multi[ig]);
iy = MIN(iy,nbiny_multi[ig]-1);
} else
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy_multi[it]) - 1;
iy = static_cast<int> ((x[1]-bboxlo[1])*bininvy_multi[ig]) - 1;
if (x[2] >= bboxhi[2])
iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz_multi[it]) + nbinz_multi[it];
iz = static_cast<int> ((x[2]-bboxhi[2])*bininvz_multi[ig]) + nbinz_multi[ig];
else if (x[2] >= bboxlo[2]) {
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz_multi[it]);
iz = MIN(iz,nbinz_multi[it]-1);
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz_multi[ig]);
iz = MIN(iz,nbinz_multi[ig]-1);
} else
iz = static_cast<int> ((x[2]-bboxlo[2])*bininvz_multi[it]) - 1;
iz = static_cast<int> ((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;
}

3
src/npair.h
View File

@ -48,6 +48,9 @@ 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

29
src/npair_full_multi.cpp
View File

@ -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,27 +78,28 @@ 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;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];

29
src/npair_half_multi_newtoff.cpp
View File

@ -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,29 +80,30 @@ 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;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];

49
src/npair_half_multi_newton.cpp
View File

@ -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,6 +109,7 @@ void NPairHalfMultiNewton::build(NeighList *list)
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
@ -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,6 +153,7 @@ void NPairHalfMultiNewton::build(NeighList *list)
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
@ -178,18 +180,19 @@ 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]) {
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];

31
src/npair_half_multi_newton_tri.cpp
View File

@ -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,6 +116,7 @@ void NPairHalfMultiNewtonTri::build(NeighList *list)
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];

29
src/npair_half_size_multi_newtoff.cpp
View File

@ -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,36 +64,37 @@ 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;
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];

49
src/npair_half_size_multi_newton.cpp
View File

@ -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,6 +98,7 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
@ -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,6 +135,7 @@ void NPairHalfSizeMultiNewton::build(NeighList *list)
}
}
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];
@ -153,18 +155,19 @@ 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]) {
jtype = type[j];
if (exclude && exclusion(i,j,itype,jtype,mask,molecule)) continue;
delx = xtmp - x[j][0];

31
src/npair_half_size_multi_newton_tri.cpp
View File

@ -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];

126
src/nstencil.cpp
View File

@ -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);
@ -156,6 +156,10 @@ void NStencil::copy_neighbor_info()
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<int> (stencil_range*bininvx_multi[bin_type]);
if (sx*binsizex_multi[bin_type] < stencil_range) sx++;
sy = static_cast<int> (stencil_range*bininvy_multi[bin_type]);
if (sy*binsizey_multi[bin_type] < stencil_range) sy++;
sz = static_cast<int> (stencil_range*bininvz_multi[bin_type]);
if (sz*binsizez_multi[bin_type] < stencil_range) sz++;
sx = static_cast<int> (stencil_range*bininvx_multi[bin_group]);
if (sx*binsizex_multi[bin_group] < stencil_range) sx++;
sy = static_cast<int> (stencil_range*bininvy_multi[bin_group]);
if (sy*binsizey_multi[bin_group] < stencil_range) sy++;
sz = static_cast<int> (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);

5
src/nstencil.h
View File

@ -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,6 +66,9 @@ class NStencil : protected Pointers {
double cutneighmaxsq;
double *cuttypesq;
double **cutneighsq;
double **cutmultisq;
int n_multi_groups;
int *map_type_multi;
// data from NBin class

36
src/nstencil_full_multi_2d.cpp
View File

@ -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;
}
}
}

40
src/nstencil_full_multi_3d.cpp
View File

@ -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;
}
}
}

50
src/nstencil_half_multi_2d.cpp
View File

@ -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;
}
}
}

50
src/nstencil_half_multi_2d_tri.cpp
View File

@ -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;
}
}
}

54
src/nstencil_half_multi_3d.cpp
View File

@ -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;
}
}
}

54
src/nstencil_half_multi_3d_tri.cpp
View File

@ -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;
}
}
}