From eb7d288360015fda9fd5c96e0b68d08f059af303 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 5 Apr 2011 19:10:43 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5896 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/ASPHERE/pair_gayberne.cpp | 9 +- src/ASPHERE/pair_resquared.cpp | 9 +- src/CLASS2/pair_lj_class2.cpp | 9 +- src/CLASS2/pair_lj_class2_coul_cut.cpp | 11 +-- src/CLASS2/pair_lj_class2_coul_long.cpp | 11 +-- src/COLLOID/pair_colloid.cpp | 9 +- src/COLLOID/pair_lubricate.cpp | 2 + src/COLLOID/pair_yukawa_colloid.cpp | 9 +- src/DIPOLE/pair_dipole_cut.cpp | 11 +-- src/GPU/pair_gayberne_gpu.cpp | 18 +--- src/GPU/pair_lj96_cut_gpu.cpp | 18 +--- src/GPU/pair_lj_charmm_coul_long_gpu.cpp | 22 ++--- src/GPU/pair_lj_cut_coul_cut_gpu.cpp | 22 ++--- src/GPU/pair_lj_cut_coul_long_gpu.cpp | 22 ++--- src/GPU/pair_lj_cut_gpu.cpp | 18 +--- src/GRANULAR/pair_gran_hertz_history.cpp | 1 + src/GRANULAR/pair_gran_hooke.cpp | 1 + src/GRANULAR/pair_gran_hooke_history.cpp | 1 + src/KSPACE/pair_born_coul_long.cpp | 11 +-- src/KSPACE/pair_buck_coul_long.cpp | 11 +-- src/KSPACE/pair_coul_long.cpp | 9 +- src/KSPACE/pair_lj_charmm_coul_long.cpp | 44 +++------- src/KSPACE/pair_lj_cut_coul_long.cpp | 48 ++++------- src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp | 11 +-- src/MANYBODY/pair_airebo.cpp | 2 + src/MANYBODY/pair_comb.cpp | 7 ++ src/MANYBODY/pair_eam.cpp | 2 + src/MANYBODY/pair_eim.cpp | 3 + src/MANYBODY/pair_sw.cpp | 3 + src/MANYBODY/pair_tersoff.cpp | 4 + src/MEAM/pair_meam.cpp | 29 +++++++ src/MEAM/pair_meam.h | 1 + src/MOLECULE/fix_bond_create.cpp | 1 + src/MOLECULE/fix_bond_swap.cpp | 1 + src/MOLECULE/pair_hbond_dreiding_lj.cpp | 9 +- src/MOLECULE/pair_hbond_dreiding_morse.cpp | 9 +- src/MOLECULE/pair_lj_charmm_coul_charmm.cpp | 11 +-- .../pair_lj_charmm_coul_charmm_implicit.cpp | 11 +-- src/OPT/pair_lj_charmm_coul_long_opt.h | 10 +-- src/OPT/pair_lj_cut_opt.h | 7 +- src/OPT/pair_morse_opt.h | 7 +- src/PERI/fix_peri_neigh.cpp | 2 + src/PERI/pair_peri_lps.cpp | 3 +- src/PERI/pair_peri_pmb.cpp | 3 +- src/REAX/pair_reax.cpp | 12 +-- src/USER-ACKLAND/compute_ackland_atom.cpp | 3 +- src/USER-CD-EAM/pair_cdeam.cpp | 3 + src/USER-CG-CMM/pair_cmm_common.cpp | 2 - src/USER-CG-CMM/pair_cmm_common.h | 56 ++++--------- src/USER-EFF/pair_eff_cut.cpp | 1 + src/USER-EWALDN/pair_buck_coul.cpp | 54 ++++++------ src/USER-EWALDN/pair_lj_coul.cpp | 56 ++++++------- src/USER-REAXC/pair_reax_c.cpp | 6 ++ src/compute.h | 4 + src/compute_centro_atom.cpp | 2 +- src/compute_cluster_atom.cpp | 2 +- src/compute_cna_atom.cpp | 3 + src/compute_coord_atom.cpp | 2 +- src/compute_group_group.cpp | 10 +-- src/compute_pair_local.cpp | 10 +-- src/compute_property_local.cpp | 10 +-- src/compute_rdf.cpp | 16 ++-- src/delete_atoms.cpp | 11 +-- src/delete_atoms.h | 4 + src/fix_orient_fcc.cpp | 1 + src/fix_shear_history.cpp | 1 + src/force.cpp | 1 + src/lmptype.h | 6 ++ src/neigh_derive.cpp | 35 ++++---- src/neigh_full.cpp | 58 +++++++------ src/neigh_half_bin.cpp | 39 ++++----- src/neigh_half_multi.cpp | 35 ++++---- src/neigh_half_nsq.cpp | 16 ++-- src/neigh_respa.cpp | 84 ++++++++++--------- src/neighbor.cpp | 14 ++-- src/neighbor.h | 66 +++++++-------- src/pair.h | 4 + src/pair_born.cpp | 9 +- src/pair_buck.cpp | 9 +- src/pair_buck_coul_cut.cpp | 11 +-- src/pair_coul_cut.cpp | 9 +- src/pair_coul_debye.cpp | 9 +- src/pair_dpd.cpp | 9 +- src/pair_dpd_tstat.cpp | 9 +- src/pair_gauss.cpp | 2 +- src/pair_lj96_cut.cpp | 37 ++------ src/pair_lj_cut.cpp | 36 ++------ src/pair_lj_cut_coul_cut.cpp | 11 +-- src/pair_lj_cut_coul_debye.cpp | 11 +-- src/pair_lj_expand.cpp | 9 +- src/pair_lj_gromacs.cpp | 9 +- src/pair_lj_gromacs_coul_gromacs.cpp | 11 +-- src/pair_lj_smooth.cpp | 9 +- src/pair_morse.cpp | 9 +- src/pair_soft.cpp | 9 +- src/pair_table.cpp | 9 +- src/pair_yukawa.cpp | 9 +- 97 files changed, 537 insertions(+), 788 deletions(-) diff --git a/src/ASPHERE/pair_gayberne.cpp b/src/ASPHERE/pair_gayberne.cpp index aa0d0c2a72..62fd2bbeae 100755 --- a/src/ASPHERE/pair_gayberne.cpp +++ b/src/ASPHERE/pair_gayberne.cpp @@ -91,7 +91,6 @@ void PairGayBerne::compute(int eflag, int vflag) double **tor = atom->torque; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -119,12 +118,8 @@ void PairGayBerne::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; // r12 = center to center vector diff --git a/src/ASPHERE/pair_resquared.cpp b/src/ASPHERE/pair_resquared.cpp index 4be78b42d6..e61f8c1f3e 100755 --- a/src/ASPHERE/pair_resquared.cpp +++ b/src/ASPHERE/pair_resquared.cpp @@ -96,7 +96,6 @@ void PairRESquared::compute(int eflag, int vflag) double **tor = atom->torque; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -120,12 +119,8 @@ void PairRESquared::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; // r12 = center to center vector diff --git a/src/CLASS2/pair_lj_class2.cpp b/src/CLASS2/pair_lj_class2.cpp index 9b6f41893f..8941239798 100644 --- a/src/CLASS2/pair_lj_class2.cpp +++ b/src/CLASS2/pair_lj_class2.cpp @@ -67,7 +67,6 @@ void PairLJClass2::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -89,12 +88,8 @@ void PairLJClass2::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/CLASS2/pair_lj_class2_coul_cut.cpp b/src/CLASS2/pair_lj_class2_coul_cut.cpp index 8c8dc7e228..2c91f9c0aa 100644 --- a/src/CLASS2/pair_lj_class2_coul_cut.cpp +++ b/src/CLASS2/pair_lj_class2_coul_cut.cpp @@ -73,7 +73,6 @@ void PairLJClass2CoulCut::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -98,13 +97,9 @@ void PairLJClass2CoulCut::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/CLASS2/pair_lj_class2_coul_long.cpp b/src/CLASS2/pair_lj_class2_coul_long.cpp index f714d1cc85..dd1bfc08f0 100644 --- a/src/CLASS2/pair_lj_class2_coul_long.cpp +++ b/src/CLASS2/pair_lj_class2_coul_long.cpp @@ -82,7 +82,6 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -107,13 +106,9 @@ void PairLJClass2CoulLong::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/COLLOID/pair_colloid.cpp b/src/COLLOID/pair_colloid.cpp index 8bd67abfac..e79891307c 100644 --- a/src/COLLOID/pair_colloid.cpp +++ b/src/COLLOID/pair_colloid.cpp @@ -84,7 +84,6 @@ void PairColloid::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -106,12 +105,8 @@ void PairColloid::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp index 43ec4031ec..8d86029221 100644 --- a/src/COLLOID/pair_lubricate.cpp +++ b/src/COLLOID/pair_lubricate.cpp @@ -117,6 +117,8 @@ void PairLubricate::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; diff --git a/src/COLLOID/pair_yukawa_colloid.cpp b/src/COLLOID/pair_yukawa_colloid.cpp index 790603df96..69c974a584 100644 --- a/src/COLLOID/pair_yukawa_colloid.cpp +++ b/src/COLLOID/pair_yukawa_colloid.cpp @@ -51,7 +51,6 @@ void PairYukawaColloid::compute(int eflag, int vflag) double **shape = atom->shape; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; int newton_pair = force->newton_pair; @@ -74,12 +73,8 @@ void PairYukawaColloid::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = 1.0; - else { - factor_coul = special_coul[j/nall]; - j %= nall; - } + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/DIPOLE/pair_dipole_cut.cpp b/src/DIPOLE/pair_dipole_cut.cpp index 20fc4076d0..a2e7c17bd4 100644 --- a/src/DIPOLE/pair_dipole_cut.cpp +++ b/src/DIPOLE/pair_dipole_cut.cpp @@ -81,7 +81,6 @@ void PairDipoleCut::compute(int eflag, int vflag) int *type = atom->type; double *dipole = atom->dipole; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -106,13 +105,9 @@ void PairDipoleCut::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j = j % nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/GPU/pair_gayberne_gpu.cpp b/src/GPU/pair_gayberne_gpu.cpp index eea45ad978..c8fc0cf35b 100644 --- a/src/GPU/pair_gayberne_gpu.cpp +++ b/src/GPU/pair_gayberne_gpu.cpp @@ -210,7 +210,6 @@ void PairGayBerneGPU::cpu_compute(int start, int eflag, int vflag) double **tor = atom->torque; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; inum = list->inum; @@ -237,12 +236,8 @@ void PairGayBerneGPU::cpu_compute(int start, int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; // r12 = center to center vector @@ -336,7 +331,6 @@ void PairGayBerneGPU::cpu_compute(int *nbors, int start, int eflag, int vflag) double **tor = atom->torque; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int stride = nlocal-start; double *special_lj = force->special_lj; @@ -360,12 +354,8 @@ void PairGayBerneGPU::cpu_compute(int *nbors, int start, int eflag, int vflag) for ( ; nbor < nbor_end; nbor += stride) { j = *nbor; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; // r12 = center to center vector diff --git a/src/GPU/pair_lj96_cut_gpu.cpp b/src/GPU/pair_lj96_cut_gpu.cpp index 8da2391068..2e059905c6 100644 --- a/src/GPU/pair_lj96_cut_gpu.cpp +++ b/src/GPU/pair_lj96_cut_gpu.cpp @@ -188,7 +188,6 @@ void PairLJ96CutGPU::cpu_compute(int start, int eflag, int vflag) { double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; inum = list->inum; @@ -209,12 +208,8 @@ void PairLJ96CutGPU::cpu_compute(int start, int eflag, int vflag) { for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -250,7 +245,6 @@ void PairLJ96CutGPU::cpu_compute(int start, int eflag, int vflag) { void PairLJ96CutGPU::cpu_compute(int *nbors, int start, int eflag, int vflag) { int i,j,itype,jtype; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int stride = nlocal-start; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double rsq,r2inv,r3inv,r6inv,forcelj,factor_lj; @@ -274,12 +268,8 @@ void PairLJ96CutGPU::cpu_compute(int *nbors, int start, int eflag, int vflag) { for (; nborq; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; double qqrd2e = force->qqrd2e; @@ -255,13 +254,9 @@ void PairLJCharmmCoulLongGPU::cpu_compute(int start, int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -365,7 +360,6 @@ void PairLJCharmmCoulLongGPU::cpu_compute(int *nbors, int start, int eflag, double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int stride = nlocal - start; double *special_coul = force->special_coul; double *special_lj = force->special_lj; @@ -386,13 +380,9 @@ void PairLJCharmmCoulLongGPU::cpu_compute(int *nbors, int start, int eflag, for (; nborq; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -221,13 +220,9 @@ void PairLJCutCoulCutGPU::cpu_compute(int start, int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -286,7 +281,6 @@ void PairLJCutCoulCutGPU::cpu_compute(int *nbors, int start, int eflag, double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int stride = nlocal-start; double *special_coul = force->special_coul; double *special_lj = force->special_lj; @@ -307,13 +301,9 @@ void PairLJCutCoulCutGPU::cpu_compute(int *nbors, int start, int eflag, for (; nborq; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; double qqrd2e = force->qqrd2e; @@ -248,13 +247,9 @@ void PairLJCutCoulLongGPU::cpu_compute(int start, int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -345,7 +340,6 @@ void PairLJCutCoulLongGPU::cpu_compute(int *nbors, int start, int eflag, double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int stride = nlocal-start; double *special_coul = force->special_coul; double *special_lj = force->special_lj; @@ -366,13 +360,9 @@ void PairLJCutCoulLongGPU::cpu_compute(int *nbors, int start, int eflag, for (; nborf; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; inum = list->inum; @@ -209,12 +208,8 @@ void PairLJCutGPU::cpu_compute(int start, int eflag, int vflag) { for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -249,7 +244,6 @@ void PairLJCutGPU::cpu_compute(int start, int eflag, int vflag) { void PairLJCutGPU::cpu_compute(int *nbors, int start, int eflag, int vflag) { int i,j,itype,jtype; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int stride = nlocal-start; double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; double rsq,r2inv,r6inv,forcelj,factor_lj; @@ -273,12 +267,8 @@ void PairLJCutGPU::cpu_compute(int *nbors, int start, int eflag, int vflag) { for (; nborq; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -114,13 +113,9 @@ void PairBornCoulLong::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/KSPACE/pair_buck_coul_long.cpp b/src/KSPACE/pair_buck_coul_long.cpp index 1cf9d0a49d..cd65fb0c51 100644 --- a/src/KSPACE/pair_buck_coul_long.cpp +++ b/src/KSPACE/pair_buck_coul_long.cpp @@ -82,7 +82,6 @@ void PairBuckCoulLong::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -107,13 +106,9 @@ void PairBuckCoulLong::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/KSPACE/pair_coul_long.cpp b/src/KSPACE/pair_coul_long.cpp index 20ad916416..1ddefad8e0 100644 --- a/src/KSPACE/pair_coul_long.cpp +++ b/src/KSPACE/pair_coul_long.cpp @@ -86,7 +86,6 @@ void PairCoulLong::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; @@ -110,12 +109,8 @@ void PairCoulLong::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = 1.0; - else { - factor_coul = special_coul[j/nall]; - j %= nall; - } + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/KSPACE/pair_lj_charmm_coul_long.cpp b/src/KSPACE/pair_lj_charmm_coul_long.cpp index b4a78d39ee..acedafe4b2 100644 --- a/src/KSPACE/pair_lj_charmm_coul_long.cpp +++ b/src/KSPACE/pair_lj_charmm_coul_long.cpp @@ -104,7 +104,6 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -129,13 +128,9 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -242,7 +237,6 @@ void PairLJCharmmCoulLong::compute_inner() double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -274,13 +268,9 @@ void PairLJCharmmCoulLong::compute_inner() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -332,7 +322,6 @@ void PairLJCharmmCoulLong::compute_middle() double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -370,14 +359,10 @@ void PairLJCharmmCoulLong::compute_middle() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } - delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; @@ -446,7 +431,6 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -478,13 +462,9 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/KSPACE/pair_lj_cut_coul_long.cpp b/src/KSPACE/pair_lj_cut_coul_long.cpp index 5273293575..26c8cb8ce5 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long.cpp @@ -79,7 +79,7 @@ PairLJCutCoulLong::~PairLJCutCoulLong() void PairLJCutCoulLong::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype,itable; + int i,ii,j,jj,inum,jnum,itype,jtype,itable; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; double fraction,table; double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj; @@ -96,7 +96,6 @@ void PairLJCutCoulLong::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -106,7 +105,7 @@ void PairLJCutCoulLong::compute(int eflag, int vflag) ilist = list->ilist; numneigh = list->numneigh; firstneigh = list->firstneigh; - + // loop over neighbors of my atoms for (ii = 0; ii < inum; ii++) { @@ -121,13 +120,9 @@ void PairLJCutCoulLong::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -222,7 +217,6 @@ void PairLJCutCoulLong::compute_inner() double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -254,13 +248,9 @@ void PairLJCutCoulLong::compute_inner() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -312,7 +302,6 @@ void PairLJCutCoulLong::compute_middle() double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -349,13 +338,9 @@ void PairLJCutCoulLong::compute_middle() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -418,7 +403,6 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -450,13 +434,9 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp index 73c180b795..1aab94395f 100644 --- a/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp @@ -88,7 +88,6 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -117,13 +116,9 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 96ef31bbb7..168c9d7136 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -351,6 +351,7 @@ void PairAIREBO::REBO_neigh() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtype = map[type[j]]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -534,6 +535,7 @@ void PairAIREBO::FLJ(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtag = tag[j]; if (itag > jtag) { diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp index ea415b968d..3e4d95fa60 100644 --- a/src/MANYBODY/pair_comb.cpp +++ b/src/MANYBODY/pair_comb.cpp @@ -182,6 +182,7 @@ void PairComb::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtag = tag[j]; if (itag > jtag) { @@ -261,6 +262,7 @@ void PairComb::compute(int eflag, int vflag) if (cor_flag) { for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtype = map[type[j]]; iparam_ij = elem2param[itype][jtype][jtype]; @@ -281,6 +283,7 @@ void PairComb::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtype = map[type[j]]; iparam_ij = elem2param[itype][jtype][jtype]; @@ -303,6 +306,7 @@ void PairComb::compute(int eflag, int vflag) for (kk = 0; kk < jnum; kk++) { if (jj == kk) continue; k = jlist[kk]; + k &= NEIGHMASK; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; @@ -347,6 +351,7 @@ void PairComb::compute(int eflag, int vflag) for (kk = 0; kk < jnum; kk++) { if (jj == kk) continue; k = jlist[kk]; + k &= NEIGHMASK; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; @@ -1678,6 +1683,7 @@ double PairComb::yasu_char(double *qf_fix, int &igroup) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtype = map[type[j]]; jq = q[j]; @@ -1719,6 +1725,7 @@ double PairComb::yasu_char(double *qf_fix, int &igroup) for (kk = 0; kk < jnum; kk++) { if (jj == kk) continue; k = jlist[kk]; + k &= NEIGHMASK; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp index 1da055a96b..b89a397bfb 100644 --- a/src/MANYBODY/pair_eam.cpp +++ b/src/MANYBODY/pair_eam.cpp @@ -177,6 +177,7 @@ void PairEAM::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -242,6 +243,7 @@ void PairEAM::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/MANYBODY/pair_eim.cpp b/src/MANYBODY/pair_eim.cpp index cc7b702da8..3a99d59a85 100644 --- a/src/MANYBODY/pair_eim.cpp +++ b/src/MANYBODY/pair_eim.cpp @@ -166,6 +166,7 @@ void PairEIM::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtype = type[j]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -205,6 +206,7 @@ void PairEIM::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtype = type[j]; delx = xtmp - x[j][0]; @@ -259,6 +261,7 @@ void PairEIM::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtype = type[j]; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/MANYBODY/pair_sw.cpp b/src/MANYBODY/pair_sw.cpp index ed062a5c76..555a7622a4 100755 --- a/src/MANYBODY/pair_sw.cpp +++ b/src/MANYBODY/pair_sw.cpp @@ -113,6 +113,7 @@ void PairSW::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtag = tag[j]; if (itag > jtag) { @@ -152,6 +153,7 @@ void PairSW::compute(int eflag, int vflag) for (jj = 0; jj < jnumm1; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtype = map[type[j]]; ijparam = elem2param[itype][jtype][jtype]; delr1[0] = x[j][0] - xtmp; @@ -162,6 +164,7 @@ void PairSW::compute(int eflag, int vflag) for (kk = jj+1; kk < jnum; kk++) { k = jlist[kk]; + k &= NEIGHMASK; ktype = map[type[k]]; ikparam = elem2param[itype][ktype][ktype]; ijkparam = elem2param[itype][jtype][ktype]; diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index 0ca96b583e..579356f817 100755 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -116,6 +116,7 @@ void PairTersoff::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtag = tag[j]; if (itag > jtag) { @@ -156,6 +157,7 @@ void PairTersoff::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtype = map[type[j]]; iparam_ij = elem2param[itype][jtype][jtype]; @@ -172,6 +174,7 @@ void PairTersoff::compute(int eflag, int vflag) for (kk = 0; kk < jnum; kk++) { if (jj == kk) continue; k = jlist[kk]; + k &= NEIGHMASK; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; @@ -203,6 +206,7 @@ void PairTersoff::compute(int eflag, int vflag) for (kk = 0; kk < jnum; kk++) { if (jj == kk) continue; k = jlist[kk]; + k &= NEIGHMASK; ktype = map[type[k]]; iparam_ijk = elem2param[itype][jtype][ktype]; diff --git a/src/MEAM/pair_meam.cpp b/src/MEAM/pair_meam.cpp index d922a75278..d35cecc853 100644 --- a/src/MEAM/pair_meam.cpp +++ b/src/MEAM/pair_meam.cpp @@ -180,6 +180,14 @@ void PairMEAM::compute(int eflag, int vflag) numneigh_full = listfull->numneigh; firstneigh_full = listfull->firstneigh; + // strip neighbor lists of any special bond flags before using with MEAM + // necessary before doing neigh_f2c and neigh_c2f conversions each step + + if (neighbor->ago == 0) { + neigh_strip(inum_half,ilist_half,numneigh_half,firstneigh_half); + neigh_strip(inum_half,ilist_half,numneigh_full,firstneigh_full); + } + // check size of scrfcn based on half neighbor list int nlocal = atom->nlocal; @@ -887,6 +895,27 @@ double PairMEAM::memory_usage() return bytes; } +/* ---------------------------------------------------------------------- + strip special bond flags from neighbor list entries + are not used with MEAM + need to do here so Fortran lib doesn't see them + done once per reneighbor so that neigh_f2c and neigh_c2f don't see them +------------------------------------------------------------------------- */ + +void PairMEAM::neigh_strip(int inum, int *ilist, + int *numneigh, int **firstneigh) +{ + int i,j,ii,jnum; + int *jlist; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + for (j = 0; j < jnum; j++) jlist[j] &= NEIGHMASK; + } +} + /* ---------------------------------------------------------------------- toggle neighbor list indices between zero- and one-based values needed for access by MEAM Fortran library diff --git a/src/MEAM/pair_meam.h b/src/MEAM/pair_meam.h index 7c993b24c1..267f0e3f34 100644 --- a/src/MEAM/pair_meam.h +++ b/src/MEAM/pair_meam.h @@ -94,6 +94,7 @@ class PairMEAM : public Pair { void allocate(); void read_files(char *, char *); + void neigh_strip(int, int *, int *, int **); void neigh_f2c(int, int *, int *, int **); void neigh_c2f(int, int *, int *, int **); }; diff --git a/src/MOLECULE/fix_bond_create.cpp b/src/MOLECULE/fix_bond_create.cpp index ec61cda059..68689c4209 100755 --- a/src/MOLECULE/fix_bond_create.cpp +++ b/src/MOLECULE/fix_bond_create.cpp @@ -328,6 +328,7 @@ void FixBondCreate::post_integrate() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; if (!(mask[j] & groupbit)) continue; jtype = type[j]; diff --git a/src/MOLECULE/fix_bond_swap.cpp b/src/MOLECULE/fix_bond_swap.cpp index cc49eb7c93..daf9b13501 100644 --- a/src/MOLECULE/fix_bond_swap.cpp +++ b/src/MOLECULE/fix_bond_swap.cpp @@ -241,6 +241,7 @@ void FixBondSwap::pre_neighbor() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; if (j >= nlocal) continue; if ((mask[j] & groupbit) == 0) continue; if (molecule[i] != molecule[j]) continue; diff --git a/src/MOLECULE/pair_hbond_dreiding_lj.cpp b/src/MOLECULE/pair_hbond_dreiding_lj.cpp index b97e2969c1..81eb85005f 100644 --- a/src/MOLECULE/pair_hbond_dreiding_lj.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_lj.cpp @@ -98,7 +98,6 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) int *type = atom->type; int **nspecial = atom->nspecial; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; inum = list->inum; @@ -123,11 +122,9 @@ void PairHbondDreidingLJ::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - if (j < nall) factor_hb = 1.0; - else { - factor_hb = special_lj[j/nall]; - j %= nall; - } + factor_hb = special_lj[sbmask(j)]; + j &= NEIGHMASK; + jtype = type[j]; if (!acceptor[jtype]) continue; diff --git a/src/MOLECULE/pair_hbond_dreiding_morse.cpp b/src/MOLECULE/pair_hbond_dreiding_morse.cpp index 084d84c5d1..2740e1b61f 100644 --- a/src/MOLECULE/pair_hbond_dreiding_morse.cpp +++ b/src/MOLECULE/pair_hbond_dreiding_morse.cpp @@ -66,7 +66,6 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) int *type = atom->type; int **nspecial = atom->nspecial; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; inum = list->inum; @@ -91,11 +90,9 @@ void PairHbondDreidingMorse::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - if (j < nall) factor_hb = 1.0; - else { - factor_hb = special_lj[j/nall]; - j %= nall; - } + factor_hb = special_lj[sbmask(j)]; + j &= NEIGHMASK; + jtype = type[j]; if (!acceptor[jtype]) continue; diff --git a/src/MOLECULE/pair_lj_charmm_coul_charmm.cpp b/src/MOLECULE/pair_lj_charmm_coul_charmm.cpp index b8f2ffa62d..bf8e713a01 100644 --- a/src/MOLECULE/pair_lj_charmm_coul_charmm.cpp +++ b/src/MOLECULE/pair_lj_charmm_coul_charmm.cpp @@ -85,7 +85,6 @@ void PairLJCharmmCoulCharmm::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -110,13 +109,9 @@ void PairLJCharmmCoulCharmm::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/MOLECULE/pair_lj_charmm_coul_charmm_implicit.cpp b/src/MOLECULE/pair_lj_charmm_coul_charmm_implicit.cpp index d81b4f9d17..990392fa04 100644 --- a/src/MOLECULE/pair_lj_charmm_coul_charmm_implicit.cpp +++ b/src/MOLECULE/pair_lj_charmm_coul_charmm_implicit.cpp @@ -47,7 +47,6 @@ void PairLJCharmmCoulCharmmImplicit::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -72,13 +71,9 @@ void PairLJCharmmCoulCharmmImplicit::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/OPT/pair_lj_charmm_coul_long_opt.h b/src/OPT/pair_lj_charmm_coul_long_opt.h index bc1171e5ae..583f2d57e0 100644 --- a/src/OPT/pair_lj_charmm_coul_long_opt.h +++ b/src/OPT/pair_lj_charmm_coul_long_opt.h @@ -82,7 +82,6 @@ void PairLJCharmmCoulLongOpt::eval() double* __restrict__ q = atom->q; int* __restrict__ type = atom->type; int nlocal = atom->nlocal; - int nall = atom->nlocal + atom->nghost; double* __restrict__ special_coul = force->special_coul; double* __restrict__ special_lj = force->special_lj; double qqrd2e = force->qqrd2e; @@ -134,7 +133,7 @@ void PairLJCharmmCoulLongOpt::eval() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - if (j < nall) { + if (j <= NEIGHMASK) { double delx = xtmp - xx[j].x; double dely = ytmp - xx[j].y; double delz = ztmp - xx[j].z; @@ -220,9 +219,10 @@ void PairLJCharmmCoulLongOpt::eval() } } else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; + double delx = xtmp - xx[j].x; double dely = ytmp - xx[j].y; double delz = ztmp - xx[j].z; diff --git a/src/OPT/pair_lj_cut_opt.h b/src/OPT/pair_lj_cut_opt.h index c5122f6d72..dce90b0249 100644 --- a/src/OPT/pair_lj_cut_opt.h +++ b/src/OPT/pair_lj_cut_opt.h @@ -61,7 +61,6 @@ void PairLJCutOpt::eval() double** __restrict__ f = atom->f; int* __restrict__ type = atom->type; int nlocal = atom->nlocal; - int nall = atom->nlocal + atom->nghost; double* __restrict__ special_lj = force->special_lj; inum = list->inum; @@ -109,7 +108,7 @@ void PairLJCutOpt::eval() j = jlist[jj]; double factor_lj; - if (j < nall) { + if (j <= NEIGHMASK) { double delx = xtmp - xx[j].x; double dely = ytmp - xx[j].y; double delz = ztmp - xx[j].z; @@ -142,8 +141,8 @@ void PairLJCutOpt::eval() } } else { - factor_lj = special_lj[j/nall]; - j = j % nall; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; double delx = xtmp - xx[j].x; double dely = ytmp - xx[j].y; diff --git a/src/OPT/pair_morse_opt.h b/src/OPT/pair_morse_opt.h index a8bfa1ec49..0e92408e1d 100644 --- a/src/OPT/pair_morse_opt.h +++ b/src/OPT/pair_morse_opt.h @@ -62,7 +62,6 @@ void PairMorseOpt::eval() double** __restrict__ f = atom->f; int* __restrict__ type = atom->type; int nlocal = atom->nlocal; - int nall = atom->nlocal + atom->nghost; double* __restrict__ special_lj = force->special_lj; inum = list->inum; @@ -110,7 +109,7 @@ void PairMorseOpt::eval() j = jlist[jj]; double factor_lj; - if (j < nall) { + if (j <= NEIGHMASK) { double delx = xtmp - xx[j].x; double dely = ytmp - xx[j].y; double delz = ztmp - xx[j].z; @@ -141,8 +140,8 @@ void PairMorseOpt::eval() } } else { - factor_lj = special_lj[j/nall]; - j = j % nall; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; double delx = xtmp - xx[j].x; double dely = ytmp - xx[j].y; diff --git a/src/PERI/fix_peri_neigh.cpp b/src/PERI/fix_peri_neigh.cpp index f006809d7d..666a276973 100644 --- a/src/PERI/fix_peri_neigh.cpp +++ b/src/PERI/fix_peri_neigh.cpp @@ -177,6 +177,7 @@ void FixPeriNeigh::setup(int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -224,6 +225,7 @@ void FixPeriNeigh::setup(int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/PERI/pair_peri_lps.cpp b/src/PERI/pair_peri_lps.cpp index a3d353c2f0..c4a5cdcae3 100644 --- a/src/PERI/pair_peri_lps.cpp +++ b/src/PERI/pair_peri_lps.cpp @@ -119,7 +119,6 @@ void PairPeriLPS::compute(int eflag, int vflag) // short-range forces - int nall = atom->nlocal + atom->nghost; int newton_pair = force->newton_pair; int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic; @@ -147,7 +146,7 @@ void PairPeriLPS::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - j %= nall; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/PERI/pair_peri_pmb.cpp b/src/PERI/pair_peri_pmb.cpp index fcbc8d1546..3ee46cb938 100644 --- a/src/PERI/pair_peri_pmb.cpp +++ b/src/PERI/pair_peri_pmb.cpp @@ -111,7 +111,6 @@ void PairPeriPMB::compute(int eflag, int vflag) // short-range forces - int nall = atom->nlocal + atom->nghost; int newton_pair = force->newton_pair; int periodic = (domain->xperiodic || domain->yperiodic || domain->zperiodic); @@ -137,7 +136,7 @@ void PairPeriPMB::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - j %= nall; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/REAX/pair_reax.cpp b/src/REAX/pair_reax.cpp index c162212ed0..7759931833 100644 --- a/src/REAX/pair_reax.cpp +++ b/src/REAX/pair_reax.cpp @@ -309,9 +309,7 @@ void PairREAX::write_reax_vlist() int nvpair, nvlself, nvpairmax; int nbond; int inum,jnum; - int *ilist; - int *jlist; - int *numneigh,**firstneigh; + int *ilist,*jlist,*numneigh,**firstneigh; double delr2; double delx, dely, delz; double rtmp[3]; @@ -343,6 +341,8 @@ void PairREAX::write_reax_vlist() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; + xjtmp = x[j][0]; yjtmp = x[j][1]; zjtmp = x[j][2]; @@ -755,9 +755,7 @@ void PairREAX::compute_charge(double &energy_charge_equilibration) double sw; double rtmp[3]; int inum,jnum; - int *ilist; - int *jlist; - int *numneigh,**firstneigh; + int *ilist,*jlist,*numneigh,**firstneigh; double **x = atom->x; double *q = atom->q; @@ -811,6 +809,8 @@ void PairREAX::compute_charge(double &energy_charge_equilibration) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; + xjtmp = x[j][0]; yjtmp = x[j][1]; zjtmp = x[j][2]; diff --git a/src/USER-ACKLAND/compute_ackland_atom.cpp b/src/USER-ACKLAND/compute_ackland_atom.cpp index f802e4a43e..fa50b5f9d2 100644 --- a/src/USER-ACKLAND/compute_ackland_atom.cpp +++ b/src/USER-ACKLAND/compute_ackland_atom.cpp @@ -126,7 +126,6 @@ void ComputeAcklandAtom::compute_peratom() double **x = atom->x; int *mask = atom->mask; - int nall = atom->nlocal + atom->nghost; double cutsq = force->pair->cutforce * force->pair->cutforce; for (ii = 0; ii < inum; ii++) { @@ -159,7 +158,7 @@ void ComputeAcklandAtom::compute_peratom() n = 0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - if (j >= nall) j %= nall; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/USER-CD-EAM/pair_cdeam.cpp b/src/USER-CD-EAM/pair_cdeam.cpp index efa829f80d..5dc2ede1a1 100644 --- a/src/USER-CD-EAM/pair_cdeam.cpp +++ b/src/USER-CD-EAM/pair_cdeam.cpp @@ -136,6 +136,7 @@ void PairCDEAM::compute(int eflag, int vflag) for(jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -225,6 +226,7 @@ void PairCDEAM::compute(int eflag, int vflag) for(jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; jtype = type[j]; if(itype == jtype) continue; @@ -304,6 +306,7 @@ void PairCDEAM::compute(int eflag, int vflag) for(jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/USER-CG-CMM/pair_cmm_common.cpp b/src/USER-CG-CMM/pair_cmm_common.cpp index fb8bf72a82..8dc7e9cace 100644 --- a/src/USER-CG-CMM/pair_cmm_common.cpp +++ b/src/USER-CG-CMM/pair_cmm_common.cpp @@ -474,5 +474,3 @@ double PairCMMCommon::eval_single(int coul_type, int i, int j, int itype, int jt fforce = factor_lj*lj_force + factor_coul*coul_force; return factor_lj*lj_erg + factor_coul*coul_erg; } - -/* ------------------------------------------------------------------------ */ diff --git a/src/USER-CG-CMM/pair_cmm_common.h b/src/USER-CG-CMM/pair_cmm_common.h index 6d5d1a5839..43132abac1 100644 --- a/src/USER-CG-CMM/pair_cmm_common.h +++ b/src/USER-CG-CMM/pair_cmm_common.h @@ -185,10 +185,10 @@ namespace LAMMPS_NS { const double * const q = atom->q; const int * const type = atom->type; const int nlocal = atom->nlocal; - const int nall = nlocal + atom->nghost; const double * const special_lj = force->special_lj; const double * const special_coul = force->special_coul; const double qqrd2e = force->qqrd2e; + double factor_lj,factor_coul; const int inum = list->inum; const int * const ilist = list->ilist; @@ -211,15 +211,9 @@ namespace LAMMPS_NS { for (jj = 0; jj < jnum; jj++) { int j2 = jlist[jj]; - - double factor_lj = 1.0; - double factor_coul = 1.0; - if (j2 >= nall) { - factor_lj = special_lj[j2/nall]; - if (COUL_TYPE != CG_COUL_NONE) factor_coul=special_coul[j2/nall]; - j2 %= nall; - } - const int j = j2; + factor_lj = special_lj[sbmask(j2)]; + factor_coul = special_coul[sbmask(j2)]; + const int j = j2 & NEIGHMASK; const double delx = xtmp - x[j][0]; const double dely = ytmp - x[j][1]; @@ -323,10 +317,10 @@ namespace LAMMPS_NS { const double * const q = atom->q; const int * const type = atom->type; const int nlocal = atom->nlocal; - const int nall = nlocal + atom->nghost; const double * const special_lj = force->special_lj; const double * const special_coul = force->special_coul; const double qqrd2e = force->qqrd2e; + double factor_lj,factor_coul; const int inum = listinner->inum; const int * const ilist = listinner->ilist; @@ -355,15 +349,9 @@ namespace LAMMPS_NS { for (jj = 0; jj < jnum; jj++) { int j2 = jlist[jj]; - - double factor_lj = 1.0; - double factor_coul = 1.0; - if (j2 >= nall) { - factor_lj = special_lj[j2/nall]; - if (COUL_TYPE != CG_COUL_NONE) factor_coul=special_coul[j2/nall]; - j2 %= nall; - } - const int j = j2; + factor_lj = special_lj[sbmask(j2)]; + factor_coul = special_coul[sbmask(j2)]; + const int j = j2 & NEIGHMASK; const double delx = xtmp - x[j][0]; const double dely = ytmp - x[j][1]; @@ -427,10 +415,10 @@ namespace LAMMPS_NS { const double * const q = atom->q; const int * const type = atom->type; const int nlocal = atom->nlocal; - const int nall = nlocal + atom->nghost; const double * const special_lj = force->special_lj; const double * const special_coul = force->special_coul; const double qqrd2e = force->qqrd2e; + double factor_lj,factor_coul; const int inum = listmiddle->inum; const int * const ilist = listmiddle->ilist; @@ -464,15 +452,9 @@ namespace LAMMPS_NS { for (jj = 0; jj < jnum; jj++) { int j2 = jlist[jj]; - - double factor_lj = 1.0; - double factor_coul = 1.0; - if (j2 >= nall) { - factor_lj = special_lj[j2/nall]; - if (COUL_TYPE != CG_COUL_NONE) factor_coul=special_coul[j2/nall]; - j2 %= nall; - } - const int j = j2; + factor_lj = special_lj[sbmask(j2)]; + factor_coul = special_coul[sbmask(j2)]; + const int j = j2 & NEIGHMASK; const double delx = xtmp - x[j][0]; const double dely = ytmp - x[j][1]; @@ -545,10 +527,10 @@ namespace LAMMPS_NS { const double * const q = atom->q; const int * const type = atom->type; const int nlocal = atom->nlocal; - const int nall = nlocal + atom->nghost; const double * const special_lj = force->special_lj; const double * const special_coul = force->special_coul; const double qqrd2e = force->qqrd2e; + double factor_lj,factor_coul; const int inum = listouter->inum; const int * const ilist = listouter->ilist; @@ -577,15 +559,9 @@ namespace LAMMPS_NS { for (jj = 0; jj < jnum; jj++) { int j2 = jlist[jj]; - - double factor_lj = 1.0; - double factor_coul = 1.0; - if (j2 >= nall) { - factor_lj = special_lj[j2/nall]; - if (COUL_TYPE != CG_COUL_NONE) factor_coul=special_coul[j2/nall]; - j2 %= nall; - } - const int j = j2; + factor_lj = special_lj[sbmask(j2)]; + factor_coul = special_coul[sbmask(j2)]; + const int j = j2 & NEIGHMASK; const double delx = xtmp - x[j][0]; const double dely = ytmp - x[j][1]; diff --git a/src/USER-EFF/pair_eff_cut.cpp b/src/USER-EFF/pair_eff_cut.cpp index 5cbebbf4ef..47c99d2838 100644 --- a/src/USER-EFF/pair_eff_cut.cpp +++ b/src/USER-EFF/pair_eff_cut.cpp @@ -169,6 +169,7 @@ void PairEffCut::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/USER-EWALDN/pair_buck_coul.cpp b/src/USER-EWALDN/pair_buck_coul.cpp index dc00608785..d2f3fdf142 100644 --- a/src/USER-EWALDN/pair_buck_coul.cpp +++ b/src/USER-EWALDN/pair_buck_coul.cpp @@ -453,7 +453,6 @@ void PairBuckCoul::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -480,8 +479,9 @@ void PairBuckCoul::compute(int eflag, int vflag) jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i]; for (; jneigh> ncoulshiftbits; register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; - if (ni < 0) { + if (ni == 0) { force_coul = qiqj*(ftable[k]+f*dftable[k]); if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); } @@ -532,7 +532,7 @@ void PairBuckCoul::compute(int eflag, int vflag) if (order6) { // long-range register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2)*buckci[typej]; - if (ni < 0) { + if (ni == 0) { force_buck = r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2; @@ -546,7 +546,7 @@ void PairBuckCoul::compute(int eflag, int vflag) } } else { // cut - if (ni < 0) { + if (ni == 0) { force_buck = r*expr*buck1i[typej]-rn*buck2i[typej]; if (eflag) evdwl = expr*buckai[typej] - rn*buckci[typej]-offseti[typej]; @@ -591,7 +591,6 @@ void PairBuckCoul::compute_inner() int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; double *special_coul = force->special_coul; double *special_lj = force->special_lj; @@ -621,8 +620,9 @@ void PairBuckCoul::compute_inner() jneighn = (jneigh = listinner->firstneigh[i])+listinner->numneigh[i]; for (; jneightype; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; double *special_coul = force->special_coul; double *special_lj = force->special_lj; @@ -711,8 +710,9 @@ void PairBuckCoul::compute_middle() jneighn = (jneigh = listmiddle->firstneigh[i])+listmiddle->numneigh[i]; for (; jneighq; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -812,8 +811,9 @@ void PairBuckCoul::compute_outer(int eflag, int vflag) jneighn = (jneigh = listouter->firstneigh[i])+listouter->numneigh[i]; for (; jneigh> ncoulshiftbits; register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; - if (ni < 0) { + if (ni == 0) { force_coul = qiqj*(ftable[k]+f*dftable[k]); if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); } @@ -870,13 +870,13 @@ void PairBuckCoul::compute_outer(int eflag, int vflag) if (rsq < cut_bucksqi[typej]) { // buckingham register double rn = r2inv*r2inv*r2inv, expr = exp(-r*rhoinvi[typej]); - if (respa_flag) respa_buck = ni<0 ? // correct for respa + if (respa_flag) respa_buck = ni == 0 ? // correct for respa frespa*(r*expr*buck1i[typej]-rn*buck2i[typej]) : frespa*(r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni]; if (order6) { // long-range form register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2)*buckci[typej]; - if (ni < 0) { + if (ni == 0) { force_buck = r*expr*buck1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2; @@ -890,7 +890,7 @@ void PairBuckCoul::compute_outer(int eflag, int vflag) } } else { // cut form - if (ni < 0) { + if (ni == 0) { force_buck = r*expr*buck1i[typej]-rn*buck2i[typej]; if (eflag) evdwl = expr*buckai[typej]-rn*buckci[typej]-offseti[typej]; diff --git a/src/USER-EWALDN/pair_lj_coul.cpp b/src/USER-EWALDN/pair_lj_coul.cpp index 60b18e2773..2fdc0c610d 100644 --- a/src/USER-EWALDN/pair_lj_coul.cpp +++ b/src/USER-EWALDN/pair_lj_coul.cpp @@ -459,7 +459,6 @@ void PairLJCoul::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -484,8 +483,9 @@ void PairLJCoul::compute(int eflag, int vflag) jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i]; for (; jneigh>ncoulshiftbits; register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; - if (ni < 0) { + if (ni == 0) { force_coul = qiqj*(ftable[k]+f*dftable[k]); if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); } @@ -533,7 +533,7 @@ void PairLJCoul::compute(int eflag, int vflag) register double rn = r2inv*r2inv*r2inv; register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2)*lj4i[typej]; - if (ni < 0) { + if (ni == 0) { force_lj = (rn*=rn)*lj1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; if (eflag) @@ -549,7 +549,7 @@ void PairLJCoul::compute(int eflag, int vflag) } else { // cut lj register double rn = r2inv*r2inv*r2inv; - if (ni < 0) { + if (ni == 0) { force_lj = rn*(rn*lj1i[typej]-lj2i[typej]); if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]; } @@ -593,7 +593,6 @@ void PairLJCoul::compute_inner() int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; double *special_coul = force->special_coul; double *special_lj = force->special_lj; @@ -623,8 +622,9 @@ void PairLJCoul::compute_inner() jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i]; for (; jneightype; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *x0 = atom->x[0], *f0 = atom->f[0], *fi = f0, *q = atom->q; double *special_coul = force->special_coul; double *special_lj = force->special_lj; @@ -711,8 +710,9 @@ void PairLJCoul::compute_middle() jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i]; for (; jneighq; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -808,9 +807,10 @@ void PairLJCoul::compute_outer(int eflag, int vflag) jneighn = (jneigh = list->firstneigh[i])+list->numneigh[i]; for (; jneigh> ncoulshiftbits; register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j]; - if (ni < 0) { + if (ni == 0) { force_coul = qiqj*(ftable[k]+f*dftable[k]); if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]); } @@ -864,13 +864,13 @@ void PairLJCoul::compute_outer(int eflag, int vflag) if (rsq < cut_ljsqi[typej]) { // lennard-jones register double rn = r2inv*r2inv*r2inv; - if (respa_flag) respa_lj = ni<0 ? // correct for respa + if (respa_flag) respa_lj = ni == 0 ? // correct for respa frespa*rn*(rn*lj1i[typej]-lj2i[typej]) : frespa*rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni]; if (order6) { // long-range form register double x2 = g2*rsq, a2 = 1.0/x2; x2 = a2*exp(-x2)*lj4i[typej]; - if (ni < 0) { + if (ni == 0) { force_lj = (rn*=rn)*lj1i[typej]-g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq; if (eflag) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2; @@ -884,7 +884,7 @@ void PairLJCoul::compute_outer(int eflag, int vflag) } } else { // cut form - if (ni < 0) { + if (ni == 0) { force_lj = rn*(rn*lj1i[typej]-lj2i[typej]); if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]; } diff --git a/src/USER-REAXC/pair_reax_c.cpp b/src/USER-REAXC/pair_reax_c.cpp index efbb368f9b..a340f9878e 100644 --- a/src/USER-REAXC/pair_reax_c.cpp +++ b/src/USER-REAXC/pair_reax_c.cpp @@ -561,6 +561,7 @@ int PairReaxC::estimate_reax_lists() for( itr_j = 0; itr_j < numneigh[i]; ++itr_j ){ j = jlist[itr_j]; + j &= NEIGHMASK; get_distance( x[j], x[i], &d_sqr, &dvec ); dist[j] = sqrt(d_sqr); @@ -571,6 +572,7 @@ int PairReaxC::estimate_reax_lists() // compute the nbrs among ghost atoms for( itr_j = 0; itr_j < numneigh[i]; ++itr_j ){ j = jlist[itr_j]; + j &= NEIGHMASK; if( j >= nlocal && !marked[j] && dist[j] <= (control->vlist_cut - control->bond_cut) ){ @@ -579,6 +581,7 @@ int PairReaxC::estimate_reax_lists() for( itr_g = 0; itr_g < numneigh[i]; ++itr_g ){ g = jlist[itr_g]; + g &= NEIGHMASK; if( g >= nlocal && !marked[g] ){ get_distance( x[g], x[j], &g_d_sqr, &g_dvec ); @@ -649,6 +652,7 @@ int PairReaxC::write_reax_lists() for( itr_j = 0; itr_j < numneigh[i]; ++itr_j ){ j = jlist[itr_j]; + j &= NEIGHMASK; get_distance( x[j], x[i], &d_sqr, &dvec ); dist[j] = sqrt( d_sqr ); @@ -662,6 +666,7 @@ int PairReaxC::write_reax_lists() // compute the nbrs among ghost atoms for( itr_j = 0; itr_j < numneigh[i]; ++itr_j ){ j = jlist[itr_j]; + j &= NEIGHMASK; if( j >= nlocal && !marked[j] && dist[j] <= (control->vlist_cut - control->bond_cut) ){ @@ -670,6 +675,7 @@ int PairReaxC::write_reax_lists() for( itr_g = 0; itr_g < numneigh[i]; ++itr_g ){ g = jlist[itr_g]; + g &= NEIGHMASK; if( g >= nlocal && !marked[g] ){ get_distance( x[g], x[j], &g_d_sqr, &g_dvec ); diff --git a/src/compute.h b/src/compute.h index ceacee82ac..9b929b345a 100644 --- a/src/compute.h +++ b/src/compute.h @@ -121,6 +121,10 @@ class Compute : protected Pointers { int *molmap; // convert molecule ID to local index int molecules_in_group(int &, int &); + + inline int sbmask(int j) { + return j >> SBBITS & 3; + } }; } diff --git a/src/compute_centro_atom.cpp b/src/compute_centro_atom.cpp index 0b3ac52a51..811d5761a6 100644 --- a/src/compute_centro_atom.cpp +++ b/src/compute_centro_atom.cpp @@ -163,7 +163,7 @@ void ComputeCentroAtom::compute_peratom() n = 0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - if (j >= nall) j %= nall; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/compute_cluster_atom.cpp b/src/compute_cluster_atom.cpp index e426047e6b..38a63d078c 100644 --- a/src/compute_cluster_atom.cpp +++ b/src/compute_cluster_atom.cpp @@ -162,7 +162,7 @@ void ComputeClusterAtom::compute_peratom() n = 0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - if (j >= nall) j %= nall; + j &= NEIGHMASK; if (!(mask[j] & groupbit)) continue; if (clusterID[i] == clusterID[j]) continue; diff --git a/src/compute_cna_atom.cpp b/src/compute_cna_atom.cpp index 6ad00a2c47..f07ffd6fc9 100644 --- a/src/compute_cna_atom.cpp +++ b/src/compute_cna_atom.cpp @@ -168,6 +168,8 @@ void ComputeCNAAtom::compute_peratom() n = 0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; + delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; @@ -247,6 +249,7 @@ void ComputeCNAAtom::compute_peratom() n = 0; for (kk = 0; kk < jnum; kk++) { k = jlist[kk]; + k &= NEIGHMASK; if (k == j) continue; delx = xtmp - x[k][0]; diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp index 46280b69b0..1461be86f8 100644 --- a/src/compute_coord_atom.cpp +++ b/src/compute_coord_atom.cpp @@ -132,7 +132,7 @@ void ComputeCoordAtom::compute_peratom() n = 0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - if (j >= nall) j %= nall; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/compute_group_group.cpp b/src/compute_group_group.cpp index bf1cb81846..e7f80f66b7 100644 --- a/src/compute_group_group.cpp +++ b/src/compute_group_group.cpp @@ -165,13 +165,9 @@ void ComputeGroupGroup::interact() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; if (!(mask[j] & othergroupbit)) continue; diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index 71828e5043..a8b6bff076 100644 --- a/src/compute_pair_local.cpp +++ b/src/compute_pair_local.cpp @@ -169,13 +169,9 @@ int ComputePairLocal::compute_pairs(int flag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; if (!(mask[j] & groupbit)) continue; if (newton_pair == 0 && j >= nlocal) continue; diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index 4f5ef51a3c..8955393a14 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -322,13 +322,9 @@ int ComputePropertyLocal::count_pairs(int allflag, int forceflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; if (!(mask[j] & groupbit)) continue; if (newton_pair == 0 && j >= nlocal) continue; diff --git a/src/compute_rdf.cpp b/src/compute_rdf.cpp index 04848b3923..e1f75014e9 100644 --- a/src/compute_rdf.cpp +++ b/src/compute_rdf.cpp @@ -178,6 +178,7 @@ void ComputeRDF::compute_array() int i,j,m,ii,jj,inum,jnum,itype,jtype,ipair,jpair,ibin,ihisto; double xtmp,ytmp,ztmp,delx,dely,delz,r; int *ilist,*jlist,*numneigh,**firstneigh; + double factor_lj,factor_coul; invoked_array = update->ntimestep; @@ -199,8 +200,6 @@ void ComputeRDF::compute_array() // tally the RDF // both atom i and j must be in fix group // itype,jtype must have been specified by user - // weighting factor must be != 0.0 for this pair - // could be 0 and still be in neigh list for long-range Coulombics // consider I,J as one interaction even if neighbor pair is stored on 2 procs // tally I,J pair each time I is central atom, and each time J is central @@ -226,12 +225,15 @@ void ComputeRDF::compute_array() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; - if (j >= nall) { - if (special_coul[j/nall] == 0.0 && special_lj[j/nall] == 0.0) - continue; - j %= nall; - } + // if both weighting factors are 0, skip this pair + // could be 0 and still be in neigh list for long-range Coulombics + // want consistency with non-charged pairs which wouldn't be in list + + if (factor_lj == 0.0 && factor_coul == 0.0) continue; if (!(mask[j] & groupbit)) continue; jtype = type[j]; diff --git a/src/delete_atoms.cpp b/src/delete_atoms.cpp index e450828c6c..56879053d6 100644 --- a/src/delete_atoms.cpp +++ b/src/delete_atoms.cpp @@ -251,6 +251,7 @@ void DeleteAtoms::delete_overlap(int narg, char **arg) int i,j,ii,jj,inum,jnum; double xtmp,ytmp,ztmp,delx,dely,delz,rsq; int *ilist,*jlist,*numneigh,**firstneigh; + double factor_lj,factor_coul; inum = list->inum; ilist = list->ilist; @@ -267,15 +268,15 @@ void DeleteAtoms::delete_overlap(int narg, char **arg) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; - // if weighting factors are 0, skip this pair + // if both weighting factors are 0, skip this pair // could be 0 and still be in neigh list for long-range Coulombics // want consistency with non-charged pairs which wouldn't be in list - if (j >= nall) { - if (special_coul[j/nall] == 0.0 && special_lj[j/nall] == 0.0) continue; - j %= nall; - } + if (factor_lj == 0.0 && factor_coul == 0.0) continue; // only consider deletion if I,J distance < cutoff diff --git a/src/delete_atoms.h b/src/delete_atoms.h index e941728fe2..840c0d65dc 100644 --- a/src/delete_atoms.h +++ b/src/delete_atoms.h @@ -38,6 +38,10 @@ class DeleteAtoms : protected Pointers { void delete_overlap(int, char **); void delete_porosity(int, char **); void options(int, char **); + + inline int sbmask(int j) { + return j >> SBBITS & 3; + } }; } diff --git a/src/fix_orient_fcc.cpp b/src/fix_orient_fcc.cpp index dfef553b0f..6630616940 100644 --- a/src/fix_orient_fcc.cpp +++ b/src/fix_orient_fcc.cpp @@ -289,6 +289,7 @@ void FixOrientFCC::post_force(int vflag) nsort = 0; for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; count++; dx = x[i][0] - x[j][0]; diff --git a/src/fix_shear_history.cpp b/src/fix_shear_history.cpp index 328756e49e..0740a2f257 100644 --- a/src/fix_shear_history.cpp +++ b/src/fix_shear_history.cpp @@ -136,6 +136,7 @@ void FixShearHistory::pre_exchange() if (touch[jj]) { shear = &allshear[3*jj]; j = jlist[jj]; + j &= NEIGHMASK; if (npartner[i] < MAXTOUCH) { m = npartner[i]; partner[i][m] = tag[j]; diff --git a/src/force.cpp b/src/force.cpp index a53ac9e25e..23a4707d5a 100644 --- a/src/force.cpp +++ b/src/force.cpp @@ -47,6 +47,7 @@ Force::Force(LAMMPS *lmp) : Pointers(lmp) { newton = newton_pair = newton_bond = 1; + special_lj[0] = special_coul[0] = 1.0; special_lj[1] = special_lj[2] = special_lj[3] = 0.0; special_coul[1] = special_coul[2] = special_coul[3] = 0.0; special_angle = special_dihedral = 0; diff --git a/src/lmptype.h b/src/lmptype.h index a1f6662fe6..3ce417c91e 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -41,6 +41,12 @@ namespace LAMMPS_NS { +// reserve 2 hi bits in molecular system neigh list for special bonds flag +// max local + ghost atoms per processor = 2^30 - 1 + +#define SBBITS 30 +#define NEIGHMASK 0x3FFFFFFF + // default settings // 32-bit smallint and tagint, 64-bit bigint diff --git a/src/neigh_derive.cpp b/src/neigh_derive.cpp index c33fac9e62..c3d3eb0245 100644 --- a/src/neigh_derive.cpp +++ b/src/neigh_derive.cpp @@ -27,7 +27,7 @@ using namespace LAMMPS_NS; void Neighbor::half_from_full_no_newton(NeighList *list) { - int i,j,ii,jj,n,jnum; + int i,j,ii,jj,n,jnum,joriginal; int *neighptr,*jlist; int *ilist = list->ilist; @@ -63,8 +63,9 @@ void Neighbor::half_from_full_no_newton(NeighList *list) jnum = numneigh_full[i]; for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - if (j > i) neighptr[n++] = j; + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; + if (j > i) neighptr[n++] = joriginal; } ilist[inum++] = i; @@ -93,7 +94,6 @@ void Neighbor::half_from_full_newton(NeighList *list) double **x = atom->x; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -132,11 +132,11 @@ void Neighbor::half_from_full_newton(NeighList *list) jnum = numneigh_full[i]; for (jj = 0; jj < jnum; jj++) { - j = joriginal = jlist[jj]; + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; if (j < nlocal) { if (i > j) continue; } else { - if (j >= nall) j %= nall; if (x[j][2] < ztmp) continue; if (x[j][2] == ztmp) { if (x[j][1] < ytmp) continue; @@ -171,7 +171,6 @@ void Neighbor::skip_from(NeighList *list) int *type = atom->type; int nlocal = atom->nlocal; - int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -214,8 +213,8 @@ void Neighbor::skip_from(NeighList *list) jnum = numneigh_skip[i]; for (jj = 0; jj < jnum; jj++) { - j = joriginal = jlist[jj]; - if (j >= nall) j %= nall; + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; if (ijskip[itype][type[j]]) continue; neighptr[n++] = joriginal; } @@ -252,7 +251,6 @@ void Neighbor::skip_from_granular(NeighList *list) double *shearptr,*shearptr_skip; int *type = atom->type; - int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -311,8 +309,8 @@ void Neighbor::skip_from_granular(NeighList *list) jnum = numneigh_skip[i]; for (jj = 0; jj < jnum; jj++) { - j = joriginal = jlist[jj]; - if (j >= nall) j %= nall; + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; if (ijskip[itype][type[j]]) continue; neighptr[n] = joriginal; touchptr[n++] = touchptr_skip[jj]; @@ -346,7 +344,6 @@ void Neighbor::skip_from_respa(NeighList *list) int *neighptr,*jlist,*neighptr_inner,*neighptr_middle; int *type = atom->type; - int nall = atom->nlocal + atom->nghost; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -431,8 +428,8 @@ void Neighbor::skip_from_respa(NeighList *list) jnum = numneigh_skip[i]; for (jj = 0; jj < jnum; jj++) { - j = joriginal = jlist[jj]; - if (j >= nall) j %= nall; + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; if (ijskip[itype][type[j]]) continue; neighptr[n++] = joriginal; } @@ -443,8 +440,8 @@ void Neighbor::skip_from_respa(NeighList *list) jnum = numneigh_inner_skip[i]; for (jj = 0; jj < jnum; jj++) { - j = joriginal = jlist[jj]; - if (j >= nall) j %= nall; + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; if (ijskip[itype][type[j]]) continue; neighptr_inner[n_inner++] = joriginal; } @@ -456,8 +453,8 @@ void Neighbor::skip_from_respa(NeighList *list) jnum = numneigh_middle_skip[i]; for (jj = 0; jj < jnum; jj++) { - j = joriginal = jlist[jj]; - if (j >= nall) j %= nall; + joriginal = jlist[jj]; + j = joriginal & NEIGHMASK; if (ijskip[itype][type[j]]) continue; neighptr_middle[n_middle++] = joriginal; } diff --git a/src/neigh_full.cpp b/src/neigh_full.cpp index a173fefeb2..c667ab6160 100644 --- a/src/neigh_full.cpp +++ b/src/neigh_full.cpp @@ -87,10 +87,10 @@ void Neighbor::full_nsq(NeighList *list) delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq <= cutneighsq[itype][jtype]) { - if (molecular) which = find_special(special[i],nspecial[i],tag[j]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } @@ -171,10 +171,10 @@ void Neighbor::full_nsq_ghost(NeighList *list) delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq <= cutneighsq[itype][jtype]) { - if (molecular) which = find_special(special[i],nspecial[i],tag[j]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } else { @@ -188,10 +188,10 @@ void Neighbor::full_nsq_ghost(NeighList *list) delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; if (rsq <= cutneighghostsq[itype][jtype]) { - if (molecular) which = find_special(special[i],nspecial[i],tag[j]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } @@ -232,7 +232,6 @@ void Neighbor::full_bin(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int molecular = atom->molecular; if (includegroup) nlocal = atom->nfirst; @@ -283,10 +282,10 @@ void Neighbor::full_bin(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } @@ -381,10 +380,10 @@ void Neighbor::full_bin_ghost(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } @@ -410,10 +409,10 @@ void Neighbor::full_bin_ghost(NeighList *list) rsq = delx*delx + dely*dely + delz*delz; if (rsq <= cutneighghostsq[itype][jtype]) { - if (molecular) which = find_special(special[i],nspecial[i],tag[j]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } @@ -459,7 +458,6 @@ void Neighbor::full_multi(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int molecular = atom->molecular; if (includegroup) nlocal = atom->nfirst; @@ -514,10 +512,10 @@ void Neighbor::full_multi(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } diff --git a/src/neigh_half_bin.cpp b/src/neigh_half_bin.cpp index 3432b4a0d1..e776b82030 100644 --- a/src/neigh_half_bin.cpp +++ b/src/neigh_half_bin.cpp @@ -46,7 +46,6 @@ void Neighbor::half_bin_no_newton(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; if (includegroup) nlocal = atom->nfirst; int molecular = atom->molecular; @@ -97,10 +96,10 @@ void Neighbor::half_bin_no_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } @@ -143,9 +142,8 @@ void Neighbor::half_bin_newton(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; - int molecular = atom->molecular; if (includegroup) nlocal = atom->nfirst; + int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -196,10 +194,10 @@ void Neighbor::half_bin_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } @@ -217,10 +215,10 @@ void Neighbor::half_bin_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } @@ -263,9 +261,8 @@ void Neighbor::half_bin_newton_tri(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; - int molecular = atom->molecular; if (includegroup) nlocal = atom->nfirst; + int molecular = atom->molecular; int *ilist = list->ilist; int *numneigh = list->numneigh; @@ -321,10 +318,10 @@ void Neighbor::half_bin_newton_tri(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } diff --git a/src/neigh_half_multi.cpp b/src/neigh_half_multi.cpp index b61a53114d..d0a8efcd2b 100644 --- a/src/neigh_half_multi.cpp +++ b/src/neigh_half_multi.cpp @@ -48,7 +48,6 @@ void Neighbor::half_multi_no_newton(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int molecular = atom->molecular; if (includegroup) nlocal = atom->nfirst; @@ -105,10 +104,10 @@ void Neighbor::half_multi_no_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } @@ -153,7 +152,6 @@ void Neighbor::half_multi_newton(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int molecular = atom->molecular; if (includegroup) nlocal = atom->nfirst; @@ -207,10 +205,10 @@ void Neighbor::half_multi_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } @@ -235,10 +233,10 @@ void Neighbor::half_multi_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } @@ -283,7 +281,6 @@ void Neighbor::half_multi_newton_tri(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int molecular = atom->molecular; if (includegroup) nlocal = atom->nfirst; @@ -349,10 +346,10 @@ void Neighbor::half_multi_newton_tri(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } } diff --git a/src/neigh_half_nsq.cpp b/src/neigh_half_nsq.cpp index 0447c3ffd4..1df4865984 100644 --- a/src/neigh_half_nsq.cpp +++ b/src/neigh_half_nsq.cpp @@ -87,10 +87,10 @@ void Neighbor::half_nsq_no_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } @@ -191,10 +191,10 @@ void Neighbor::half_nsq_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; } } diff --git a/src/neigh_respa.cpp b/src/neigh_respa.cpp index a1e0194a33..a5190316ff 100644 --- a/src/neigh_respa.cpp +++ b/src/neigh_respa.cpp @@ -128,19 +128,20 @@ void Neighbor::respa_nsq_no_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; - else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; + else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; - else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); } } } @@ -303,20 +304,21 @@ void Neighbor::respa_nsq_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; - else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; + else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; - else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); } } } @@ -381,7 +383,6 @@ void Neighbor::respa_bin_no_newton(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int molecular = atom->molecular; if (includegroup) nlocal = atom->nfirst; @@ -471,20 +472,22 @@ void Neighbor::respa_bin_no_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; - else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; + else if (which > 0) + neighptr_inner[n_inner++] = j ^ (which << SBBITS); } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; - else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); } } } @@ -549,7 +552,6 @@ void Neighbor::respa_bin_newton(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int molecular = atom->molecular; if (includegroup) nlocal = atom->nfirst; @@ -642,20 +644,21 @@ void Neighbor::respa_bin_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; - else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; + else if (which > 0) neighptr_inner[n_inner++] = j ^ (which << SBBITS); } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; - else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); } } } @@ -674,20 +677,22 @@ void Neighbor::respa_bin_newton(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; - else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; + else if (which > 0) + neighptr_inner[n_inner++] = j ^ (which << SBBITS); } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; - else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); } } } @@ -752,7 +757,6 @@ void Neighbor::respa_bin_newton_tri(NeighList *list) int *mask = atom->mask; int *molecule = atom->molecule; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int molecular = atom->molecular; if (includegroup) nlocal = atom->nfirst; @@ -850,20 +854,22 @@ void Neighbor::respa_bin_newton_tri(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]); - else which = 0; - if (which == 0) neighptr[n++] = j; - else if (which > 0) neighptr[n++] = which*nall + j; + if (molecular) { + which = find_special(special[i],nspecial[i],tag[j]); + if (which >= 0) neighptr[n++] = j ^ (which << SBBITS); + } else neighptr[n++] = j; if (rsq < cut_inner_sq) { if (which == 0) neighptr_inner[n_inner++] = j; - else if (which > 0) neighptr_inner[n_inner++] = which*nall + j; + else if (which > 0) + neighptr_inner[n_inner++] = j ^ (which << SBBITS); } if (respamiddle && rsq < cut_middle_sq && rsq > cut_middle_inside_sq) { if (which == 0) neighptr_middle[n_middle++] = j; - else if (which > 0) neighptr_middle[n_middle++] = which*nall + j; + else if (which > 0) + neighptr_middle[n_middle++] = j ^ (which << SBBITS); } } } diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 86cad4b069..fc4f1125e5 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -1172,13 +1172,10 @@ void Neighbor::build() memory->create(bins,maxbin,"bins"); } - // check that pairwise lists with special bond weighting will not overflow + // check that neighbor list with special bond flags will not overflow - if (atom->molecular && maxwt && nblist) { - bigint max = maxwt * static_cast (atom->nlocal + atom->nghost); - if (max > MAXSMALLINT) - error->one("Weighted neighbor list values are too big"); - } + if (atom->nlocal+atom->nghost > NEIGHMASK) + error->one("Too many local+ghost atoms for neighbor list"); // invoke building of pair and molecular neighbor lists // only for pairwise lists with buildflag set @@ -1218,6 +1215,11 @@ void Neighbor::build_one(int i) memory->create(bins,maxbin,"bins"); } + // check that neighbor list with special bond flags will not overflow + + if (atom->nlocal+atom->nghost > NEIGHMASK) + error->one("Too many local+ghost atoms for neighbor list"); + // when occasional list built, LAMMPS can crash if atoms have moved too far // why is this?, give warning if this is the case // no easy workaround b/c all neighbor lists really need to be rebuilt diff --git a/src/neighbor.h b/src/neighbor.h index 64d12788a6..97e0f354c4 100644 --- a/src/neighbor.h +++ b/src/neighbor.h @@ -165,39 +165,6 @@ class Neighbor : protected Pointers { void choose_build(int, class NeighRequest *); void choose_stencil(int, class NeighRequest *); - // find_special: determine if atom j is in special list of atom i - // if it is not, return 0 - // if it is and special flag is 0 (both coeffs are 0.0), return -1 - // if it is and special flag is 1 (both coeffs are 1.0), return 0 - // if it is and special flag is 2 (otherwise), return 1,2,3 - // for which neighbor it is (and which coeff it maps to) - - inline int find_special(const int *list, const int *nspecial, - const int tag) const { - const int n1 = nspecial[0]; - const int n2 = nspecial[1]; - const int n3 = nspecial[2]; - - for (int i = 0; i < n3; i++) { - if (list[i] == tag) { - if (i < n1) { - if (special_flag[1] == 0) return -1; - else if (special_flag[1] == 1) return 0; - else return 1; - } else if (i < n2) { - if (special_flag[2] == 0) return -1; - else if (special_flag[2] == 1) return 0; - else return 2; - } else { - if (special_flag[3] == 0) return -1; - else if (special_flag[3] == 1) return 0; - else return 3; - } - } - } - return 0; - }; - // pairwise build functions typedef void (Neighbor::*PairPtr)(class NeighList *); @@ -284,6 +251,39 @@ class Neighbor : protected Pointers { BondPtr improper_build; // ptr to improper list functions void improper_all(); // improper list with all impropers void improper_partial(); // exclude certain impropers + + // find_special: determine if atom j is in special list of atom i + // if it is not, return 0 + // if it is and special flag is 0 (both coeffs are 0.0), return -1 + // if it is and special flag is 1 (both coeffs are 1.0), return 0 + // if it is and special flag is 2 (otherwise), return 1,2,3 + // for which neighbor it is (and which coeff it maps to) + + inline int find_special(const int *list, const int *nspecial, + const int tag) const { + const int n1 = nspecial[0]; + const int n2 = nspecial[1]; + const int n3 = nspecial[2]; + + for (int i = 0; i < n3; i++) { + if (list[i] == tag) { + if (i < n1) { + if (special_flag[1] == 0) return -1; + else if (special_flag[1] == 1) return 0; + else return 1; + } else if (i < n2) { + if (special_flag[2] == 0) return -1; + else if (special_flag[2] == 1) return 0; + else return 2; + } else { + if (special_flag[3] == 0) return -1; + else if (special_flag[3] == 1) return 0; + else return 3; + } + } + } + return 0; + }; }; } diff --git a/src/pair.h b/src/pair.h index 871c0688e4..1ce8fea1a6 100644 --- a/src/pair.h +++ b/src/pair.h @@ -147,6 +147,10 @@ class Pair : protected Pointers { void v_tally_tensor(int, int, int, int, double, double, double, double, double, double); void virial_compute(); + + inline int sbmask(int j) { + return j >> SBBITS & 3; + } }; } diff --git a/src/pair_born.cpp b/src/pair_born.cpp index ba9953d8b1..99100d5c23 100644 --- a/src/pair_born.cpp +++ b/src/pair_born.cpp @@ -76,7 +76,6 @@ void PairBorn::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -98,12 +97,8 @@ void PairBorn::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_buck.cpp b/src/pair_buck.cpp index 6b91e23056..cc2a376114 100644 --- a/src/pair_buck.cpp +++ b/src/pair_buck.cpp @@ -69,7 +69,6 @@ void PairBuck::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -91,12 +90,8 @@ void PairBuck::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_buck_coul_cut.cpp b/src/pair_buck_coul_cut.cpp index 9fa00e13b6..f34956afd1 100644 --- a/src/pair_buck_coul_cut.cpp +++ b/src/pair_buck_coul_cut.cpp @@ -77,7 +77,6 @@ void PairBuckCoulCut::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -102,13 +101,9 @@ void PairBuckCoulCut::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_coul_cut.cpp b/src/pair_coul_cut.cpp index cc92998f02..130b9efd2e 100644 --- a/src/pair_coul_cut.cpp +++ b/src/pair_coul_cut.cpp @@ -64,7 +64,6 @@ void PairCoulCut::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; @@ -88,12 +87,8 @@ void PairCoulCut::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = 1.0; - else { - factor_coul = special_coul[j/nall]; - j %= nall; - } + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_coul_debye.cpp b/src/pair_coul_debye.cpp index 7e3b60ee19..ca3d3b09af 100644 --- a/src/pair_coul_debye.cpp +++ b/src/pair_coul_debye.cpp @@ -51,7 +51,6 @@ void PairCoulDebye::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; int newton_pair = force->newton_pair; double qqrd2e = force->qqrd2e; @@ -75,12 +74,8 @@ void PairCoulDebye::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = 1.0; - else { - factor_coul = special_coul[j/nall]; - j %= nall; - } + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_dpd.cpp b/src/pair_dpd.cpp index 6a42decede..f3db674f85 100644 --- a/src/pair_dpd.cpp +++ b/src/pair_dpd.cpp @@ -80,7 +80,6 @@ void PairDPD::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double dtinvsqrt = 1.0/sqrt(update->dt); @@ -106,12 +105,8 @@ void PairDPD::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_dpd = 1.0; - else { - factor_dpd = special_lj[j/nall]; - j %= nall; - } + factor_dpd = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_dpd_tstat.cpp b/src/pair_dpd_tstat.cpp index 3a22e08eaf..b64ff6b0f2 100644 --- a/src/pair_dpd_tstat.cpp +++ b/src/pair_dpd_tstat.cpp @@ -65,7 +65,6 @@ void PairDPDTstat::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; double dtinvsqrt = 1.0/sqrt(update->dt); @@ -91,12 +90,8 @@ void PairDPDTstat::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_dpd = 1.0; - else { - factor_dpd = special_lj[j/nall]; - j %= nall; - } + factor_dpd = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_gauss.cpp b/src/pair_gauss.cpp index ad3d15f4ae..cd76d6cf8e 100644 --- a/src/pair_gauss.cpp +++ b/src/pair_gauss.cpp @@ -78,7 +78,6 @@ void PairGauss::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; int newton_pair = force->newton_pair; inum = list->inum; @@ -99,6 +98,7 @@ void PairGauss::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_lj96_cut.cpp b/src/pair_lj96_cut.cpp index 03b9a4455b..4e20798ec8 100644 --- a/src/pair_lj96_cut.cpp +++ b/src/pair_lj96_cut.cpp @@ -80,11 +80,9 @@ void PairLJ96Cut::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; - inum = list->inum; ilist = list->ilist; numneigh = list->numneigh; @@ -103,12 +101,8 @@ void PairLJ96Cut::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -160,7 +154,6 @@ void PairLJ96Cut::compute_inner() double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -189,12 +182,8 @@ void PairLJ96Cut::compute_inner() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -239,7 +228,6 @@ void PairLJ96Cut::compute_middle() double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -273,12 +261,8 @@ void PairLJ96Cut::compute_middle() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -331,7 +315,6 @@ void PairLJ96Cut::compute_outer(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -360,12 +343,8 @@ void PairLJ96Cut::compute_outer(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index 2c36c5401e..24ba2dad2c 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -80,7 +80,6 @@ void PairLJCut::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -102,12 +101,8 @@ void PairLJCut::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -158,7 +153,6 @@ void PairLJCut::compute_inner() double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -187,12 +181,8 @@ void PairLJCut::compute_inner() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -236,7 +226,6 @@ void PairLJCut::compute_middle() double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -270,12 +259,8 @@ void PairLJCut::compute_middle() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; @@ -327,7 +312,6 @@ void PairLJCut::compute_outer(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -356,12 +340,8 @@ void PairLJCut::compute_outer(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_lj_cut_coul_cut.cpp b/src/pair_lj_cut_coul_cut.cpp index c55c50ea6a..1ebf1e2e2e 100644 --- a/src/pair_lj_cut_coul_cut.cpp +++ b/src/pair_lj_cut_coul_cut.cpp @@ -73,7 +73,6 @@ void PairLJCutCoulCut::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -98,13 +97,9 @@ void PairLJCutCoulCut::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_lj_cut_coul_debye.cpp b/src/pair_lj_cut_coul_debye.cpp index e4fb8ce63e..f54332a42d 100644 --- a/src/pair_lj_cut_coul_debye.cpp +++ b/src/pair_lj_cut_coul_debye.cpp @@ -45,7 +45,6 @@ void PairLJCutCoulDebye::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -70,13 +69,9 @@ void PairLJCutCoulDebye::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_lj_expand.cpp b/src/pair_lj_expand.cpp index 7f255258d1..3b80e5771c 100644 --- a/src/pair_lj_expand.cpp +++ b/src/pair_lj_expand.cpp @@ -69,7 +69,6 @@ void PairLJExpand::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -91,12 +90,8 @@ void PairLJExpand::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_lj_gromacs.cpp b/src/pair_lj_gromacs.cpp index b2c567d949..5f0584b513 100644 --- a/src/pair_lj_gromacs.cpp +++ b/src/pair_lj_gromacs.cpp @@ -80,7 +80,6 @@ void PairLJGromacs::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -102,12 +101,8 @@ void PairLJGromacs::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_lj_gromacs_coul_gromacs.cpp b/src/pair_lj_gromacs_coul_gromacs.cpp index 83332d9dc5..23abf053ed 100644 --- a/src/pair_lj_gromacs_coul_gromacs.cpp +++ b/src/pair_lj_gromacs_coul_gromacs.cpp @@ -78,7 +78,6 @@ void PairLJGromacsCoulGromacs::compute(int eflag, int vflag) double *q = atom->q; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -103,13 +102,9 @@ void PairLJGromacsCoulGromacs::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = factor_lj = 1.0; - else { - factor_coul = special_coul[j/nall]; - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_lj_smooth.cpp b/src/pair_lj_smooth.cpp index d49a23bbfb..a880735fed 100644 --- a/src/pair_lj_smooth.cpp +++ b/src/pair_lj_smooth.cpp @@ -79,7 +79,6 @@ void PairLJSmooth::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -101,13 +100,9 @@ void PairLJSmooth::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } - delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; diff --git a/src/pair_morse.cpp b/src/pair_morse.cpp index d01e227fac..efe169fa9a 100644 --- a/src/pair_morse.cpp +++ b/src/pair_morse.cpp @@ -65,7 +65,6 @@ void PairMorse::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -87,12 +86,8 @@ void PairMorse::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_soft.cpp b/src/pair_soft.cpp index c58290e35e..1d4fa8060e 100644 --- a/src/pair_soft.cpp +++ b/src/pair_soft.cpp @@ -66,7 +66,6 @@ void PairSoft::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -88,12 +87,8 @@ void PairSoft::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; diff --git a/src/pair_table.cpp b/src/pair_table.cpp index a77dd6b5a6..30f1c8eb23 100644 --- a/src/pair_table.cpp +++ b/src/pair_table.cpp @@ -86,7 +86,6 @@ void PairTable::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_lj = force->special_lj; int newton_pair = force->newton_pair; @@ -108,13 +107,9 @@ void PairTable::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; + factor_lj = special_lj[sbmask(j)]; + j &= NEIGHMASK; - if (j < nall) factor_lj = 1.0; - else { - factor_lj = special_lj[j/nall]; - j %= nall; - } - delx = xtmp - x[j][0]; dely = ytmp - x[j][1]; delz = ztmp - x[j][2]; diff --git a/src/pair_yukawa.cpp b/src/pair_yukawa.cpp index b755aefb6e..430a9604bf 100644 --- a/src/pair_yukawa.cpp +++ b/src/pair_yukawa.cpp @@ -61,7 +61,6 @@ void PairYukawa::compute(int eflag, int vflag) double **f = atom->f; int *type = atom->type; int nlocal = atom->nlocal; - int nall = nlocal + atom->nghost; double *special_coul = force->special_coul; int newton_pair = force->newton_pair; @@ -83,12 +82,8 @@ void PairYukawa::compute(int eflag, int vflag) for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - - if (j < nall) factor_coul = 1.0; - else { - factor_coul = special_coul[j/nall]; - j %= nall; - } + factor_coul = special_coul[sbmask(j)]; + j &= NEIGHMASK; delx = xtmp - x[j][0]; dely = ytmp - x[j][1];