diff --git a/src/USER-OMP/neigh_full_omp.cpp b/src/USER-OMP/neigh_full_omp.cpp index ea1eb7b0ff..f8036727dd 100644 --- a/src/USER-OMP/neigh_full_omp.cpp +++ b/src/USER-OMP/neigh_full_omp.cpp @@ -15,6 +15,8 @@ #include "neighbor_omp.h" #include "neigh_list.h" #include "atom.h" +#include "atom_vec.h" +#include "molecule.h" #include "comm.h" #include "domain.h" #include "group.h" @@ -32,6 +34,8 @@ void Neighbor::full_nsq_omp(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -39,20 +43,23 @@ void Neighbor::full_nsq_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,n,itype,jtype,which; + int i,j,n,itype,jtype,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + int nall = atom->nlocal + atom->nghost; - int molecular = atom->molecular; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -73,6 +80,11 @@ void Neighbor::full_nsq_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms, owned and ghost // skip i = j @@ -89,7 +101,13 @@ void Neighbor::full_nsq_omp(NeighList *list) rsq = delx*delx + dely*dely + delz*delz; if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -120,6 +138,8 @@ void Neighbor::full_nsq_ghost_omp(NeighList *list) { const int nlocal = atom->nlocal; const int nall = nlocal + atom->nghost; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -127,19 +147,22 @@ void Neighbor::full_nsq_ghost_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nall); - int i,j,n,itype,jtype,which; + int i,j,n,itype,jtype,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -160,6 +183,11 @@ void Neighbor::full_nsq_ghost_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms, owned and ghost // skip i = j @@ -177,7 +205,13 @@ void Neighbor::full_nsq_ghost_omp(NeighList *list) rsq = delx*delx + dely*dely + delz*delz; if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -224,6 +258,8 @@ void Neighbor::full_bin_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -231,19 +267,21 @@ void Neighbor::full_bin_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which; + int i,j,k,n,itype,jtype,ibin,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -266,6 +304,11 @@ void Neighbor::full_bin_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in surrounding bins in stencil including self // skip i = j @@ -286,7 +329,13 @@ void Neighbor::full_bin_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -322,6 +371,8 @@ void Neighbor::full_bin_ghost_omp(NeighList *list) const int nlocal = atom->nlocal; const int nall = nlocal + atom->nghost; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -329,20 +380,23 @@ void Neighbor::full_bin_ghost_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nall); - int i,j,k,n,itype,jtype,ibin,which; + int i,j,k,n,itype,jtype,ibin,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int xbin,ybin,zbin,xbin2,ybin2,zbin2; int *neighptr; - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -366,6 +420,11 @@ void Neighbor::full_bin_ghost_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in surrounding bins in stencil including self // when i is a ghost atom, must check if stencil bin is out of bounds @@ -388,7 +447,13 @@ void Neighbor::full_bin_ghost_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -448,6 +513,8 @@ void Neighbor::full_multi_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -455,22 +522,25 @@ void Neighbor::full_multi_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which,ns; + int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; double *cutsq,*distsq; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -492,6 +562,11 @@ void Neighbor::full_multi_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in other bins in stencil, including self // skip if i,j neighbor cutoff is less than bin distance @@ -517,6 +592,13 @@ void Neighbor::full_multi_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; which = find_special(special[i],nspecial[i],tag[j]); if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) diff --git a/src/USER-OMP/neigh_half_bin_omp.cpp b/src/USER-OMP/neigh_half_bin_omp.cpp index a8f515ab15..3974d9f6ab 100644 --- a/src/USER-OMP/neigh_half_bin_omp.cpp +++ b/src/USER-OMP/neigh_half_bin_omp.cpp @@ -15,6 +15,8 @@ #include "neighbor_omp.h" #include "neigh_list.h" #include "atom.h" +#include "atom_vec.h" +#include "molecule.h" #include "comm.h" #include "domain.h" #include "my_page.h" @@ -36,6 +38,8 @@ void Neighbor::half_bin_no_newton_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -43,21 +47,24 @@ void Neighbor::half_bin_no_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which; + int i,j,k,n,itype,jtype,ibin,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -78,6 +85,11 @@ void Neighbor::half_bin_no_newton_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in other bins in stencil including self // only store pair if i < j @@ -100,7 +112,13 @@ void Neighbor::half_bin_no_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -132,8 +150,14 @@ void Neighbor::half_bin_no_newton_omp(NeighList *list) void Neighbor::half_bin_no_newton_ghost_omp(NeighList *list) { + // bin local & ghost atoms + + bin_atoms(); + const int nlocal = atom->nlocal; const int nall = nlocal + atom->nghost; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -141,26 +165,25 @@ void Neighbor::half_bin_no_newton_ghost_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nall); - int i,j,k,n,itype,jtype,ibin,which; + int i,j,k,n,itype,jtype,ibin,which,imol,iatom; + tagint tagprev; int xbin,ybin,zbin,xbin2,ybin2,zbin2; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - // bin local & ghost atoms - - bin_atoms(); - // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -182,6 +205,11 @@ void Neighbor::half_bin_no_newton_ghost_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in other bins in stencil including self // when i is a ghost atom, must check if stencil bin is out of bounds @@ -208,7 +236,13 @@ void Neighbor::half_bin_no_newton_ghost_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -268,6 +302,8 @@ void Neighbor::half_bin_newton_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -275,21 +311,24 @@ void Neighbor::half_bin_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which; + int i,j,k,n,itype,jtype,ibin,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -310,6 +349,11 @@ void Neighbor::half_bin_newton_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // 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 @@ -334,7 +378,13 @@ void Neighbor::half_bin_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -359,7 +409,13 @@ void Neighbor::half_bin_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -394,6 +450,8 @@ void Neighbor::half_bin_newton_tri_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -401,21 +459,24 @@ void Neighbor::half_bin_newton_tri_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which; + int i,j,k,n,itype,jtype,ibin,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -436,6 +497,11 @@ void Neighbor::half_bin_newton_tri_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in bins in stencil // pairs for atoms j "below" i are excluded @@ -465,7 +531,13 @@ void Neighbor::half_bin_newton_tri_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; diff --git a/src/USER-OMP/neigh_half_multi_omp.cpp b/src/USER-OMP/neigh_half_multi_omp.cpp index c503b99cf6..ee5c3cf6f4 100644 --- a/src/USER-OMP/neigh_half_multi_omp.cpp +++ b/src/USER-OMP/neigh_half_multi_omp.cpp @@ -15,6 +15,8 @@ #include "neighbor_omp.h" #include "neigh_list.h" #include "atom.h" +#include "atom_vec.h" +#include "molecule.h" #include "comm.h" #include "domain.h" #include "my_page.h" @@ -37,6 +39,8 @@ void Neighbor::half_multi_no_newton_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -44,22 +48,25 @@ void Neighbor::half_multi_no_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which,ns; + int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; double *cutsq,*distsq; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -81,6 +88,11 @@ void Neighbor::half_multi_no_newton_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in other bins in stencil including self // only store pair if i < j @@ -108,7 +120,13 @@ void Neighbor::half_multi_no_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -143,6 +161,8 @@ void Neighbor::half_multi_newton_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -150,22 +170,25 @@ void Neighbor::half_multi_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which,ns; + int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; double *cutsq,*distsq; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -187,6 +210,11 @@ void Neighbor::half_multi_newton_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // 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 @@ -211,7 +239,13 @@ void Neighbor::half_multi_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -242,7 +276,13 @@ void Neighbor::half_multi_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -277,6 +317,8 @@ void Neighbor::half_multi_newton_tri_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -284,22 +326,25 @@ void Neighbor::half_multi_newton_tri_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,which,ns; + int i,j,k,n,itype,jtype,ibin,which,ns,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*s; double *cutsq,*distsq; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -321,6 +366,11 @@ void Neighbor::half_multi_newton_tri_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in bins, including self, in stencil // skip if i,j neighbor cutoff is less than bin distance @@ -357,7 +407,13 @@ void Neighbor::half_multi_newton_tri_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; diff --git a/src/USER-OMP/neigh_half_nsq_omp.cpp b/src/USER-OMP/neigh_half_nsq_omp.cpp index 6068ea174e..ab6e8e5d8c 100644 --- a/src/USER-OMP/neigh_half_nsq_omp.cpp +++ b/src/USER-OMP/neigh_half_nsq_omp.cpp @@ -15,6 +15,8 @@ #include "neighbor_omp.h" #include "neigh_list.h" #include "atom.h" +#include "atom_vec.h" +#include "molecule.h" #include "comm.h" #include "domain.h" #include "group.h" @@ -34,6 +36,8 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list) const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; const int nall = atom->nlocal + atom->nghost; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -41,19 +45,22 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,n,itype,jtype,which; + int i,j,n,itype,jtype,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -74,6 +81,11 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over remaining atoms, owned and ghost // only store pair if i < j @@ -90,7 +102,13 @@ void Neighbor::half_nsq_no_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -123,6 +141,8 @@ void Neighbor::half_nsq_no_newton_ghost_omp(NeighList *list) const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; const int nall = nlocal + atom->nghost; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -130,19 +150,22 @@ void Neighbor::half_nsq_no_newton_ghost_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nall); - int i,j,n,itype,jtype,which; + int i,j,n,itype,jtype,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -163,6 +186,11 @@ void Neighbor::half_nsq_no_newton_ghost_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over remaining atoms, owned and ghost // only store pair if i < j @@ -184,7 +212,13 @@ void Neighbor::half_nsq_no_newton_ghost_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -230,6 +264,8 @@ void Neighbor::half_nsq_newton_omp(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; #if defined(_OPENMP) @@ -237,22 +273,26 @@ void Neighbor::half_nsq_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,n,itype,jtype,itag,jtag,which; + int i,j,n,itype,jtype,itag,jtag,which,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + int nall = atom->nlocal + atom->nghost; - int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -272,6 +312,11 @@ void Neighbor::half_nsq_newton_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over remaining atoms, owned and ghost // itag = jtag is possible for long cutoffs that include images of self @@ -304,7 +349,13 @@ void Neighbor::half_nsq_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; diff --git a/src/USER-OMP/neigh_respa_omp.cpp b/src/USER-OMP/neigh_respa_omp.cpp index 2494c3abd9..6ac70cf9d7 100644 --- a/src/USER-OMP/neigh_respa_omp.cpp +++ b/src/USER-OMP/neigh_respa_omp.cpp @@ -15,6 +15,8 @@ #include "neighbor_omp.h" #include "neigh_list.h" #include "atom.h" +#include "atom_vec.h" +#include "molecule.h" #include "comm.h" #include "domain.h" #include "group.h" @@ -34,6 +36,8 @@ void Neighbor::respa_nsq_no_newton_omp(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; @@ -46,22 +50,26 @@ void Neighbor::respa_nsq_no_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,n,itype,jtype,n_inner,n_middle; + int i,j,n,itype,jtype,n_inner,n_middle,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + int nall = atom->nlocal + atom->nghost; - int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -107,6 +115,11 @@ void Neighbor::respa_nsq_no_newton_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over remaining atoms, owned and ghost @@ -122,7 +135,13 @@ void Neighbor::respa_nsq_no_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (minchange = domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -185,6 +204,8 @@ void Neighbor::respa_nsq_newton_omp(NeighList *list) { const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; const int bitmask = (includegroup) ? group->bitmask[includegroup] : 0; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; @@ -197,22 +218,26 @@ void Neighbor::respa_nsq_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,n,itype,jtype,itag,jtag,n_inner,n_middle; + int i,j,n,itype,jtype,itag,jtag,n_inner,n_middle,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + int nall = atom->nlocal + atom->nghost; - int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -259,6 +284,11 @@ void Neighbor::respa_nsq_newton_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over remaining atoms, owned and ghost @@ -290,7 +320,13 @@ void Neighbor::respa_nsq_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (minchange = domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -357,6 +393,8 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; @@ -369,21 +407,24 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,n_inner,n_middle; + int i,j,k,n,itype,jtype,ibin,n_inner,n_middle,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -432,6 +473,11 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list) ytmp = x[i][1]; ztmp = x[i][2]; ibin = coord2bin(x[i]); + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in surrounding bins in stencil including self // only store pair if i < j @@ -452,7 +498,13 @@ void Neighbor::respa_bin_no_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (minchange = domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -520,6 +572,8 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; @@ -532,21 +586,24 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,n_inner,n_middle; + int i,j,k,n,itype,jtype,ibin,n_inner,n_middle,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -594,6 +651,11 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // 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 @@ -618,7 +680,13 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (minchange = domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -656,7 +724,13 @@ void Neighbor::respa_bin_newton_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (minchange = domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j; @@ -724,6 +798,8 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list) bin_atoms(); const int nlocal = (includegroup) ? atom->nfirst : atom->nlocal; + const int molecular = atom->molecular; + const int moltemplate = (molecular == 2) ? 1 : 0; NEIGH_OMP_INIT; @@ -736,21 +812,24 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list) #endif NEIGH_OMP_SETUP(nlocal); - int i,j,k,n,itype,jtype,ibin,n_inner,n_middle; + int i,j,k,n,itype,jtype,ibin,n_inner,n_middle,imol,iatom; + tagint tagprev; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *neighptr,*neighptr_inner,*neighptr_middle; // loop over each atom, storing neighbors - tagint **special = atom->special; - int **nspecial = atom->nspecial; - tagint *tag = atom->tag; - double **x = atom->x; int *type = atom->type; int *mask = atom->mask; + tagint *tag = atom->tag; tagint *molecule = atom->molecule; - int molecular = atom->molecular; + tagint **special = atom->special; + int **nspecial = atom->nspecial; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -798,6 +877,11 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list) xtmp = x[i][0]; ytmp = x[i][1]; ztmp = x[i][2]; + if (moltemplate) { + imol = molindex[i]; + iatom = molatom[i]; + tagprev = tag[i] - iatom - 1; + } // loop over all atoms in bins in stencil // pairs for atoms j "below" i are excluded @@ -827,7 +911,13 @@ void Neighbor::respa_bin_newton_tri_omp(NeighList *list) if (rsq <= cutneighsq[itype][jtype]) { if (molecular) { - which = find_special(special[i],nspecial[i],tag[j]); + if (!moltemplate) + which = find_special(special[i],nspecial[i],tag[j]); + else if (imol >=0) + which = find_special(onemols[imol]->special[iatom], + onemols[imol]->nspecial[iatom], + tag[j]-tagprev); + else which = 0; if (which == 0) neighptr[n++] = j; else if (minchange = domain->minimum_image_check(delx,dely,delz)) neighptr[n++] = j;