From cace3e3530cbb6895a222f772fee2679d908614a Mon Sep 17 00:00:00 2001 From: Stefan Paquay Date: Tue, 30 May 2017 16:08:32 -0400 Subject: [PATCH 01/22] Added missing :pre to doc/src/fix_adapt.txt --- doc/src/fix_adapt.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/fix_adapt.txt b/doc/src/fix_adapt.txt index d7c32bef3d..19d1009b8a 100644 --- a/doc/src/fix_adapt.txt +++ b/doc/src/fix_adapt.txt @@ -47,7 +47,7 @@ keyword = {scale} or {reset} :l fix 1 all adapt 1 pair soft a 1 1 v_prefactor fix 1 all adapt 1 pair soft a 2* 3 v_prefactor fix 1 all adapt 1 pair lj/cut epsilon * * v_scale1 coul/cut scale 3 3 v_scale2 scale yes reset yes -fix 1 all adapt 10 atom diameter v_size +fix 1 all adapt 10 atom diameter v_size :pre variable ramp_up equal "ramp(0.01,0.5)" fix stretch all adapt 1 bond harmonic r0 1 v_ramp_up :pre From 03ab8d0f48aa432ac444adc49a1d92567f7010c3 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 30 May 2017 17:04:48 -0400 Subject: [PATCH 02/22] major neighbor list style whitespace cleanup --- src/KOKKOS/npair_kokkos.cpp | 10 +- src/KOKKOS/npair_kokkos.h | 2 +- src/USER-INTEL/nbin_intel.cpp | 14 +- src/USER-INTEL/npair_full_bin_intel.cpp | 222 +++++------ src/USER-INTEL/npair_full_bin_intel.h | 4 +- .../npair_half_bin_newtoff_intel.cpp | 80 ++-- src/USER-INTEL/npair_half_bin_newtoff_intel.h | 2 +- .../npair_half_bin_newton_intel.cpp | 356 +++++++++--------- src/USER-INTEL/npair_half_bin_newton_intel.h | 4 +- .../npair_half_bin_newton_tri_intel.cpp | 178 ++++----- .../npair_half_bin_newton_tri_intel.h | 2 +- src/USER-INTEL/npair_intel.cpp | 2 +- src/npair.cpp | 12 +- src/npair.h | 2 +- src/npair_full_bin.h | 2 +- src/npair_full_bin_atomonly.h | 2 +- src/npair_full_bin_ghost.h | 2 +- src/npair_full_nsq_ghost.h | 2 +- src/npair_half_bin_atomonly_newton.cpp | 2 +- src/npair_half_respa_bin_newton_tri.cpp | 2 +- src/npair_half_size_bin_newton.cpp | 2 +- src/npair_half_size_bin_newton_tri.cpp | 2 +- src/npair_skip_respa.h | 4 +- src/npair_skip_size.cpp | 2 +- src/npair_skip_size.h | 2 +- src/npair_skip_size_off2on.cpp | 2 +- src/npair_skip_size_off2on.h | 4 +- src/npair_skip_size_off2on_oneside.cpp | 4 +- src/npair_skip_size_off2on_oneside.h | 4 +- src/nstencil_full_bin_2d.cpp | 2 +- src/nstencil_full_bin_2d.h | 2 +- src/nstencil_full_bin_3d.h | 2 +- src/nstencil_full_ghost_bin_2d.h | 2 +- src/nstencil_full_ghost_bin_3d.h | 2 +- src/nstencil_full_multi_2d.h | 2 +- src/nstencil_half_bin_2d_newtoff.cpp | 2 +- src/nstencil_half_bin_2d_newton_tri.cpp | 2 +- src/nstencil_half_bin_3d_newtoff.cpp | 2 +- src/nstencil_half_bin_3d_newton_tri.cpp | 2 +- src/nstencil_half_ghost_bin_2d_newtoff.h | 2 +- 40 files changed, 475 insertions(+), 475 deletions(-) diff --git a/src/KOKKOS/npair_kokkos.cpp b/src/KOKKOS/npair_kokkos.cpp index 3614a82cfe..b7b550369d 100644 --- a/src/KOKKOS/npair_kokkos.cpp +++ b/src/KOKKOS/npair_kokkos.cpp @@ -287,12 +287,12 @@ int NeighborKokkosExecute::exclusion(const int &i,const int &j, if (nex_mol) { for (m = 0; m < nex_mol; m++) - if (ex_mol_intra[m]) { // intra-chain: exclude i-j pair if on same molecule + if (ex_mol_intra[m]) { // intra-chain: exclude i-j pair if on same molecule if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] && - molecule[i] == molecule[j]) return 1; - } else // exclude i-j pair if on different molecules - if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] && - molecule[i] != molecule[j]) return 1; + molecule[i] == molecule[j]) return 1; + } else // exclude i-j pair if on different molecules + if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] && + molecule[i] != molecule[j]) return 1; } return 0; diff --git a/src/KOKKOS/npair_kokkos.h b/src/KOKKOS/npair_kokkos.h index 52cdfe0d53..8e81c57618 100644 --- a/src/KOKKOS/npair_kokkos.h +++ b/src/KOKKOS/npair_kokkos.h @@ -385,7 +385,7 @@ struct NPairKokkosBuildFunctor { } #ifdef KOKKOS_HAVE_CUDA __device__ inline - + void operator() (typename Kokkos::TeamPolicy::member_type dev) const { c.template build_ItemCuda(dev); } diff --git a/src/USER-INTEL/nbin_intel.cpp b/src/USER-INTEL/nbin_intel.cpp index 194b9a5f97..c3335b2c26 100644 --- a/src/USER-INTEL/nbin_intel.cpp +++ b/src/USER-INTEL/nbin_intel.cpp @@ -55,7 +55,7 @@ NBinIntel::~NBinIntel() { nocopy(binhead,bins,_atombin,_binpacked:alloc_if(0) free_if(1)) } #endif -} +} /* ---------------------------------------------------------------------- setup for bin_atoms() @@ -71,7 +71,7 @@ void NBinIntel::bin_atoms_setup(int nall) if (_offload_alloc) { const int * binhead = this->binhead; #pragma offload_transfer target(mic:_cop) \ - nocopy(binhead:alloc_if(0) free_if(1)) + nocopy(binhead:alloc_if(0) free_if(1)) } #endif @@ -99,7 +99,7 @@ void NBinIntel::bin_atoms_setup(int nall) const int * _atombin = this->_atombin; const int * _binpacked = this->_binpacked; #pragma offload_transfer target(mic:_cop) \ - nocopy(bins,_atombin,_binpacked:alloc_if(0) free_if(1)) + nocopy(bins,_atombin,_binpacked:alloc_if(0) free_if(1)) } #endif memory->destroy(bins); @@ -157,10 +157,10 @@ void NBinIntel::bin_atoms(IntelBuffers * buffers) { const flt_t dx = (INTEL_BIGP - bboxhi[0]); const flt_t dy = (INTEL_BIGP - bboxhi[1]); const flt_t dz = (INTEL_BIGP - bboxhi[2]); - if (dx * dx + dy * dy + dz * dz < - static_cast(neighbor->cutneighmaxsq)) + if (dx * dx + dy * dy + dz * dz < + static_cast(neighbor->cutneighmaxsq)) error->one(FLERR, - "Intel package expects no atoms within cutoff of {1e15,1e15,1e15}."); + "Intel package expects no atoms within cutoff of {1e15,1e15,1e15}."); } // ---------- Grow and cast/pack buffers ------------- @@ -181,7 +181,7 @@ void NBinIntel::bin_atoms(IntelBuffers * buffers) { { int ifrom, ito, tid; IP_PRE_omp_range_id_align(ifrom, ito, tid, nall, nthreads, - sizeof(ATOM_T)); + sizeof(ATOM_T)); buffers->thr_pack(ifrom, ito, 0); } _fix->stop_watch(TIME_PACK); diff --git a/src/USER-INTEL/npair_full_bin_intel.cpp b/src/USER-INTEL/npair_full_bin_intel.cpp index 1ec93bf113..7e0d2abdcb 100644 --- a/src/USER-INTEL/npair_full_bin_intel.cpp +++ b/src/USER-INTEL/npair_full_bin_intel.cpp @@ -70,12 +70,12 @@ fbi(NeighList *list, IntelBuffers *buffers) { #endif buffers->grow_list(list, atom->nlocal, comm->nthreads, off_end, - _fix->nbor_pack_width()); + _fix->nbor_pack_width()); int need_ic = 0; if (atom->molecular) dminimum_image_check(need_ic, neighbor->cutneighmax, neighbor->cutneighmax, - neighbor->cutneighmax); + neighbor->cutneighmax); #ifdef _LMP_INTEL_OFFLOAD if (need_ic) { @@ -167,7 +167,7 @@ fbi(const int offload, NeighList *list, IntelBuffers *buffers, overflow = _fix->get_off_overflow_flag(); _fix->stop_watch(TIME_HOST_NEIGHBOR); _fix->start_watch(TIME_OFFLOAD_LATENCY); - } else + } else #endif { tnum = comm->nthreads; @@ -255,8 +255,8 @@ fbi(const int offload, NeighList *list, IntelBuffers *buffers, ito += astart; int e_ito = ito; if (ito == num) { - int imod = ito % pack_width; - if (imod) e_ito += pack_width - imod; + int imod = ito % pack_width; + if (imod) e_ito += pack_width - imod; } const int list_size = (e_ito + tid * 2 + 2) * maxnbors; int which; @@ -276,91 +276,91 @@ fbi(const int offload, NeighList *list, IntelBuffers *buffers, const int ioffset = ntypes * itype; const int ibin = atombin[i]; - int raw_count = pack_offset; + int raw_count = pack_offset; // loop over all atoms in surrounding bins in stencil including self // skip i = j - if (exclude) { - for (int k = 0; k < nstencilp; k++) { - const int bstart = binhead[ibin + binstart[k]]; - const int bend = binhead[ibin + binend[k]]; + if (exclude) { + for (int k = 0; k < nstencilp; k++) { + const int bstart = binhead[ibin + binstart[k]]; + const int bend = binhead[ibin + binend[k]]; #ifndef _LMP_INTEL_OFFLOAD - #ifdef INTEL_VMASK + #ifdef INTEL_VMASK #pragma simd - #endif #endif - for (int jj = bstart; jj < bend; jj++) { - int j = binpacked[jj]; + #endif + for (int jj = bstart; jj < bend; jj++) { + int j = binpacked[jj]; + + if (i == j) j=e_nall; - if (i == j) j=e_nall; - #ifdef _LMP_INTEL_OFFLOAD - if (offload_noghost) { + if (offload_noghost) { if (j < nlocal) { if (i < offload_end) continue; } else if (offload) continue; } - #endif + #endif #ifndef _LMP_INTEL_OFFLOAD - const int jtype = x[j].w; - if (exclusion(i,j,itype,jtype,mask,molecule)) continue; - #endif + const int jtype = x[j].w; + if (exclusion(i,j,itype,jtype,mask,molecule)) continue; + #endif - neighptr[raw_count++] = j; + neighptr[raw_count++] = j; } } - } else { - for (int k = 0; k < nstencilp; k++) { - const int bstart = binhead[ibin + binstart[k]]; - const int bend = binhead[ibin + binend[k]]; + } else { + for (int k = 0; k < nstencilp; k++) { + const int bstart = binhead[ibin + binstart[k]]; + const int bend = binhead[ibin + binend[k]]; #ifndef _LMP_INTEL_OFFLOAD - #ifdef INTEL_VMASK + #ifdef INTEL_VMASK #pragma simd #endif #endif - for (int jj = bstart; jj < bend; jj++) { - int j = binpacked[jj]; + for (int jj = bstart; jj < bend; jj++) { + int j = binpacked[jj]; + + if (i == j) j=e_nall; - if (i == j) j=e_nall; - #ifdef _LMP_INTEL_OFFLOAD - if (offload_noghost) { + if (offload_noghost) { if (j < nlocal) { if (i < offload_end) continue; } else if (offload) continue; } - #endif + #endif - neighptr[raw_count++] = j; + neighptr[raw_count++] = j; } } - } + } - if (raw_count > obound) *overflow = 1; + if (raw_count > obound) *overflow = 1; #if defined(LMP_SIMD_COMPILER) - #ifdef _LMP_INTEL_OFFLOAD - int vlmin = lmin, vlmax = lmax, vgmin = gmin, vgmax = gmax; - #if __INTEL_COMPILER+0 > 1499 - #pragma vector aligned + #ifdef _LMP_INTEL_OFFLOAD + int vlmin = lmin, vlmax = lmax, vgmin = gmin, vgmax = gmax; + #if __INTEL_COMPILER+0 > 1499 + #pragma vector aligned #pragma simd reduction(max:vlmax,vgmax) reduction(min:vlmin, vgmin) - #endif - #else - #pragma vector aligned + #endif + #else + #pragma vector aligned #pragma simd - #endif - #endif - for (int u = pack_offset; u < raw_count; u++) { + #endif + #endif + for (int u = pack_offset; u < raw_count; u++) { int j = neighptr[u]; const flt_t delx = xtmp - x[j].x; const flt_t dely = ytmp - x[j].y; const flt_t delz = ztmp - x[j].z; const int jtype = x[j].w; const flt_t rsq = delx * delx + dely * dely + delz * delz; - if (rsq > cutneighsq[ioffset + jtype]) - neighptr[u] = e_nall; - else { + if (rsq > cutneighsq[ioffset + jtype]) + neighptr[u] = e_nall; + else { if (need_ic) { int no_special; ominimum_image_check(no_special, delx, dely, delz); @@ -376,73 +376,73 @@ fbi(const int offload, NeighList *list, IntelBuffers *buffers, if (j > vgmax) vgmax = j; } #endif - } - } + } + } #ifdef _LMP_INTEL_OFFLOAD - lmin = MIN(lmin,vlmin); - gmin = MIN(gmin,vgmin); - lmax = MAX(lmax,vlmax); - gmax = MAX(gmax,vgmax); + lmin = MIN(lmin,vlmin); + gmin = MIN(gmin,vgmin); + lmax = MAX(lmax,vlmax); + gmax = MAX(gmax,vgmax); #endif int n = lane, n2 = pack_offset; - for (int u = pack_offset; u < raw_count; u++) { - const int j = neighptr[u]; - int pj = j; - if (pj < e_nall) { - if (need_ic) - if (pj < 0) pj = -pj - 1; - - const int jtag = tag[pj]; - int flist = 0; - if (itag > jtag) { - if ((itag+jtag) % 2 == 0) flist = 1; - } else if (itag < jtag) { - if ((itag+jtag) % 2 == 1) flist = 1; - } else { + for (int u = pack_offset; u < raw_count; u++) { + const int j = neighptr[u]; + int pj = j; + if (pj < e_nall) { + if (need_ic) + if (pj < 0) pj = -pj - 1; + + const int jtag = tag[pj]; + int flist = 0; + if (itag > jtag) { + if ((itag+jtag) % 2 == 0) flist = 1; + } else if (itag < jtag) { + if ((itag+jtag) % 2 == 1) flist = 1; + } else { if (x[pj].z < ztmp) flist = 1; - else if (x[pj].z == ztmp && x[pj].y < ytmp) flist = 1; - else if (x[pj].z == ztmp && x[pj].y == ytmp && x[pj].x < xtmp) - flist = 1; - } - if (flist) { - neighptr[n2++] = j; - } else { - neighptr[n] = j; - n += pack_width; - } + else if (x[pj].z == ztmp && x[pj].y < ytmp) flist = 1; + else if (x[pj].z == ztmp && x[pj].y == ytmp && x[pj].x < xtmp) + flist = 1; + } + if (flist) { + neighptr[n2++] = j; + } else { + neighptr[n] = j; + n += pack_width; + } } } - int ns = (n - lane) / pack_width; - atombin[i] = ns; - for (int u = pack_offset; u < n2; u++) { - neighptr[n] = neighptr[u]; - n += pack_width; - } + int ns = (n - lane) / pack_width; + atombin[i] = ns; + for (int u = pack_offset; u < n2; u++) { + neighptr[n] = neighptr[u]; + n += pack_width; + } ilist[i] = i; cnumneigh[i] = ct + lane; - ns += n2 - pack_offset; + ns += n2 - pack_offset; numneigh[i] = ns; - if (ns > max_chunk) max_chunk = ns; - lane++; - if (lane == pack_width) { - ct += max_chunk * pack_width; - const int alignb = (INTEL_DATA_ALIGN / sizeof(int)); - const int edge = (ct % alignb); - if (edge) ct += alignb - edge; - neighptr = firstneigh + ct; - max_chunk = 0; - pack_offset = maxnbors * pack_width; - lane = 0; - if (ct + obound > list_size) { - if (i < ito - 1) { - *overflow = 1; - ct = (ifrom + tid * 2) * maxnbors; - } + if (ns > max_chunk) max_chunk = ns; + lane++; + if (lane == pack_width) { + ct += max_chunk * pack_width; + const int alignb = (INTEL_DATA_ALIGN / sizeof(int)); + const int edge = (ct % alignb); + if (edge) ct += alignb - edge; + neighptr = firstneigh + ct; + max_chunk = 0; + pack_offset = maxnbors * pack_width; + lane = 0; + if (ct + obound > list_size) { + if (i < ito - 1) { + *overflow = 1; + ct = (ifrom + tid * 2) * maxnbors; + } } - } + } } if (*overflow == 1) @@ -482,13 +482,13 @@ fbi(const int offload, NeighList *list, IntelBuffers *buffers, int * _noalias jlist = firstneigh + cnumneigh[i]; const int jnum = numneigh[i]; - const int trip = jnum * pack_width; + const int trip = jnum * pack_width; for (int jj = 0; jj < trip; jj+=pack_width) { const int j = jlist[jj]; - if (need_ic && j < 0) { - which = 0; - jlist[jj] = -j - 1; - } else + if (need_ic && j < 0) { + which = 0; + jlist[jj] = -j - 1; + } else ofind_special(which, special, nspecial, i, tag[j]); #ifdef _LMP_INTEL_OFFLOAD if (j >= nlocal) { @@ -511,9 +511,9 @@ fbi(const int offload, NeighList *list, IntelBuffers *buffers, int jj = 0; for (jj = 0; jj < jnum; jj++) { if (jlist[jj] >= nlocal) { - if (jlist[jj] == e_nall) jlist[jj] = nall_offset; - else jlist[jj] -= ghost_offset; - } + if (jlist[jj] == e_nall) jlist[jj] = nall_offset; + else jlist[jj] -= ghost_offset; + } } } } diff --git a/src/USER-INTEL/npair_full_bin_intel.h b/src/USER-INTEL/npair_full_bin_intel.h index 608bd0f5dd..f1be71abbc 100644 --- a/src/USER-INTEL/npair_full_bin_intel.h +++ b/src/USER-INTEL/npair_full_bin_intel.h @@ -15,7 +15,7 @@ NPairStyle(full/bin/intel, NPairFullBinIntel, - NP_FULL | NP_BIN | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | + NP_FULL | NP_BIN | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI | NP_INTEL) #else @@ -38,7 +38,7 @@ class NPairFullBinIntel : public NPairIntel { void fbi(NeighList *, IntelBuffers *); template void fbi(const int, NeighList *, IntelBuffers *, const int, - const int, const int offload_end = 0); + const int, const int offload_end = 0); }; } diff --git a/src/USER-INTEL/npair_half_bin_newtoff_intel.cpp b/src/USER-INTEL/npair_half_bin_newtoff_intel.cpp index 1fcc3f0759..9a40e2a07c 100644 --- a/src/USER-INTEL/npair_half_bin_newtoff_intel.cpp +++ b/src/USER-INTEL/npair_half_bin_newtoff_intel.cpp @@ -26,7 +26,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfBinNewtoffIntel::NPairHalfBinNewtoffIntel(LAMMPS *lmp) : +NPairHalfBinNewtoffIntel::NPairHalfBinNewtoffIntel(LAMMPS *lmp) : NPairIntel(lmp) {} /* ---------------------------------------------------------------------- @@ -75,7 +75,7 @@ hbnni(NeighList *list, IntelBuffers *buffers) { int need_ic = 0; if (atom->molecular) dminimum_image_check(need_ic, neighbor->cutneighmax, neighbor->cutneighmax, - neighbor->cutneighmax); + neighbor->cutneighmax); #ifdef _LMP_INTEL_OFFLOAD if (need_ic) { @@ -159,7 +159,7 @@ hbnni(const int offload, NeighList *list, IntelBuffers *buffers, overflow = _fix->get_off_overflow_flag(); _fix->stop_watch(TIME_HOST_NEIGHBOR); _fix->start_watch(TIME_OFFLOAD_LATENCY); - } else + } else #endif { tnum = comm->nthreads; @@ -294,13 +294,13 @@ hbnni(const int offload, NeighList *list, IntelBuffers *buffers, ominimum_image_check(no_special, delx, dely, delz); if (no_special) neighptr[n++] = -j - 1; - else + else neighptr[n++] = j; } else neighptr[n++] = j; #ifdef _LMP_INTEL_OFFLOAD - if (j < lmin) lmin = j; - if (j > lmax) lmax = j; + if (j < lmin) lmin = j; + if (j > lmax) lmax = j; #endif } else { if (need_ic) { @@ -308,16 +308,16 @@ hbnni(const int offload, NeighList *list, IntelBuffers *buffers, ominimum_image_check(no_special, delx, dely, delz); if (no_special) neighptr[n2++] = -j - 1; - else + else neighptr[n2++] = j; } else neighptr[n2++] = j; - #ifdef _LMP_INTEL_OFFLOAD - if (j < gmin) gmin = j; - if (j > gmax) gmax = j; + #ifdef _LMP_INTEL_OFFLOAD + if (j < gmin) gmin = j; + if (j > gmax) gmax = j; #endif - } - } + } + } } } ilist[i] = i; @@ -341,13 +341,13 @@ hbnni(const int offload, NeighList *list, IntelBuffers *buffers, neighptr += n; if (ct + n + maxnbors > list_size) { *overflow = 1; - ct = (ifrom + tid) * maxnbors; + ct = (ifrom + tid) * maxnbors; } } if (*overflow == 1) - for (int i = ifrom; i < ito; i++) - numneigh[i] = 0; + for (int i = ifrom; i < ito; i++) + numneigh[i] = 0; #ifdef _LMP_INTEL_OFFLOAD if (separate_buffers) { @@ -370,7 +370,7 @@ hbnni(const int offload, NeighList *list, IntelBuffers *buffers, if (offload) { ghost_offset = overflow[LMP_GHOST_MIN] - overflow[LMP_LOCAL_MAX] - 1; nall_offset = overflow[LMP_LOCAL_MAX] + 1 + nghost; - } else { + } else { ghost_offset = overflow[LMP_GHOST_MIN] - nlocal; nall_offset = nlocal + nghost; } @@ -383,38 +383,38 @@ hbnni(const int offload, NeighList *list, IntelBuffers *buffers, const int jnum = numneigh[i]; for (int jj = 0; jj < jnum; jj++) { const int j = jlist[jj]; - if (need_ic && j < 0) { - which = 0; - jlist[jj] = -j - 1; - } else + if (need_ic && j < 0) { + which = 0; + jlist[jj] = -j - 1; + } else ofind_special(which, special, nspecial, i, tag[j]); #ifdef _LMP_INTEL_OFFLOAD - if (j >= nlocal) { - if (j == nall) - jlist[jj] = nall_offset; - else if (which) - jlist[jj] = (j-ghost_offset) ^ (which << SBBITS); - else jlist[jj]-=ghost_offset; + if (j >= nlocal) { + if (j == nall) + jlist[jj] = nall_offset; + else if (which) + jlist[jj] = (j-ghost_offset) ^ (which << SBBITS); + else jlist[jj]-=ghost_offset; } else #endif - if (which) jlist[jj] = j ^ (which << SBBITS); + if (which) jlist[jj] = j ^ (which << SBBITS); } } } #ifdef _LMP_INTEL_OFFLOAD else if (separate_buffers) { - for (int i = ifrom; i < ito; ++i) { + for (int i = ifrom; i < ito; ++i) { int * _noalias jlist = firstneigh + cnumneigh[i]; const int jnum = numneigh[i]; - int jj = 0; - for (jj = 0; jj < jnum; jj++) - if (jlist[jj] >= nlocal) break; - while (jj < jnum) { - if (jlist[jj] == nall) jlist[jj] = nall_offset; - else jlist[jj] -= ghost_offset; - jj++; - } - } + int jj = 0; + for (jj = 0; jj < jnum; jj++) + if (jlist[jj] >= nlocal) break; + while (jj < jnum) { + if (jlist[jj] == nall) jlist[jj] = nall_offset; + else jlist[jj] -= ghost_offset; + jj++; + } + } } #endif } // end omp @@ -438,9 +438,9 @@ hbnni(const int offload, NeighList *list, IntelBuffers *buffers, _fix->start_watch(TIME_PACK); _fix->set_neighbor_host_sizes(); buffers->pack_sep_from_single(_fix->host_min_local(), - _fix->host_used_local(), - _fix->host_min_ghost(), - _fix->host_used_ghost()); + _fix->host_used_local(), + _fix->host_min_ghost(), + _fix->host_used_ghost()); _fix->stop_watch(TIME_PACK); } } diff --git a/src/USER-INTEL/npair_half_bin_newtoff_intel.h b/src/USER-INTEL/npair_half_bin_newtoff_intel.h index ccb4560909..49482f8b3e 100644 --- a/src/USER-INTEL/npair_half_bin_newtoff_intel.h +++ b/src/USER-INTEL/npair_half_bin_newtoff_intel.h @@ -38,7 +38,7 @@ class NPairHalfBinNewtoffIntel : public NPairIntel { void hbnni(NeighList *, IntelBuffers *); template void hbnni(const int, NeighList *, IntelBuffers *, const int, - const int); + const int); }; } diff --git a/src/USER-INTEL/npair_half_bin_newton_intel.cpp b/src/USER-INTEL/npair_half_bin_newton_intel.cpp index 5584f962e9..6313ab944f 100644 --- a/src/USER-INTEL/npair_half_bin_newton_intel.cpp +++ b/src/USER-INTEL/npair_half_bin_newton_intel.cpp @@ -26,7 +26,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfBinNewtonIntel::NPairHalfBinNewtonIntel(LAMMPS *lmp) : +NPairHalfBinNewtonIntel::NPairHalfBinNewtonIntel(LAMMPS *lmp) : NPairIntel(lmp) {} /* ---------------------------------------------------------------------- @@ -75,7 +75,7 @@ hbni(NeighList *list, IntelBuffers *buffers) { int need_ic = 0; if (atom->molecular) dminimum_image_check(need_ic, neighbor->cutneighmax, neighbor->cutneighmax, - neighbor->cutneighmax); + neighbor->cutneighmax); #ifdef _LMP_INTEL_OFFLOAD if (need_ic) { @@ -96,7 +96,7 @@ hbni(NeighList *list, IntelBuffers *buffers) { } } #else - if (need_ic) + if (need_ic) hbni(0, list, buffers, host_start, nlocal); else hbni(0, list, buffers, host_start, nlocal); @@ -119,7 +119,7 @@ hbni(const int offload, NeighList *list, IntelBuffers *buffers, if (offload) { if (INTEL_MIC_NBOR_PAD > 1) pad = INTEL_MIC_NBOR_PAD * sizeof(float) / sizeof(flt_t); - } else + } else #endif if (INTEL_NBOR_PAD > 1) pad = INTEL_NBOR_PAD * sizeof(float) / sizeof(flt_t); @@ -172,7 +172,7 @@ hbni(const int offload, NeighList *list, IntelBuffers *buffers, overflow = _fix->get_off_overflow_flag(); _fix->stop_watch(TIME_HOST_NEIGHBOR); _fix->start_watch(TIME_OFFLOAD_LATENCY); - } else + } else #endif { tnum = comm->nthreads; @@ -235,8 +235,8 @@ hbni(const int offload, NeighList *list, IntelBuffers *buffers, int end = stencil[k] + 1; for (int kk = k + 1; kk < nstencil; kk++) { if (stencil[kk-1]+1 == stencil[kk]) { - end++; - k++; + end++; + k++; } else break; } binend[nstencilp] = end; @@ -262,8 +262,8 @@ hbni(const int offload, NeighList *list, IntelBuffers *buffers, ito += astart; int e_ito = ito; if (ito == num) { - int imod = ito % swidth; - if (imod) e_ito += swidth - imod; + int imod = ito % swidth; + if (imod) e_ito += swidth - imod; } const int list_size = (e_ito + tid * 2 + 2) * maxnbors; #else @@ -294,118 +294,118 @@ hbni(const int offload, NeighList *list, IntelBuffers *buffers, // if j is owned atom, store it, since j is beyond i in linked list // if j is ghost, only store if j coords are "above/to the right" of i - int raw_count = pack_offset; + int raw_count = pack_offset; for (int j = bins[i]; j >= 0; j = bins[j]) { if (j >= nlocal) { - #ifdef _LMP_INTEL_OFFLOAD + #ifdef _LMP_INTEL_OFFLOAD if (offload_noghost && offload) continue; - #endif + #endif if (x[j].z < ztmp) continue; if (x[j].z == ztmp) { if (x[j].y < ytmp) continue; if (x[j].y == ytmp && x[j].x < xtmp) continue; } - } + } #ifdef _LMP_INTEL_OFFLOAD else if (offload_noghost && i < offload_end) continue; - #endif + #endif #ifndef _LMP_INTEL_OFFLOAD if (exclude) { - const int jtype = x[j].w; - if (exclusion(i,j,itype,jtype,mask,molecule)) continue; - } - #endif + const int jtype = x[j].w; + if (exclusion(i,j,itype,jtype,mask,molecule)) continue; + } + #endif - neighptr[raw_count++] = j; - } + neighptr[raw_count++] = j; + } // loop over all atoms in other bins in stencil, store every pair const int ibin = atombin[i]; - if (exclude) { - for (int k = 0; k < nstencilp; k++) { - const int bstart = binhead[ibin + binstart[k]]; - const int bend = binhead[ibin + binend[k]]; - #ifndef _LMP_INTEL_OFFLOAD - #ifdef INTEL_VMASK - #pragma simd - #endif - #endif - for (int jj = bstart; jj < bend; jj++) { - const int j = binpacked[jj]; - + if (exclude) { + for (int k = 0; k < nstencilp; k++) { + const int bstart = binhead[ibin + binstart[k]]; + const int bend = binhead[ibin + binend[k]]; + #ifndef _LMP_INTEL_OFFLOAD + #ifdef INTEL_VMASK + #pragma simd + #endif + #endif + for (int jj = bstart; jj < bend; jj++) { + const int j = binpacked[jj]; + #ifdef _LMP_INTEL_OFFLOAD if (offload_noghost) { if (j < nlocal) { if (i < offload_end) continue; } else if (offload) continue; } - #endif + #endif #ifndef _LMP_INTEL_OFFLOAD - const int jtype = x[j].w; - if (exclusion(i,j,itype,jtype,mask,molecule)) continue; - #endif + const int jtype = x[j].w; + if (exclusion(i,j,itype,jtype,mask,molecule)) continue; + #endif + + neighptr[raw_count++] = j; + } + } + } else { + for (int k = 0; k < nstencilp; k++) { + const int bstart = binhead[ibin + binstart[k]]; + const int bend = binhead[ibin + binend[k]]; + #ifndef _LMP_INTEL_OFFLOAD + #ifdef INTEL_VMASK + #pragma simd + #endif + #endif + for (int jj = bstart; jj < bend; jj++) { + const int j = binpacked[jj]; - neighptr[raw_count++] = j; - } - } - } else { - for (int k = 0; k < nstencilp; k++) { - const int bstart = binhead[ibin + binstart[k]]; - const int bend = binhead[ibin + binend[k]]; - #ifndef _LMP_INTEL_OFFLOAD - #ifdef INTEL_VMASK - #pragma simd - #endif - #endif - for (int jj = bstart; jj < bend; jj++) { - const int j = binpacked[jj]; - #ifdef _LMP_INTEL_OFFLOAD if (offload_noghost) { if (j < nlocal) { if (i < offload_end) continue; } else if (offload) continue; } - #endif + #endif - neighptr[raw_count++] = j; - } - } - } + neighptr[raw_count++] = j; + } + } + } - if (raw_count > obound) *overflow = 1; + if (raw_count > obound) *overflow = 1; #if defined(LMP_SIMD_COMPILER) - #ifdef _LMP_INTEL_OFFLOAD - int vlmin = lmin, vlmax = lmax, vgmin = gmin, vgmax = gmax; - #if __INTEL_COMPILER+0 > 1499 - #pragma vector aligned + #ifdef _LMP_INTEL_OFFLOAD + int vlmin = lmin, vlmax = lmax, vgmin = gmin, vgmax = gmax; + #if __INTEL_COMPILER+0 > 1499 + #pragma vector aligned #pragma simd reduction(max:vlmax,vgmax) reduction(min:vlmin, vgmin) - #endif - #else - #pragma vector aligned + #endif + #else + #pragma vector aligned #pragma simd - #endif - #endif - for (int u = pack_offset; u < raw_count; u++) { + #endif + #endif + for (int u = pack_offset; u < raw_count; u++) { int j = neighptr[u]; const flt_t delx = xtmp - x[j].x; const flt_t dely = ytmp - x[j].y; const flt_t delz = ztmp - x[j].z; - const int jtype = x[j].w; + const int jtype = x[j].w; const flt_t rsq = delx * delx + dely * dely + delz * delz; - if (rsq > cutneighsq[ioffset + jtype]) - neighptr[u] = e_nall; - else { - if (need_ic) { - int no_special; - ominimum_image_check(no_special, delx, dely, delz); - if (no_special) - neighptr[u] = -j - 1; - } + if (rsq > cutneighsq[ioffset + jtype]) + neighptr[u] = e_nall; + else { + if (need_ic) { + int no_special; + ominimum_image_check(no_special, delx, dely, delz); + if (no_special) + neighptr[u] = -j - 1; + } #ifdef _LMP_INTEL_OFFLOAD if (j < nlocal) { if (j < vlmin) vlmin = j; @@ -415,40 +415,40 @@ hbni(const int offload, NeighList *list, IntelBuffers *buffers, if (j > vgmax) vgmax = j; } #endif - } - } + } + } #ifdef _LMP_INTEL_OFFLOAD - lmin = MIN(lmin,vlmin); - gmin = MIN(gmin,vgmin); - lmax = MAX(lmax,vlmax); - gmax = MAX(gmax,vgmax); + lmin = MIN(lmin,vlmin); + gmin = MIN(gmin,vgmin); + lmax = MAX(lmax,vlmax); + gmax = MAX(gmax,vgmax); #endif int n = lane, n2 = pack_offset; - for (int u = pack_offset; u < raw_count; u++) { - const int j = neighptr[u]; - int pj = j; - if (pj < e_nall) { - if (need_ic) - if (pj < 0) pj = -pj - 1; + for (int u = pack_offset; u < raw_count; u++) { + const int j = neighptr[u]; + int pj = j; + if (pj < e_nall) { + if (need_ic) + if (pj < 0) pj = -pj - 1; - if (pj < nlocal) { - neighptr[n] = j; - n += swidth; - } else - neighptr[n2++] = j; - } - } - int ns = (n - lane) / swidth; - for (int u = pack_offset; u < n2; u++) { - neighptr[n] = neighptr[u]; - n += swidth; - } + if (pj < nlocal) { + neighptr[n] = j; + n += swidth; + } else + neighptr[n2++] = j; + } + } + int ns = (n - lane) / swidth; + for (int u = pack_offset; u < n2; u++) { + neighptr[n] = neighptr[u]; + n += swidth; + } ilist[i] = i; cnumneigh[i] = ct + lane; - ns += n2 - pack_offset; - #ifndef OUTER_CHUNK + ns += n2 - pack_offset; + #ifndef OUTER_CHUNK int edge = (ns % pad_width); if (edge) { const int pad_end = ns + (pad_width - edge); @@ -458,41 +458,41 @@ hbni(const int offload, NeighList *list, IntelBuffers *buffers, for ( ; ns < pad_end; ns++) neighptr[ns] = e_nall; } - #endif + #endif numneigh[i] = ns; - #ifdef OUTER_CHUNK - if (ns > max_chunk) max_chunk = ns; - lane++; - if (lane == swidth) { - ct += max_chunk * swidth; - const int alignb = (INTEL_DATA_ALIGN / sizeof(int)); - int edge = (ct % alignb); - if (edge) ct += alignb - edge; - neighptr = firstneigh + ct; - max_chunk = 0; - pack_offset = maxnbors * swidth; - lane = 0; - if (ct + obound > list_size) { + #ifdef OUTER_CHUNK + if (ns > max_chunk) max_chunk = ns; + lane++; + if (lane == swidth) { + ct += max_chunk * swidth; + const int alignb = (INTEL_DATA_ALIGN / sizeof(int)); + int edge = (ct % alignb); + if (edge) ct += alignb - edge; + neighptr = firstneigh + ct; + max_chunk = 0; + pack_offset = maxnbors * swidth; + lane = 0; + if (ct + obound > list_size) { if (i < ito - 1) { - *overflow = 1; - ct = (ifrom + tid * 2) * maxnbors; + *overflow = 1; + ct = (ifrom + tid * 2) * maxnbors; } - } - } - #else - ct += ns; - const int alignb = (INTEL_DATA_ALIGN / sizeof(int)); - edge = (ct % alignb); - if (edge) ct += alignb - edge; - neighptr = firstneigh + ct; - if (ct + obound > list_size) { - if (i < ito - 1) { - *overflow = 1; - ct = (ifrom + tid * 2) * maxnbors; - } - } - #endif + } + } + #else + ct += ns; + const int alignb = (INTEL_DATA_ALIGN / sizeof(int)); + edge = (ct % alignb); + if (edge) ct += alignb - edge; + neighptr = firstneigh + ct; + if (ct + obound > list_size) { + if (i < ito - 1) { + *overflow = 1; + ct = (ifrom + tid * 2) * maxnbors; + } + } + #endif } if (*overflow == 1) @@ -505,25 +505,25 @@ hbni(const int offload, NeighList *list, IntelBuffers *buffers, #pragma omp critical #endif { - if (lmin < overflow[LMP_LOCAL_MIN]) overflow[LMP_LOCAL_MIN] = lmin; - if (lmax > overflow[LMP_LOCAL_MAX]) overflow[LMP_LOCAL_MAX] = lmax; - if (gmin < overflow[LMP_GHOST_MIN]) overflow[LMP_GHOST_MIN] = gmin; - if (gmax > overflow[LMP_GHOST_MAX]) overflow[LMP_GHOST_MAX] = gmax; + if (lmin < overflow[LMP_LOCAL_MIN]) overflow[LMP_LOCAL_MIN] = lmin; + if (lmax > overflow[LMP_LOCAL_MAX]) overflow[LMP_LOCAL_MAX] = lmax; + if (gmin < overflow[LMP_GHOST_MIN]) overflow[LMP_GHOST_MIN] = gmin; + if (gmax > overflow[LMP_GHOST_MAX]) overflow[LMP_GHOST_MAX] = gmax; } - #pragma omp barrier + #pragma omp barrier } int ghost_offset = 0, nall_offset = e_nall; if (separate_buffers) { - int nghost = overflow[LMP_GHOST_MAX] + 1 - overflow[LMP_GHOST_MIN]; - if (nghost < 0) nghost = 0; - if (offload) { - ghost_offset = overflow[LMP_GHOST_MIN] - overflow[LMP_LOCAL_MAX] - 1; - nall_offset = overflow[LMP_LOCAL_MAX] + 1 + nghost; - } else { - ghost_offset = overflow[LMP_GHOST_MIN] - nlocal; - nall_offset = nlocal + nghost; - } + int nghost = overflow[LMP_GHOST_MAX] + 1 - overflow[LMP_GHOST_MIN]; + if (nghost < 0) nghost = 0; + if (offload) { + ghost_offset = overflow[LMP_GHOST_MIN] - overflow[LMP_LOCAL_MAX] - 1; + nall_offset = overflow[LMP_LOCAL_MAX] + 1 + nghost; + } else { + ghost_offset = overflow[LMP_GHOST_MIN] - nlocal; + nall_offset = nlocal + nghost; + } } #endif @@ -531,49 +531,49 @@ hbni(const int offload, NeighList *list, IntelBuffers *buffers, for (int i = ifrom; i < ito; ++i) { int * _noalias jlist = firstneigh + cnumneigh[i]; const int jnum = numneigh[i]; - #ifndef OUTER_CHUNK + #ifndef OUTER_CHUNK #if defined(LMP_SIMD_COMPILER) - #pragma vector aligned + #pragma vector aligned #pragma simd - #endif + #endif for (int jj = 0; jj < jnum; jj++) { - #else - const int trip = jnum * swidth; + #else + const int trip = jnum * swidth; for (int jj = 0; jj < trip; jj+= swidth) { - #endif + #endif const int j = jlist[jj]; - if (need_ic && j < 0) { - which = 0; - jlist[jj] = -j - 1; + if (need_ic && j < 0) { + which = 0; + jlist[jj] = -j - 1; } else ofind_special(which, special, nspecial, i, tag[j]); - #ifdef _LMP_INTEL_OFFLOAD - if (j >= nlocal) { - if (j == e_nall) - jlist[jj] = nall_offset; - else if (which) - jlist[jj] = (j-ghost_offset) ^ (which << SBBITS); - else jlist[jj]-=ghost_offset; + #ifdef _LMP_INTEL_OFFLOAD + if (j >= nlocal) { + if (j == e_nall) + jlist[jj] = nall_offset; + else if (which) + jlist[jj] = (j-ghost_offset) ^ (which << SBBITS); + else jlist[jj]-=ghost_offset; } else - #endif + #endif if (which) jlist[jj] = j ^ (which << SBBITS); } } } #ifdef _LMP_INTEL_OFFLOAD else if (separate_buffers) { - for (int i = ifrom; i < ito; ++i) { + for (int i = ifrom; i < ito; ++i) { int * _noalias jlist = firstneigh + cnumneigh[i]; const int jnum = numneigh[i]; - int jj = 0; - for (jj = 0; jj < jnum; jj++) - if (jlist[jj] >= nlocal) break; - while (jj < jnum) { - if (jlist[jj] == e_nall) jlist[jj] = nall_offset; - else jlist[jj] -= ghost_offset; - jj++; - } - } + int jj = 0; + for (jj = 0; jj < jnum; jj++) + if (jlist[jj] >= nlocal) break; + while (jj < jnum) { + if (jlist[jj] == e_nall) jlist[jj] = nall_offset; + else jlist[jj] -= ghost_offset; + jj++; + } + } } #endif } // end omp @@ -597,9 +597,9 @@ hbni(const int offload, NeighList *list, IntelBuffers *buffers, _fix->start_watch(TIME_PACK); _fix->set_neighbor_host_sizes(); buffers->pack_sep_from_single(_fix->host_min_local(), - _fix->host_used_local(), - _fix->host_min_ghost(), - _fix->host_used_ghost()); + _fix->host_used_local(), + _fix->host_min_ghost(), + _fix->host_used_ghost()); _fix->stop_watch(TIME_PACK); } } diff --git a/src/USER-INTEL/npair_half_bin_newton_intel.h b/src/USER-INTEL/npair_half_bin_newton_intel.h index 4e496986b4..9b5d0780a1 100644 --- a/src/USER-INTEL/npair_half_bin_newton_intel.h +++ b/src/USER-INTEL/npair_half_bin_newton_intel.h @@ -37,8 +37,8 @@ class NPairHalfBinNewtonIntel : public NPairIntel { template void hbni(NeighList *, IntelBuffers *); template - void hbni(const int, NeighList *, IntelBuffers *, const int, - const int, const int offload_end = 0); + void hbni(const int, NeighList *, IntelBuffers *, const int, + const int, const int offload_end = 0); }; } diff --git a/src/USER-INTEL/npair_half_bin_newton_tri_intel.cpp b/src/USER-INTEL/npair_half_bin_newton_tri_intel.cpp index 3b6d68d4de..5f191e0797 100644 --- a/src/USER-INTEL/npair_half_bin_newton_tri_intel.cpp +++ b/src/USER-INTEL/npair_half_bin_newton_tri_intel.cpp @@ -26,7 +26,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfBinNewtonTriIntel::NPairHalfBinNewtonTriIntel(LAMMPS *lmp) : +NPairHalfBinNewtonTriIntel::NPairHalfBinNewtonTriIntel(LAMMPS *lmp) : NPairIntel(lmp) {} /* ---------------------------------------------------------------------- @@ -75,7 +75,7 @@ hbnti(NeighList *list, IntelBuffers *buffers) { int need_ic = 0; if (atom->molecular) dminimum_image_check(need_ic, neighbor->cutneighmax, neighbor->cutneighmax, - neighbor->cutneighmax); + neighbor->cutneighmax); #ifdef _LMP_INTEL_OFFLOAD if (need_ic) { @@ -171,7 +171,7 @@ hbnti(const int offload, NeighList *list, IntelBuffers *buffers, overflow = _fix->get_off_overflow_flag(); _fix->stop_watch(TIME_HOST_NEIGHBOR); _fix->start_watch(TIME_OFFLOAD_LATENCY); - } else + } else #endif { tnum = comm->nthreads; @@ -279,20 +279,20 @@ hbnti(const int offload, NeighList *list, IntelBuffers *buffers, const int ibin = atombin[i]; - int raw_count = maxnbors; + int raw_count = maxnbors; for (int k = 0; k < nstencilp; k++) { const int bstart = binhead[ibin + binstart[k]]; const int bend = binhead[ibin + binend[k]]; for (int jj = bstart; jj < bend; jj++) { const int j = binpacked[jj]; - #ifdef _LMP_INTEL_OFFLOAD - if (offload_noghost) { + #ifdef _LMP_INTEL_OFFLOAD + if (offload_noghost) { if (j < nlocal) { if (i < offload_end) continue; } else if (offload) continue; } - #endif + #endif if (x[j].z < ztmp) continue; if (x[j].z == ztmp) { @@ -305,45 +305,45 @@ hbnti(const int offload, NeighList *list, IntelBuffers *buffers, #ifndef _LMP_INTEL_OFFLOAD if (exclude) { - const int jtype = x[j].w; - if (exclusion(i,j,itype,jtype,mask,molecule)) continue; - } - #endif + const int jtype = x[j].w; + if (exclusion(i,j,itype,jtype,mask,molecule)) continue; + } + #endif - neighptr[raw_count++] = j; - } - } - if (raw_count > obound) - *overflow = 1; + neighptr[raw_count++] = j; + } + } + if (raw_count > obound) + *overflow = 1; #if defined(LMP_SIMD_COMPILER) - #ifdef _LMP_INTEL_OFFLOAD - int vlmin = lmin, vlmax = lmax, vgmin = gmin, vgmax = gmax; - #if __INTEL_COMPILER+0 > 1499 - #pragma vector aligned + #ifdef _LMP_INTEL_OFFLOAD + int vlmin = lmin, vlmax = lmax, vgmin = gmin, vgmax = gmax; + #if __INTEL_COMPILER+0 > 1499 + #pragma vector aligned #pragma simd reduction(max:vlmax,vgmax) reduction(min:vlmin, vgmin) - #endif - #else - #pragma vector aligned + #endif + #else + #pragma vector aligned #pragma simd - #endif - #endif - for (int u = maxnbors; u < raw_count; u++) { + #endif + #endif + for (int u = maxnbors; u < raw_count; u++) { int j = neighptr[u]; const flt_t delx = xtmp - x[j].x; const flt_t dely = ytmp - x[j].y; const flt_t delz = ztmp - x[j].z; - const int jtype = x[j].w; + const int jtype = x[j].w; const flt_t rsq = delx * delx + dely * dely + delz * delz; - if (rsq > cutneighsq[ioffset + jtype]) - neighptr[u] = e_nall; - else { + if (rsq > cutneighsq[ioffset + jtype]) + neighptr[u] = e_nall; + else { if (need_ic) { - int no_special; - ominimum_image_check(no_special, delx, dely, delz); - if (no_special) - neighptr[u] = -j - 1; - } + int no_special; + ominimum_image_check(no_special, delx, dely, delz); + if (no_special) + neighptr[u] = -j - 1; + } #ifdef _LMP_INTEL_OFFLOAD if (j < nlocal) { @@ -358,26 +358,26 @@ hbnti(const int offload, NeighList *list, IntelBuffers *buffers, } int n = 0, n2 = maxnbors; - for (int u = maxnbors; u < raw_count; u++) { - const int j = neighptr[u]; - int pj = j; - if (pj < e_nall) { - if (need_ic) - if (pj < 0) pj = -pj - 1; + for (int u = maxnbors; u < raw_count; u++) { + const int j = neighptr[u]; + int pj = j; + if (pj < e_nall) { + if (need_ic) + if (pj < 0) pj = -pj - 1; - if (pj < nlocal) - neighptr[n++] = j; - else - neighptr[n2++] = j; - } - } - int ns = n; - for (int u = maxnbors; u < n2; u++) - neighptr[n++] = neighptr[u]; + if (pj < nlocal) + neighptr[n++] = j; + else + neighptr[n2++] = j; + } + } + int ns = n; + for (int u = maxnbors; u < n2; u++) + neighptr[n++] = neighptr[u]; ilist[i] = i; cnumneigh[i] = ct; - ns += n2 - maxnbors; + ns += n2 - maxnbors; int edge = (ns % pad_width); if (edge) { @@ -390,17 +390,17 @@ hbnti(const int offload, NeighList *list, IntelBuffers *buffers, } numneigh[i] = ns; - ct += ns; - const int alignb = (INTEL_DATA_ALIGN / sizeof(int)); - edge = (ct % alignb); - if (edge) ct += alignb - edge; - neighptr = firstneigh + ct; - if (ct + obound > list_size) { - if (i < ito - 1) { - *overflow = 1; - ct = (ifrom + tid * 2) * maxnbors; - } - } + ct += ns; + const int alignb = (INTEL_DATA_ALIGN / sizeof(int)); + edge = (ct % alignb); + if (edge) ct += alignb - edge; + neighptr = firstneigh + ct; + if (ct + obound > list_size) { + if (i < ito - 1) { + *overflow = 1; + ct = (ifrom + tid * 2) * maxnbors; + } + } } if (*overflow == 1) @@ -428,7 +428,7 @@ hbnti(const int offload, NeighList *list, IntelBuffers *buffers, if (offload) { ghost_offset = overflow[LMP_GHOST_MIN] - overflow[LMP_LOCAL_MAX] - 1; nall_offset = overflow[LMP_LOCAL_MAX] + 1 + nghost; - } else { + } else { ghost_offset = overflow[LMP_GHOST_MIN] - nlocal; nall_offset = nlocal + nghost; } @@ -440,43 +440,43 @@ hbnti(const int offload, NeighList *list, IntelBuffers *buffers, int * _noalias jlist = firstneigh + cnumneigh[i]; const int jnum = numneigh[i]; #if defined(LMP_SIMD_COMPILER) - #pragma vector aligned + #pragma vector aligned #pragma simd - #endif + #endif for (int jj = 0; jj < jnum; jj++) { const int j = jlist[jj]; - if (need_ic && j < 0) { - which = 0; - jlist[jj] = -j - 1; - } else + if (need_ic && j < 0) { + which = 0; + jlist[jj] = -j - 1; + } else ofind_special(which, special, nspecial, i, tag[j]); #ifdef _LMP_INTEL_OFFLOAD - if (j >= nlocal) { - if (j == e_nall) - jlist[jj] = nall_offset; - else if (which) - jlist[jj] = (j-ghost_offset) ^ (which << SBBITS); - else jlist[jj]-=ghost_offset; + if (j >= nlocal) { + if (j == e_nall) + jlist[jj] = nall_offset; + else if (which) + jlist[jj] = (j-ghost_offset) ^ (which << SBBITS); + else jlist[jj]-=ghost_offset; } else #endif - if (which) jlist[jj] = j ^ (which << SBBITS); + if (which) jlist[jj] = j ^ (which << SBBITS); } } } #ifdef _LMP_INTEL_OFFLOAD else if (separate_buffers) { - for (int i = ifrom; i < ito; ++i) { + for (int i = ifrom; i < ito; ++i) { int * _noalias jlist = firstneigh + cnumneigh[i]; const int jnum = numneigh[i]; - int jj = 0; - for (jj = 0; jj < jnum; jj++) - if (jlist[jj] >= nlocal) break; - while (jj < jnum) { - if (jlist[jj] == e_nall) jlist[jj] = nall_offset; - else jlist[jj] -= ghost_offset; - jj++; - } - } + int jj = 0; + for (jj = 0; jj < jnum; jj++) + if (jlist[jj] >= nlocal) break; + while (jj < jnum) { + if (jlist[jj] == e_nall) jlist[jj] = nall_offset; + else jlist[jj] -= ghost_offset; + jj++; + } + } } #endif } // end omp @@ -500,9 +500,9 @@ hbnti(const int offload, NeighList *list, IntelBuffers *buffers, _fix->start_watch(TIME_PACK); _fix->set_neighbor_host_sizes(); buffers->pack_sep_from_single(_fix->host_min_local(), - _fix->host_used_local(), - _fix->host_min_ghost(), - _fix->host_used_ghost()); + _fix->host_used_local(), + _fix->host_min_ghost(), + _fix->host_used_ghost()); _fix->stop_watch(TIME_PACK); } } diff --git a/src/USER-INTEL/npair_half_bin_newton_tri_intel.h b/src/USER-INTEL/npair_half_bin_newton_tri_intel.h index d1b9ee9cd1..d144c2fc52 100644 --- a/src/USER-INTEL/npair_half_bin_newton_tri_intel.h +++ b/src/USER-INTEL/npair_half_bin_newton_tri_intel.h @@ -38,7 +38,7 @@ class NPairHalfBinNewtonTriIntel : public NPairIntel { void hbnti(NeighList *, IntelBuffers *); template void hbnti(const int, NeighList *, IntelBuffers *, const int, - const int, const int offload_end = 0); + const int, const int offload_end = 0); }; } diff --git a/src/USER-INTEL/npair_intel.cpp b/src/USER-INTEL/npair_intel.cpp index bffb31b710..c92ed88774 100644 --- a/src/USER-INTEL/npair_intel.cpp +++ b/src/USER-INTEL/npair_intel.cpp @@ -62,6 +62,6 @@ void NPairIntel::grow_stencil() const int maxstencil = ns->get_maxstencil(); #pragma offload_transfer target(mic:_cop) \ in(stencil:length(maxstencil) alloc_if(1) free_if(0)) - } + } } #endif diff --git a/src/npair.cpp b/src/npair.cpp index 3451cd6eae..9fbb4d219d 100644 --- a/src/npair.cpp +++ b/src/npair.cpp @@ -123,7 +123,7 @@ void NPair::copy_bin_info() mbinxlo = nb->mbinxlo; mbinylo = nb->mbinylo; mbinzlo = nb->mbinzlo; - + bininvx = nb->bininvx; bininvy = nb->bininvy; bininvz = nb->bininvz; @@ -183,15 +183,15 @@ int NPair::exclusion(int i, int j, int itype, int jtype, if (nex_mol) { for (m = 0; m < nex_mol; m++) - // intra-chain: exclude i-j pair if in same molecule - // inter-chain: exclude i-j pair if in different molecules + // intra-chain: exclude i-j pair if in same molecule + // inter-chain: exclude i-j pair if in different molecules if (ex_mol_intra[m]) { if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] && - molecule[i] == molecule[j]) return 1; + molecule[i] == molecule[j]) return 1; } else { - if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] && - molecule[i] != molecule[j]) return 1; + if (mask[i] & ex_mol_bit[m] && mask[j] & ex_mol_bit[m] && + molecule[i] != molecule[j]) return 1; } } diff --git a/src/npair.h b/src/npair.h index 4e5e3f5dfd..6941b86164 100644 --- a/src/npair.h +++ b/src/npair.h @@ -79,7 +79,7 @@ class NPair : protected Pointers { double bininvx,bininvy,bininvz; int *atom2bin,*bins; int *binhead; - + // data from NStencil class int nstencil; diff --git a/src/npair_full_bin.h b/src/npair_full_bin.h index 432fb3cbf8..56c338e360 100644 --- a/src/npair_full_bin.h +++ b/src/npair_full_bin.h @@ -15,7 +15,7 @@ NPairStyle(full/bin, NPairFullBin, - NP_FULL | NP_BIN | NP_MOLONLY | + NP_FULL | NP_BIN | NP_MOLONLY | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) #else diff --git a/src/npair_full_bin_atomonly.h b/src/npair_full_bin_atomonly.h index f8c33d9558..0845d1ecef 100644 --- a/src/npair_full_bin_atomonly.h +++ b/src/npair_full_bin_atomonly.h @@ -16,7 +16,7 @@ NPairStyle(full/bin/atomonly, NPairFullBinAtomonly, NP_FULL | NP_BIN | NP_ATOMONLY | - NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) + NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) #else diff --git a/src/npair_full_bin_ghost.h b/src/npair_full_bin_ghost.h index a09aab8512..c5a86e68af 100644 --- a/src/npair_full_bin_ghost.h +++ b/src/npair_full_bin_ghost.h @@ -15,7 +15,7 @@ NPairStyle(full/bin/ghost, NPairFullBinGhost, - NP_FULL | NP_BIN | NP_GHOST | NP_NEWTON | NP_NEWTOFF | + NP_FULL | NP_BIN | NP_GHOST | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) #else diff --git a/src/npair_full_nsq_ghost.h b/src/npair_full_nsq_ghost.h index 3e259ed098..58cd73c392 100644 --- a/src/npair_full_nsq_ghost.h +++ b/src/npair_full_nsq_ghost.h @@ -15,7 +15,7 @@ NPairStyle(full/nsq/ghost, NPairFullNsqGhost, - NP_FULL | NP_NSQ | NP_GHOST | NP_NEWTON | NP_NEWTOFF | + NP_FULL | NP_NSQ | NP_GHOST | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) #else diff --git a/src/npair_half_bin_atomonly_newton.cpp b/src/npair_half_bin_atomonly_newton.cpp index 6bbef0700a..6da44b4a5c 100644 --- a/src/npair_half_bin_atomonly_newton.cpp +++ b/src/npair_half_bin_atomonly_newton.cpp @@ -25,7 +25,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfBinAtomonlyNewton::NPairHalfBinAtomonlyNewton(LAMMPS *lmp) : +NPairHalfBinAtomonlyNewton::NPairHalfBinAtomonlyNewton(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- diff --git a/src/npair_half_respa_bin_newton_tri.cpp b/src/npair_half_respa_bin_newton_tri.cpp index 38621224c4..4ec6685e1d 100644 --- a/src/npair_half_respa_bin_newton_tri.cpp +++ b/src/npair_half_respa_bin_newton_tri.cpp @@ -25,7 +25,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfRespaBinNewtonTri::NPairHalfRespaBinNewtonTri(LAMMPS *lmp) : +NPairHalfRespaBinNewtonTri::NPairHalfRespaBinNewtonTri(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- diff --git a/src/npair_half_size_bin_newton.cpp b/src/npair_half_size_bin_newton.cpp index a8be3ce691..4f4ecccb16 100644 --- a/src/npair_half_size_bin_newton.cpp +++ b/src/npair_half_size_bin_newton.cpp @@ -190,7 +190,7 @@ void NPairHalfSizeBinNewton::build(NeighList *list) nn += dnum; } } - + n++; } } diff --git a/src/npair_half_size_bin_newton_tri.cpp b/src/npair_half_size_bin_newton_tri.cpp index 1107f73026..559eb09a7a 100644 --- a/src/npair_half_size_bin_newton_tri.cpp +++ b/src/npair_half_size_bin_newton_tri.cpp @@ -27,7 +27,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairHalfSizeBinNewtonTri::NPairHalfSizeBinNewtonTri(LAMMPS *lmp) : +NPairHalfSizeBinNewtonTri::NPairHalfSizeBinNewtonTri(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- diff --git a/src/npair_skip_respa.h b/src/npair_skip_respa.h index 62077f85df..deff301909 100644 --- a/src/npair_skip_respa.h +++ b/src/npair_skip_respa.h @@ -15,8 +15,8 @@ NPairStyle(skip/half/respa, NPairSkipRespa, - NP_SKIP | NP_RESPA | NP_HALF | NP_FULL | - NP_NSQ | NP_BIN | NP_MULTI | + NP_SKIP | NP_RESPA | NP_HALF | NP_FULL | + NP_NSQ | NP_BIN | NP_MULTI | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) #else diff --git a/src/npair_skip_size.cpp b/src/npair_skip_size.cpp index 98e757e5c7..e8d19dedca 100644 --- a/src/npair_skip_size.cpp +++ b/src/npair_skip_size.cpp @@ -32,7 +32,7 @@ NPairSkipSize::NPairSkipSize(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- build skip list for subset of types from parent list iskip and ijskip flag which atom types and type pairs to skip - if list requests it, preserve shear history via fix shear/history + if list requests it, preserve shear history via fix shear/history ------------------------------------------------------------------------- */ void NPairSkipSize::build(NeighList *list) diff --git a/src/npair_skip_size.h b/src/npair_skip_size.h index 9573396641..b462c9dc97 100644 --- a/src/npair_skip_size.h +++ b/src/npair_skip_size.h @@ -15,7 +15,7 @@ NPairStyle(skip/half/size, NPairSkipSize, - NP_SKIP | NP_SIZE | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | + NP_SKIP | NP_SIZE | NP_HALF | NP_FULL | NP_NSQ | NP_BIN | NP_MULTI | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) #else diff --git a/src/npair_skip_size_off2on.cpp b/src/npair_skip_size_off2on.cpp index 996e9939df..da9dd57047 100644 --- a/src/npair_skip_size_off2on.cpp +++ b/src/npair_skip_size_off2on.cpp @@ -33,7 +33,7 @@ NPairSkipSizeOff2on::NPairSkipSizeOff2on(LAMMPS *lmp) : NPair(lmp) {} build skip list for subset of types from parent list iskip and ijskip flag which atom types and type pairs to skip parent non-skip list used newton off, this skip list is newton on - if list requests it, preserve shear history via fix shear/history + if list requests it, preserve shear history via fix shear/history ------------------------------------------------------------------------- */ void NPairSkipSizeOff2on::build(NeighList *list) diff --git a/src/npair_skip_size_off2on.h b/src/npair_skip_size_off2on.h index 4b4e9a9c29..dab32f04ff 100644 --- a/src/npair_skip_size_off2on.h +++ b/src/npair_skip_size_off2on.h @@ -15,8 +15,8 @@ NPairStyle(skip/size/off2on, NPairSkipSizeOff2on, - NP_SKIP | NP_SIZE | NP_OFF2ON | NP_HALF | - NP_NSQ | NP_BIN | NP_MULTI | + NP_SKIP | NP_SIZE | NP_OFF2ON | NP_HALF | + NP_NSQ | NP_BIN | NP_MULTI | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) #else diff --git a/src/npair_skip_size_off2on_oneside.cpp b/src/npair_skip_size_off2on_oneside.cpp index a4c1625590..7377feec5b 100644 --- a/src/npair_skip_size_off2on_oneside.cpp +++ b/src/npair_skip_size_off2on_oneside.cpp @@ -27,7 +27,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NPairSkipSizeOff2onOneside::NPairSkipSizeOff2onOneside(LAMMPS *lmp) : +NPairSkipSizeOff2onOneside::NPairSkipSizeOff2onOneside(LAMMPS *lmp) : NPair(lmp) {} /* ---------------------------------------------------------------------- @@ -35,7 +35,7 @@ NPairSkipSizeOff2onOneside::NPairSkipSizeOff2onOneside(LAMMPS *lmp) : iskip and ijskip flag which atom types and type pairs to skip parent non-skip list used newton off and was not onesided, this skip list is newton on and onesided - if list requests it, preserve shear history via fix shear/history + if list requests it, preserve shear history via fix shear/history ------------------------------------------------------------------------- */ void NPairSkipSizeOff2onOneside::build(NeighList *list) diff --git a/src/npair_skip_size_off2on_oneside.h b/src/npair_skip_size_off2on_oneside.h index 9f3c06e7bc..73448ca279 100644 --- a/src/npair_skip_size_off2on_oneside.h +++ b/src/npair_skip_size_off2on_oneside.h @@ -15,8 +15,8 @@ NPairStyle(skip/size/off2on/oneside, NPairSkipSizeOff2onOneside, - NP_SKIP | NP_SIZE | NP_OFF2ON | NP_ONESIDE | NP_HALF | - NP_NSQ | NP_BIN | NP_MULTI | NP_NEWTON | NP_NEWTOFF | + NP_SKIP | NP_SIZE | NP_OFF2ON | NP_ONESIDE | NP_HALF | + NP_NSQ | NP_BIN | NP_MULTI | NP_NEWTON | NP_NEWTOFF | NP_ORTHO | NP_TRI) #else diff --git a/src/nstencil_full_bin_2d.cpp b/src/nstencil_full_bin_2d.cpp index 1f2b666dfb..0986c40dd4 100644 --- a/src/nstencil_full_bin_2d.cpp +++ b/src/nstencil_full_bin_2d.cpp @@ -28,7 +28,7 @@ NStencilFullBin2d::NStencilFullBin2d(LAMMPS *lmp) : NStencil(lmp) {} void NStencilFullBin2d::create() { int i,j; - + nstencil = 0; for (j = -sy; j <= sy; j++) diff --git a/src/nstencil_full_bin_2d.h b/src/nstencil_full_bin_2d.h index 18f848f275..d85063596f 100644 --- a/src/nstencil_full_bin_2d.h +++ b/src/nstencil_full_bin_2d.h @@ -15,7 +15,7 @@ NStencilStyle(full/bin/2d, NStencilFullBin2d, - NS_FULL | NS_BIN | NS_2D | + NS_FULL | NS_BIN | NS_2D | NS_NEWTON | NS_NEWTOFF | NS_ORTHO | NS_TRI) #else diff --git a/src/nstencil_full_bin_3d.h b/src/nstencil_full_bin_3d.h index d9acc9c535..facddd8ead 100644 --- a/src/nstencil_full_bin_3d.h +++ b/src/nstencil_full_bin_3d.h @@ -15,7 +15,7 @@ NStencilStyle(full/bin/3d, NStencilFullBin3d, - NS_FULL | NS_BIN | NS_3D | + NS_FULL | NS_BIN | NS_3D | NS_NEWTON | NS_NEWTOFF | NS_ORTHO | NS_TRI) #else diff --git a/src/nstencil_full_ghost_bin_2d.h b/src/nstencil_full_ghost_bin_2d.h index af47913e7f..531c7d2eb1 100644 --- a/src/nstencil_full_ghost_bin_2d.h +++ b/src/nstencil_full_ghost_bin_2d.h @@ -15,7 +15,7 @@ NStencilStyle(full/ghost/bin/2d, NStencilFullGhostBin2d, - NS_FULL | NS_GHOST | NS_BIN | NS_2D | + NS_FULL | NS_GHOST | NS_BIN | NS_2D | NS_NEWTON | NS_NEWTOFF | NS_ORTHO | NS_TRI) #else diff --git a/src/nstencil_full_ghost_bin_3d.h b/src/nstencil_full_ghost_bin_3d.h index beca6573de..ed4ca6c4d6 100644 --- a/src/nstencil_full_ghost_bin_3d.h +++ b/src/nstencil_full_ghost_bin_3d.h @@ -15,7 +15,7 @@ NStencilStyle(full/ghost/bin/3d, NStencilFullGhostBin3d, - NS_FULL | NS_GHOST | NS_BIN | NS_3D | + NS_FULL | NS_GHOST | NS_BIN | NS_3D | NS_NEWTON | NS_NEWTOFF | NS_ORTHO | NS_TRI) #else diff --git a/src/nstencil_full_multi_2d.h b/src/nstencil_full_multi_2d.h index 8154144eda..f78eecc55f 100644 --- a/src/nstencil_full_multi_2d.h +++ b/src/nstencil_full_multi_2d.h @@ -15,7 +15,7 @@ NStencilStyle(full/multi/2d, NStencilFullMulti2d, - NS_FULL | NS_MULTI | NS_2D | + NS_FULL | NS_MULTI | NS_2D | NS_NEWTON | NS_NEWTOFF | NS_ORTHO | NS_TRI) #else diff --git a/src/nstencil_half_bin_2d_newtoff.cpp b/src/nstencil_half_bin_2d_newtoff.cpp index be5bc81dbf..e51db6fe7a 100644 --- a/src/nstencil_half_bin_2d_newtoff.cpp +++ b/src/nstencil_half_bin_2d_newtoff.cpp @@ -19,7 +19,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NStencilHalfBin2dNewtoff::NStencilHalfBin2dNewtoff(LAMMPS *lmp) : +NStencilHalfBin2dNewtoff::NStencilHalfBin2dNewtoff(LAMMPS *lmp) : NStencil(lmp) {} /* ---------------------------------------------------------------------- diff --git a/src/nstencil_half_bin_2d_newton_tri.cpp b/src/nstencil_half_bin_2d_newton_tri.cpp index 3a645a7434..4f89b1c326 100644 --- a/src/nstencil_half_bin_2d_newton_tri.cpp +++ b/src/nstencil_half_bin_2d_newton_tri.cpp @@ -19,7 +19,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NStencilHalfBin2dNewtonTri::NStencilHalfBin2dNewtonTri(LAMMPS *lmp) : +NStencilHalfBin2dNewtonTri::NStencilHalfBin2dNewtonTri(LAMMPS *lmp) : NStencil(lmp) {} /* ---------------------------------------------------------------------- diff --git a/src/nstencil_half_bin_3d_newtoff.cpp b/src/nstencil_half_bin_3d_newtoff.cpp index 44678b05df..433de400c2 100644 --- a/src/nstencil_half_bin_3d_newtoff.cpp +++ b/src/nstencil_half_bin_3d_newtoff.cpp @@ -19,7 +19,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NStencilHalfBin3dNewtoff::NStencilHalfBin3dNewtoff(LAMMPS *lmp) : +NStencilHalfBin3dNewtoff::NStencilHalfBin3dNewtoff(LAMMPS *lmp) : NStencil(lmp) {} /* ---------------------------------------------------------------------- diff --git a/src/nstencil_half_bin_3d_newton_tri.cpp b/src/nstencil_half_bin_3d_newton_tri.cpp index 9e8c41f97a..691ce0bb80 100644 --- a/src/nstencil_half_bin_3d_newton_tri.cpp +++ b/src/nstencil_half_bin_3d_newton_tri.cpp @@ -19,7 +19,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -NStencilHalfBin3dNewtonTri::NStencilHalfBin3dNewtonTri(LAMMPS *lmp) : +NStencilHalfBin3dNewtonTri::NStencilHalfBin3dNewtonTri(LAMMPS *lmp) : NStencil(lmp) {} /* ---------------------------------------------------------------------- diff --git a/src/nstencil_half_ghost_bin_2d_newtoff.h b/src/nstencil_half_ghost_bin_2d_newtoff.h index 3b70f0042a..3286810c1c 100644 --- a/src/nstencil_half_ghost_bin_2d_newtoff.h +++ b/src/nstencil_half_ghost_bin_2d_newtoff.h @@ -15,7 +15,7 @@ NStencilStyle(half/ghost/bin/2d/newtoff, NStencilHalfGhostBin2dNewtoff, - NS_HALF | NS_GHOST | NS_BIN | NS_2D | + NS_HALF | NS_GHOST | NS_BIN | NS_2D | NS_NEWTOFF | NS_ORTHO | NS_TRI) #else From 167a51538e1689c45d045afa457a3514c50e866b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 30 May 2017 21:52:32 -0400 Subject: [PATCH 03/22] support atom style variables for assigning image flags with the set command --- doc/src/set.txt | 1 + src/set.cpp | 12 +++++++++--- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/doc/src/set.txt b/doc/src/set.txt index 6b59bf1332..0b428d56ed 100644 --- a/doc/src/set.txt +++ b/doc/src/set.txt @@ -80,6 +80,7 @@ keyword = {type} or {type/fraction} or {mol} or {x} or {y} or {z} or \ value can be an atom-style variable (see below) {image} nx ny nz nx,ny,nz = which periodic image of the simulation box the atom is in + any of nx,ny,nz can be an atom-style variable (see below) {bond} value = bond type for all bonds between selected atoms {angle} value = angle type for all angles between selected atoms {dihedral} value = dihedral type for all dihedrals between selected atoms diff --git a/src/set.cpp b/src/set.cpp index 4ed07d423b..b97a96167c 100644 --- a/src/set.cpp +++ b/src/set.cpp @@ -327,15 +327,18 @@ void Set::command(int narg, char **arg) ximageflag = yimageflag = zimageflag = 0; if (strcmp(arg[iarg+1],"NULL") != 0) { ximageflag = 1; - ximage = force->inumeric(FLERR,arg[iarg+1]); + if (strstr(arg[iarg+1],"v_") == arg[iarg+1]) varparse(arg[iarg+1],1); + else ximage = force->inumeric(FLERR,arg[iarg+1]); } if (strcmp(arg[iarg+2],"NULL") != 0) { yimageflag = 1; - yimage = force->inumeric(FLERR,arg[iarg+2]); + if (strstr(arg[iarg+2],"v_") == arg[iarg+2]) varparse(arg[iarg+2],2); + else yimage = force->inumeric(FLERR,arg[iarg+2]); } if (strcmp(arg[iarg+3],"NULL") != 0) { zimageflag = 1; - zimage = force->inumeric(FLERR,arg[iarg+3]); + if (strstr(arg[iarg+3],"v_") == arg[iarg+3]) varparse(arg[iarg+3],3); + else zimage = force->inumeric(FLERR,arg[iarg+3]); } if (ximageflag && ximage && !domain->xperiodic) error->all(FLERR, @@ -789,6 +792,9 @@ void Set::set(int keyword) int xbox = (atom->image[i] & IMGMASK) - IMGMAX; int ybox = (atom->image[i] >> IMGBITS & IMGMASK) - IMGMAX; int zbox = (atom->image[i] >> IMG2BITS) - IMGMAX; + if (varflag1) ximage = static_cast(xvalue); + if (varflag2) yimage = static_cast(yvalue); + if (varflag3) zimage = static_cast(zvalue); if (ximageflag) xbox = ximage; if (yimageflag) ybox = yimage; if (zimageflag) zbox = zimage; From 066123007cd9277d4c03669707cf34ab6c05bc06 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 30 May 2017 21:54:16 -0400 Subject: [PATCH 04/22] avoid undesired negative forces for high particle velocities in granular models --- src/GRANULAR/pair_gran_hertz_history.cpp | 2 ++ src/GRANULAR/pair_gran_hooke.cpp | 2 ++ src/GRANULAR/pair_gran_hooke_history.cpp | 2 ++ 3 files changed, 6 insertions(+) diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index e52aac10db..14cbc8eb0a 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -183,6 +183,7 @@ void PairGranHertzHistory::compute(int eflag, int vflag) ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; + if (ccel < 0.0) ccel = 0.0; // relative velocities @@ -390,6 +391,7 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype, ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; + if (ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 66e1c3dd7a..2d6d113ab7 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -161,6 +161,7 @@ void PairGranHooke::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (ccel < 0.0) ccel = 0.0; // relative velocities @@ -302,6 +303,7 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index e9662c9e73..1378709cb5 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -223,6 +223,7 @@ void PairGranHookeHistory::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (ccel < 0.0) ccel = 0.0; // relative velocities @@ -687,6 +688,7 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; + if (ccel < 0.0) ccel = 0.0; // relative velocities From 0be2cd3d432dcfd122651ab7ac68f3f4e24d1fac Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 30 May 2017 23:58:56 -0400 Subject: [PATCH 05/22] fix bug reported on lammps-users, when not using the first molecule template --- src/MC/fix_gcmc.cpp | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/src/MC/fix_gcmc.cpp b/src/MC/fix_gcmc.cpp index e0a269b88c..fbe6b6bb62 100644 --- a/src/MC/fix_gcmc.cpp +++ b/src/MC/fix_gcmc.cpp @@ -501,7 +501,7 @@ void FixGCMC::init() if (ifix < 0) error->all(FLERR,"Fix gcmc rigid fix does not exist"); fixrigid = modify->fix[ifix]; int tmp; - if (onemols != (Molecule **) fixrigid->extract("onemol",tmp)) + if (&onemols[imol] != (Molecule **) fixrigid->extract("onemol",tmp)) error->all(FLERR, "Fix gcmc and fix rigid/small not using " "same molecule template ID"); @@ -516,7 +516,7 @@ void FixGCMC::init() if (ifix < 0) error->all(FLERR,"Fix gcmc shake fix does not exist"); fixshake = modify->fix[ifix]; int tmp; - if (onemols != (Molecule **) fixshake->extract("onemol",tmp)) + if (&onemols[imol] != (Molecule **) fixshake->extract("onemol",tmp)) error->all(FLERR,"Fix gcmc and fix shake not using " "same molecule template ID"); } @@ -1397,12 +1397,13 @@ void FixGCMC::attempt_molecule_insertion() // FixRigidSmall::set_molecule stores rigid body attributes // FixShake::set_molecule stores shake info for molecule - - if (rigidflag) - fixrigid->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat); - else if (shakeflag) - fixshake->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat); + for (int submol = 0; submol < nmol; ++submol) { + if (rigidflag) + fixrigid->set_molecule(nlocalprev,maxtag_all,submol,com_coord,vnew,quat); + else if (shakeflag) + fixshake->set_molecule(nlocalprev,maxtag_all,submol,com_coord,vnew,quat); + } atom->natoms += natoms_per_molecule; if (atom->natoms < 0) error->all(FLERR,"Too many total atoms"); @@ -2058,11 +2059,12 @@ void FixGCMC::attempt_molecule_insertion_full() // FixRigidSmall::set_molecule stores rigid body attributes // FixShake::set_molecule stores shake info for molecule - if (rigidflag) - fixrigid->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat); - else if (shakeflag) - fixshake->set_molecule(nlocalprev,maxtag_all,imol,com_coord,vnew,quat); - + for (int submol = 0; submol < nmol; ++submol) { + if (rigidflag) + fixrigid->set_molecule(nlocalprev,maxtag_all,submol,com_coord,vnew,quat); + else if (shakeflag) + fixshake->set_molecule(nlocalprev,maxtag_all,submol,com_coord,vnew,quat); + } atom->natoms += natoms_per_molecule; if (atom->natoms < 0) error->all(FLERR,"Too many total atoms"); From f57f1efdff0cea43fcad22b4061a3c944fa91110 Mon Sep 17 00:00:00 2001 From: Anders Hafreager Date: Wed, 31 May 2017 00:34:26 -0700 Subject: [PATCH 06/22] Setting lattice to NULL before creating --- src/domain.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/domain.cpp b/src/domain.cpp index 8ead12cd4e..aadd1b751f 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -1697,6 +1697,7 @@ int Domain::ownatom(int id, double *x, imageint *image, int shrinkexceed) void Domain::set_lattice(int narg, char **arg) { if (lattice) delete lattice; + lattice = NULL; lattice = new Lattice(lmp,narg,arg); } From af5f19604cceb62f547132c723de6ea49d9aaf6e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 31 May 2017 23:36:39 -0400 Subject: [PATCH 07/22] remove no longer correct sentence from set command docs --- doc/src/set.txt | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/doc/src/set.txt b/doc/src/set.txt index 0b428d56ed..14460c9741 100644 --- a/doc/src/set.txt +++ b/doc/src/set.txt @@ -364,9 +364,8 @@ A value of -1 means subtract 1 box length to get the true value. LAMMPS updates these flags as atoms cross periodic boundaries during the simulation. The flags can be output with atom snapshots via the "dump"_dump.html command. If a value of NULL is specified for any of -nx,ny,nz, then the current image value for that dimension is -unchanged. For non-periodic dimensions only a value of 0 can be -specified. This keyword does not allow use of atom-style variables. +nx,ny,nz, then the current image value for that dimension is unchanged. +For non-periodic dimensions only a value of 0 can be specified. This command can be useful after a system has been equilibrated and atoms have diffused one or more box lengths in various directions. This command can then reset the image values for atoms so that they From f59ee5bd62c7a830991ae63c54c9b45cab5ef4a2 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 2 Jun 2017 08:45:15 -0400 Subject: [PATCH 08/22] enable support for dynamic groups in fix planeforce and fix lineforce --- src/fix_lineforce.cpp | 2 ++ src/fix_planeforce.cpp | 2 ++ 2 files changed, 4 insertions(+) diff --git a/src/fix_lineforce.cpp b/src/fix_lineforce.cpp index f82ed957f7..1e78bf3ec3 100644 --- a/src/fix_lineforce.cpp +++ b/src/fix_lineforce.cpp @@ -29,6 +29,8 @@ using namespace FixConst; FixLineForce::FixLineForce(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { + dynamic_group_allow = 1; + if (narg != 6) error->all(FLERR,"Illegal fix lineforce command"); xdir = force->numeric(FLERR,arg[3]); ydir = force->numeric(FLERR,arg[4]); diff --git a/src/fix_planeforce.cpp b/src/fix_planeforce.cpp index 872bd98716..5e999c888c 100644 --- a/src/fix_planeforce.cpp +++ b/src/fix_planeforce.cpp @@ -29,6 +29,8 @@ using namespace FixConst; FixPlaneForce::FixPlaneForce(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { + dynamic_group_allow = 1; + if (narg != 6) error->all(FLERR,"Illegal fix planeforce command"); xdir = force->numeric(FLERR,arg[3]); ydir = force->numeric(FLERR,arg[4]); From ff58ccac2870341b1ca911f056853c150d510394 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sun, 4 Jun 2017 21:21:32 -0400 Subject: [PATCH 09/22] add clarification to impact of special bonds to manybody potentials --- doc/src/special_bonds.txt | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/doc/src/special_bonds.txt b/doc/src/special_bonds.txt index 6924b321a0..6a661015bd 100644 --- a/doc/src/special_bonds.txt +++ b/doc/src/special_bonds.txt @@ -65,7 +65,13 @@ sense to define permanent bonds between atoms that interact via these potentials, though such bonds may exist elsewhere in your system, e.g. when using the "pair_style hybrid"_pair_hybrid.html command. Thus LAMMPS ignores special_bonds settings when manybody potentials -are calculated. +are calculated. Please note, that the existence of explicit bonds +for atoms that are described by a manybody potential will alter the +neigborlist and thus can render the computation of those interactions +invalid, since those pairs are not only used to determine direct +pairwise interactions but also neighbors of neighbors and more. +The recommended course of action is to remove such bonds, or - if +that is not possible - use a special bonds setting of 1.0 1.0 1.0. NOTE: Unlike some commands in LAMMPS, you cannot use this command multiple times in an incremental fashion: e.g. to first set the LJ From d437650c776f218c139c2900e7fa8efcad6a2750 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 6 Jun 2017 08:08:10 -0400 Subject: [PATCH 10/22] make certain Domain::box_change is initialized before use --- src/domain.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/src/domain.cpp b/src/domain.cpp index aadd1b751f..427f7785e8 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -60,6 +60,7 @@ enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files Domain::Domain(LAMMPS *lmp) : Pointers(lmp) { box_exist = 0; + box_change = 0; dimension = 3; nonperiodic = 0; From cd67eaa5f42a4b6506dd95b871de872689757e17 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 6 Jun 2017 16:26:57 -0400 Subject: [PATCH 11/22] update e-mail and affiliation for stefan paquay in USER-MANIFOLD related files --- doc/src/Section_packages.txt | 4 ++-- src/USER-MANIFOLD/README | 8 ++++++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/doc/src/Section_packages.txt b/doc/src/Section_packages.txt index cc44c05906..d4276fdf82 100644 --- a/doc/src/Section_packages.txt +++ b/doc/src/Section_packages.txt @@ -2027,8 +2027,8 @@ algorithm to formulate single-particle constraint functions g(xi,yi,zi) = 0 and their derivative (i.e. the normal of the manifold) n = grad(g). -[Author:] Stefan Paquay (Eindhoven University of Technology (TU/e), The -Netherlands) +[Author:] Stefan Paquay (until 2017: Eindhoven University of Technology (TU/e), The +Netherlands; since 2017: Brandeis University, Waltham, MA, USA) [Install or un-install:] diff --git a/src/USER-MANIFOLD/README b/src/USER-MANIFOLD/README index f55a9bb8e3..eb83cfc5ab 100644 --- a/src/USER-MANIFOLD/README +++ b/src/USER-MANIFOLD/README @@ -7,10 +7,14 @@ box). It achieves this using the RATTLE constraint algorithm applied to single-particle constraint functions g(xi,yi,zi) = 0 and their derivative (i.e. the normal of the manifold) n = grad(g). -Stefan Paquay, s.paquay@tue.nl -Applied Physics/Theory of Polymers and Soft Matter, +Stefan Paquay, stefanpaquay@gmail.com + +until 2017: Applied Physics/Theory of Polymers and Soft Matter, Eindhoven University of Technology (TU/e), The Netherlands +since 2017: Brandeis University, Waltham, MA, USA. + + Thanks to Remy Kusters at TU/e for testing. This software is distributed under the GNU General Public License. From 5cb56796a22fcc615d613a612a0d9edc2a248026 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 6 Jun 2017 17:26:06 -0400 Subject: [PATCH 12/22] alias pair style lj/sf to lj/smooth/linear and remove/update related files --- doc/src/Eqs/pair_lj_sf.jpg | Bin 15535 -> 0 bytes doc/src/Eqs/pair_lj_sf.tex | 11 - doc/src/Section_commands.txt | 1 - doc/src/pair_lj_sf.txt | 114 -------- doc/src/pair_lj_smooth_linear.txt | 31 +- src/Purge.list | 5 + src/USER-MISC/pair_lj_sf.cpp | 355 ----------------------- src/USER-MISC/pair_lj_sf.h | 53 ---- src/USER-OMP/pair_lj_sf_omp.cpp | 164 ----------- src/USER-OMP/pair_lj_sf_omp.h | 48 --- src/USER-OMP/pair_lj_smooth_linear_omp.h | 1 + src/pair_lj_smooth_linear.h | 1 + 12 files changed, 23 insertions(+), 761 deletions(-) delete mode 100644 doc/src/Eqs/pair_lj_sf.jpg delete mode 100644 doc/src/Eqs/pair_lj_sf.tex delete mode 100644 doc/src/pair_lj_sf.txt delete mode 100644 src/USER-MISC/pair_lj_sf.cpp delete mode 100644 src/USER-MISC/pair_lj_sf.h delete mode 100644 src/USER-OMP/pair_lj_sf_omp.cpp delete mode 100644 src/USER-OMP/pair_lj_sf_omp.h diff --git a/doc/src/Eqs/pair_lj_sf.jpg b/doc/src/Eqs/pair_lj_sf.jpg deleted file mode 100644 index a702240003cd4e1bd5043fc5e6a523998fadf8de..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 15535 zcmd6O1yEamv}Pz?tk9w%NwLzR#U0wxq6JD@+})uNJXoPP1S{@d3c=ke65QP-xCf`m z@_(~CJMZnjeLM4Z-t3o|Br}n7@9)U(e&0Fwe)@hH@I+2ZRtkWDfdPwIfX90Ns7?J?o^}YX3;FETC!m%yKR&sa6IG3r7oA=X~ zW~zfe!0X7bo=ar4`+zkH?{=-@%Ec?-`A?xhLA@x#b34@qvRGSrS?y|mnIG9+Km>9ao=;CzSpLnALPW2 zaF*oqddx@VxSi)BFKy`Q{oLEZ7qO~@S)o$Bi%zkgSX|8Oqiy>5?b=tw;^T4JTR+rv zTg=*(W}@7OgNPx-OPFx-4!_u+m=6#f-Q2MKFfng#jh&7ewkE_aEz^=c+-0}cW`V&S z@U%GZOoo2BA+A84T^?jE`}EL>zcIHZ+anyLtEX1c_x2m^IN$6kYZK{z2!nfoc@Y8# zW_H!+D4B;hF(kx@xD83LYAA9W5?=?3T~CdBJ-$^F2s)-D>3uyNOkV|uGp{jKZ^TP& zFOhQdWPBf|n(C5W47LZ=!{#jCf;3Af{F}Jmb#=9!0?av8PFQE2*+xs_l=;*RGM8eG z|6V%22L#*$ig;UAC9XNAuMSC~rMd>+Dh%3VWIUa}TQ_;H?mf%ODuoPQq_AhbU?`E< zU>sb84)soc2}|&6&K&ZpI)nJOs0^%%erQLB@zcv7jIx6-du-}*N9tWVK89^Qcf#R2 zj&{k?mbOhhpN*7OE#fn<`}nD>tke@gLoo6A#=ONaA1f~W6{-Zv<4_F7H)0g$v2=d7 zSqWdRih#9&Q$^0dBoe3yZgBXs@CqkakLRnBP4ePM=^$!-Cac2i5BC;`qA?KqqRO3S zzmGivI0;gcfjwaA&m#qyz{VQUD?gZ(w5)|o-r1-&zw}gavF%+lY4+=u_QC?5#>+yqE+^VBCZJ^VQ9L+9}eVnK-1toz{$nc zt^0CN7|gRW1Zp5w>}30*_3B6Dv{f6>h&u*BX3SNi6THKG&Fs9$AfjCzi(InSfNG{b^U@*HYVz_UBeoJ-2PjH%#TTx@6<)Zno9x zeX7ph^sD&ZU-wVyop#Z$%U_H6(TK&TmOR^P(~7o=FwJl13ibA9dHa;n(O|e-^|Y>- z9TuvKD=$medMLRe&>KU*XHYsR3p6p(1v0p6FNAMJ7+ta~g;LFwr-3%I@-Cz!OEf08 zfSVBhstQ6m-b)AE=jTYpOKUTWkaRv|m8J;fZ&8%)EK4sK*J^cGMnw}m$WTF)k z@M!k1LM+kam!VbCViwJ2wDv0I8!K8_355r56fFN=dE>*`2M1kgzcnYLvgDqroSF0y zOt1^%;jPRx|2F7W)z=&QQGX3L>OCPCl;Q$bf?dT?YPen_RdO8t&SLkTUw17RTW0on zyy?{E`FULPYB8#Yo2W~Ymyl$qdP|4LXbZPc2VgV1jXi*MyE;S6)RUZE75~egn$2;# z@LgrrxTf`wvXXgQtUw=H42J5q5r+3eSr66h-@#=WFJEDgVMkdz;YZGBSz4M&UZQB( zJpgzYe|bJI0_|O;_Hu}JLukb4^~7seof*^OL|U`N;ft$2ZA^UIdw|3I%VC6^yu+JO zh%HofklAEt+sfr0KtLn$yk6N$#FIJ4$7wpMO7oJ#`>w)1{mRSxB=qcO-IgmSs?og0 zv$a(bAClh1%EB4?R8q4sj!mBSb;qA8#;7)g(VH2@v{sCR(lnxO*Lwi-UloV2<3XcT z3Lc^_xc&=qmfKNXZmA^RDOIhn;|?#_rIS~V2bOO##EYN=Dn*(Ja>Ga4Si;1I2bjQ; zJm2sm3Yu6^(iN!mGx9qI%y({!lL1%J^_8AV_80czITtKX_X*yG2fIG2ex({GD+CvB zSiDSH-HJ$%$7g6zsy+R@C@4j@CipQzF{8`jqe`OB9+0vV2~54>l_*48{{8y~mEYki z*2AeB*x%ipm@8A&05dD!ml;~FOaL%p{hP>mrHBgQb*0d$eg!(b% zKp-1VsVB`a;2gu%`$Y=$mZ|?9uqj$|73*M5WNyCbGX<^kbG_=mJTU7dA-eRq+Ou|i zWu_P3FZY*?bka$tFDu0CXrG9Pt6!YHvS##7q!YK>H2Y=s@8S4shyGw?Q5=?xj18EcYLw zZh)&~vILHM>{9YiyRO0ANmIOs3&j?Ygim^=lj56LO?{?Di2^MB^Zgd>@84e$Jcq=` z(7NA@PaHr+IlRfVZmO}llgGV5LBG7s zX@6cu;MOtiyuM1L9OEsZsIjv?*PcC&yCwOGP=zxDlpqC*nIk5qqQY~DWKPAjfS7OZ z)?fbmawK01SBYG2l)E)EK5Y0s;n8U7EFpWsj@LY043%-Kb$favx-7H_xWWNV! z<;Fd(Jt0=Yg~UIVtx9)!l~{M=dW~>VN;pEwTnNpIu^CZqKD-3#%K}c1`*Uyn&`*rN zofzfknsx7Ra0axp8FfU7R^<*{UCVCMK($47O*Q_b-CPPC=9y9hZVD~5WGUDBenf5j zqnBpt(Wtv*az*x;3auIT(CA8|T@YO5qv6p&_Ok#{W$Ve6-$vuM&^)NtBm+i=z;LRC zr`{2IQJVeT#|Y*sD7AE?h#0wBHgtpRrc3aBE{dOmGL@?cXWE8`NF25y244sO>eh1T(v|)9;bNjO@d{qIpox2 zPZVfviJxHXoCAQqJs@<|JS&Z|o!2^Az;~)_AG;WL{E6qgzaMlrv~}5T z!^EY4QgK=Y)^o^N^i5aUVE9N=ge~$5EwV1>(b(2otXCfuU`O0)#&+c>v?vHk5=*NO z5Y`;sd?G^Tyy~n}?X&-%xQ_o{e)=K7b-VXUf~KLv!anrWsIHsmo8cUGl9X+nM8W7# z>G`?QoBhkX?KDaGbhST=$-sX57klIjYzk!obNGNr&=^N6qm|JG$U?)lDYW`>l4DXv zic>+6nk{Hl2$NO#nnqt%lG*3Zw4z*FMA>N%*w_aq&UO99b9RdN+Lv8{ zZh+VVi?md_{`S4z(`xCr31uK5O8z-jGStx#@G zR0i?I=^^xk`et2dEHNu)P<(Hn89U9e!YJrygfI6*DSP^~EJr}!_o54LgaxTeF;z(J zEG}LB0tZ(HeG!RMy$MF$${wg?-9%;XL20!yzvJe1K&GX?{hnuLx4P8h>L`zChFwCR zMXAT=s9N-cP`&%!G8D2vctb+ZU7UY`v;%FYC^=s7;)!Shhf3aJz}{6HGaT<8V5W_j zIR~pFaz_*fOv=d--ur<1noEbDP6+*qUmq?Rt{Pqk%Tl>D08dC=iq^8 z0pxS=g8RA3P^BLG;5TogBDA=eX`utD0wH;!lG*@8Zk!ZX5roI7IA8Gk{ISI}+}*;W zv5#s_2jaRSoJ}${Hi#8|$6C`KT2~a0-EJ&||J$M3vRvy$1%V-li)SU749(QS1K8I?WBdKkZDCFBP^BLwN0Lw+7s)L3#(2YZRtn-CF|A4!lu#FtjllPX)Y^BdY~SBsX5)1KuA@ z|4S|KdzI9&y7vVU)jd%tw*7;gKKczrWaw0VNexomgbST<$@;e!knMGyFpu_q73Ax2 z2a;R}X{&l~v-A5IOlu9-??vq2o;PGr62TW2CU(OodhH9Yz-n%H}p?;c-X>0sH{ zP4qEXe+sCO^zx`6O4miCV!Hq^GI$?x24DZ+WQQ5K$d1kP(>U!&kC?V}&3a#sv`eg% z+ygF3YTZnjKX%;%*iwWQYKa4?n-`Xbz0xNvWI~BI$@fo4s?=!PmAl1%dPrj!v0TkO zH6|)ES;Z->YV|$v(dHL{xu_KjsJ-gSN4e>0#K*JjYT&W1tvix-ZK|<_m%XDIAi6^L z^Lx5iJ5ttw*JlXRgJuq%@3G4ruh13pPA0nza+-X96_K0v=|D7loS?zlJl48)m#ygH zO8|8~Cd#-xVb)s!r#505uHyt z?XJX9&n%l!AyL)raIbsc|X(b)M?Qs7rZk&m}`<@cd2k#)9zwBp@ zvN>>4qJOy>=Pc9_PboKpN6d(8)Bb+nQZ(cT^r#5VFl$$hcJT~5SEf5 zvt6m!cM#HZ53gk{p5Y$w$II*+p91G*bw0W%JrQmAos0alZ8Y-Sx}7oX z5D_kT&fqTcn6qQ>{UFoI=5y2_J1#XmtzZ8hAiHnCX1DJgxUa$KFAhTYyZEFik;q5S z_{Vk13kkm{VlF0l%CBXgNGp7?-t&C1XFqeVQjhEkM(;18+095v#lyqgcGr6EN{$yo{6+FJ&NiLC zYt!1nqpDzc#1j@N{1!DIdy6{&6$6fqx@#_Aj63bPSbqhY_56Lmra8Do*-=dG4SBI( z8*z7&-lT)(=~Mk#{XEz3Gz@6tO54Vm?!yVyY8hxV)n!N1HqNWYFfKQS)Y1mxEnhKb zyk+lIp=)7}c3#Dlwy_S6a#YzOY+fmIl>zEq8lrNjF*kr1vr~ zH@+RO7$iK=AE(%#yHKCwxjpZ4X}ycaP!sedOf%~wHu#Y&5iM9G&hIax0B^5nL**F9 z$V)i5JI*j*`zc7gd>-IadB+qSkJ)8EZGn#G9vv&5)^T2|&rA7UoS*KX#>_ahZm-7R z8;C!&neFIl>i$&Oh{{>|?bGStxjhnnY5OJEBhJh9?t8#ddg49cL_PhacfhZ~j3F#! zp;2s4bckc8!2YOS{Zt0{x`#nz*GN$;Xc7KrpwTUcN8?iPSjg&`gNDbV2z8jJ8!eID zn6MRhoafX5RiDf^h8)scdZHoBPB_C)y?$g$?~L!ddFCbXX(E029>)vi=3CpHk@I1I z^8+`{IDWeBx-w+NQ)!1ZWUuwUgi~ZQ-%)A&jg{QT2c;>j+WFr+3J-X<&0W4()V+1n z)iG-yFQvNx&#OC}tjUsg@>mnJmBr=nUFVr`HlZ`2vRZSue7ygoZpd-%?I6bvyPxr0 zan)Qyd$SfSW(mdj;?D%K-M6PFBx+>Au}fBi_OlIwix8H6q=2BoQ}bcWd$yy5Hdy z{r-KJ#PMrZG(A0A-#)~OF6?NIJucVQ)iflz&}gkhWG!;3jPcu4Bt_H$-0iswP-Lq2 zap|g{f4xU2Ftw})J(*wf;q2heHg1fA+9`yu|I|)NcH7DkS9z*nq0FSpP1hyX4ArXH zJP0^uNa$6iC+XET7v2dz1it|%bGhARZ&>UVH(CWPyPl3}X3Ce)@km-5v&(Gn4~j=5 zq8cNzT}H)x9y=L-mW-=gThJBeQOhnNCYzM!pFP*57EJ2X_zck-Ipu%FyWv>tb2ue3 zT&wy>S0KQPU%T_BdT&6)w02pys~f#SHI=PNA_IB*_u1^f-tjm?ADO8|RiWTQp?Klw z*G&lL2D0I|3|@1 z%N|zo;66pg3&ACXdW5Vgw~hRkv?%3iEZa8gn&{2qj9$2@t+YVp3r?g((%q_%}q87=B7R&LQ zw<*R&vnaKx%EOl@YEuV78;$P)SaKwV*+Cdsen%eSU+w{m4CCZOCkwa0sL4QU%#@P8fyLo`;4P%6iQ$-WYynBeWa zn5;KP!}kEbOav1G%kglk#c9cke3kHLx@HT1X96-@DLi9$h{b{nS2wRM znK6sjFEY`KYT=_6?Q_Pnp4ykWt7*yioL`kC0uBR++svJ>CIOqJ2_JP;1LV$C`?9T_ zHp#B+%bB~cKqJjjdFD?uePZg*@3Mm*oYEu37;&0NFXGA)MUJt`jUsOwA)!ZgfOAiMg(6-_3|c81P7tN8%jic zC93xOfmCc1Ipoh4$hWEcGHieAQak60Mo*SyKK3M#O$n9ciTj6?s|sQ@-Z4fPc;T}VsvM)c6NV#kkj z4#sd8J6#7~?%;C(zwEjCrr%9|xM8ttyd6V&mx}*Q^FYD^QyvL7r~hx^xup%Ps#emb zNa^3Zo8ae04qz%1E8P;qA2YUe98!fF=Ek3|j+@F_2eW9m|JK+C>^TX{`{WlKYe2mK zlPif7!Z{a$9j@ES&&M*OCj5wz6UFxcw2teH?TI)3=I_5G9VR7frkOUk7zf$V%n^C? zo)d-axb!zXOlfILmk@*B$Ag?9Yo`s9U1K3R@59Odkuf(j$`2^@K#}>j;#LlpHw!6n z4P6FTo6@caO7k4c$%7}86mtw6SFe~|jz%=Ueg_q_q!q?-i`lf3q|Kh{D|i+j*B^L< zo*O3F?H&5v@IW+*P7W9cHD8ksxlI(VR;ZxMI~aB-JR(=ShwMj%`ykEt0I23+ndyhl zNa=YPhIwD6k1g{=QB_gH;%klc_d23ZY z;L<@&;PGkA1{GoK!L!U>w_%8FW-OdVN`Q4!kM8gMz7YD*ZlXbWP^9rkK**2IZ8S2E zd73E)4>I%RHXV`EUcotW6*~2jPGwp7%Kux5jq){ej9S_Rh^b{i90XyN#Xj~Mr&$&e zd7oRI9l3mV%q_YqI2$(&+sSzqs#Z~r5i=&D&xHYi%U zuyz_&?s4@!D#5cw{lefv}T z5Pl?(hNg?xoR_#2V~=u_o_+%nk*~R57`tr!MVsBpcl1rv0En)`-N|vw$;y~7BbtV9 zG&D)nlizr~ryA$RxQw4!PL^ZKf|1v5*g7T?k+%Mi(SINLmo?oy+Ad*VE0Fg9zb>f_ z-2i^(MSpZVq+gu#%GgBt^tc z{bW^k)1GB(P*4k_izMw`yhzT|yzrQ-mtfXjD+j_fD1@A#a`Vp|91$=>h&|{*+CRD& z5AKYvnR*ZMeH0wn5hk*`jFzfOU;cReYsbhT%+W4DmT@{rkgv0!gh7GPJrXpk7(4M# zY&_%8i)@qU)52A)oI+RAIeph98(DuhiV-a>=reitWMoHn@jFI`lx|a>8xAJr&y_f? zTJIo(XxBQN7cRXRlgYy~hHTIMnhrgyE^l0AAI?2AOsk?>!f#M2P4|vpZN&drlDSAZu->81Eku;6FsD#PQTs!^0|8O-me8QTP;s>ChInw=|B&PBT!6_|jLa zzgs*xfBJ69+`5oTwQe@v0yF3RS-XzQSB^%$R>joBl=5Hg&w4Ge!9}1|TV03lMvjx= z01>vsDA$5+T=HIr1L`)7Zr5{up&S9`zQ?xJ-7XvA4OpKkOC1{&?akB_XYrmVdj}>C zN_Lkbdd*&|&wr&b46&4oQ<)OoT3js7(jn9Vg*u8{|FvrL2K6kg&3AtbTj)WwIWW(* z>ZW#%%5_-{AXtYMqs0RGkPT~M=Nn>Od#f}LswDu(_?K<$Pk%MbQg9k$b)I6ZpBa~` z4Az&bMJr(VZ=_JY#A{-!B$-ooh8eZ|OCLkWoK|ELRX%G$@tSlYYTjmSCwHM?dkJ4d1Ze-f-U95>gilz zQ2C#JCgOaq;vyxh&own))`y{z)l`G27K08u>y3Xt+B(!ww;AfBpj55u$1DjR;bnX; zCOR?w={G#OHoCw!Xup`8vC1mS+yRDl3#F;S=a23IBw>PD-`b)h(E@Vk(yUlvar_|} zG6b^292V2hx0Nn~%~+9`Zh695surIjG(60i_Tw@bo6VE=d0`y;7qhgZ_*?AN%2k?Y zFFD52gm%@#%iGy(=X8-oRjq%0@mM(q*2+#j9D7Z#q&|!vQ`WAVzv^K}s{QC8JgHZ8 z!o?ULVF2UY19Dbzdu7TlAuzG1qY8YN=unA2+}sNFeZnRl^lCx6$Jc{FQfrsdwKtFN zY|$k9&d?hvr;L9mm;e19VY8)riEWYUZzb}Olh@*}3|*n%QGQyyy|>^Ms9E{tJOzt& zo-WF9`;BmJl2ZN1Wl7PCFOcMpuYIcsX5I2rKb-@;+Tu=akAdPntteV@N3n&akPIj!mTHUTJuM(uzrARmB%0xhFv1R_g z{aSOAJd0BQbzOdAsA~2yYDY)UHYz=4<-~|#UQDNS=5d{e7#o85Y@$I|!g(O? z-72N#1Gya9$2PU9wWh=@WS9|n-S;6prOWYr)SP>c;Uq+W=k^{jTlvcfPL=pfxms{~ zc#cCh+CV1efjf>=u+{3X-ycqpd zkFBK z3<;;gaIBE?QVq9jIncdK@n@@M$x{Ov#93&6k7Mm0YGAZ8hiXea!eh-PB_(%}7<3K5 zWCru_b(fuT^uB@P=MuqsfWzQTol{Eo2r^fY*`Geq5=xA;x6gE~+fwbA-)4K|q?=fO z*tF~*x9qd*A8E&8<9=XEL{>C7eZk^B#D8_Etbe#^Yg1Dhls)4IPCYe?$Xrjp($<94 z7`?~L!TZ&_Wu1HcTdH5lo#ER3(8_m`eheEyvdUZ|1q3qvvjFS$d^m(vw#h9 zhh{pnOHWJHj|swOxmDJ?y(i4XP(#qgI2FxdUh4jwv$67JW7EgfzO(JCX#TL~_Z5zJ zVSTqRc<%u+HKooCkDQtUF`m~-`2QxnIAY9_K?sus1tRLL_iz=r#Mg6be!0W*(j}`F zgF-S_l|)P0jJgJ35gX-vXi9^uo~;{+_g~}=g68PrOeuUFg{wyS`QKgnrAgfbNLQE3 zkHhw1=ubj8Ge0LTp@44&&T~R;bSPe*QK`1cAl4-0tBA;B#ZUc#095qa4!@JwJXT}+AZ44HC_X>m$Yn6w7aw? zZ*8;$sIVBqiWb5TG;%0Xyjc@{TX!YI-W1-EFlvjntL*p&TS#d z@M?h7i!97bOLw#@5;NF-5@}yv-L*QXy963Su}W*+;-%At<`;}KtNZD;9)7}YnUdVf zkNW>^a>=9j`~xs2s;x+$0~)*sIAf`~@m^`!CE~n{xC=cOhqa@k$|@qiKvUkBV_itp zX(&@Lul+3F$XgYXl4r1YbnBvTlB$}uN%Z0hkNzw+yFejud_o+=#`w~YOE`MH5uE*ENej`777e$2pvsr5n&O~*e z?>09?!oZ8k`L4S+3u_NWX%BK4M({>Mdd8ybFi->h?%-e$En57h!5ikBPBpditM&X; zvz~slXFQM&Uh%GVp5MTOo1{=1PhwsbdalyMFS&Qojw*Fwd$yp^4x5T%U76ZLpZ%ps zd2sJb)SwV&l>iYkZNL4v5iUb?O>G=NlwM*J;r+SF@x6v3*Cyn=ynifWJ%Y2{N++ge zF`B^PLz5R#Y)i|9Wp~9r;89C*qu*bhl!r)9HOC7Is)JI~sE2*%uaNF5_up z2oL9gMl1uERM$2gXQJEGTBIx|-0-(M7L7*+u!gu2fxr|bI9@0mnp~Wh7L#xH=48V| zGL&+aN%5%wu)lfzI(jBsgb;`I9e3+_&1tf2{keQQuF>wCgQ$0ca?d)nHq>do6ZM*& zc?sS=Gd!JZlunCqWu$mAbs_w`Yz|aymExDuo^&)Ml-jk7j?rgCYH%4lQ9YGn@ZO$* zs3jQgmpVNV(NR1Wwbu08?r#!%iz4@cF{xtnprF!Xx75y=zv}kxLLYzd{W{hrN{LCO zyHZHzcvAJlt2B&~Bts%(ct^X_HMQeKT((g>@l6c2#o}V7Nmaf9tX&2qkpq|VTRI4Y zZEHsL>}VP3m9qo<6^naaM}i%pMr)aitE z*EZ1t6tzy^4iM>9vqZ#=z^y^Ld$ShT?lPBiBVuJgMLZnhwrOy{PcQN zoB}9nD<1BAV7tm66&R~Lk9?XJJrvEmRk@w-E|EZ>gshJbqSeAkC*`Zsv1LP0it~(b zC2a~*1HNJBj>I+W2nXpdg6)F^_dPm8%R9xDcw)<<>e!6qZVtACeIkhy2@K8t=KL)b zF1}3ehID#}plIO~fkP2U1yMo(MT{wZ8cuvwpFEd8 z9_UfHgv1Ahuau!{@omzC+zCUqB@X=B*?8}WzoY&3=OHGr_a9V`?LvtyR;pdXF~`Ju z;g@QMx;g_!1YE+tErB6Fs-b*E?W&8SJfEAx2FPW#Xl#7S_%Lm5Jv zG|z;-a}DV#Gun+j-)=7(!GdIBapb)J@by?NVXao_ajz-6E`q4LK#Bb^S_UgF&j{a& z!Qud|Gs-L8GzneSd0+>E_kfuT+)k)RLj1TdnxJ?j&mccRChCE#{aC>*#NEQhG612DeYMG6X#g*aD2K*)K zjTFw0lO_bw36qogb=um&s;BJ1SLRs?e<0b>E9=h)wrD7=r%Tf@U*`zqbsU+rbGrRj zZsWHw$*ymS;MBKasDwfujA?$(P`H2y+T2?WHF0tNy+@hub-mSbqniCR9w9kVR>})GP_d&?xB0T_I(@dOwVr^K5*vyV z=# zLi_#j2`eOjr|$?I{R=nG9*m}-=`W3eRRUCO)8V~Ez%rM3D%-g?8XUt+pn@Z9BDO|p5SP$c%FMe3VLZnWx#mw>TX*3H2pfAWNGTs2?ZRx2h2H4Dx2M5dkeVXo;L{HF$fPo zNG>Xkn%zqsmo}+QVV%g>%3in?OAkh31TkGPXP66&W{OfMdz?4I>F4E-Jy8uiVvaR8 zwz_)P=C{}>xf^pN1*)8{j=2tv-F*Q%&xyvg=DgqKmTQ;b8|@SQv7YH3%V(`!@6Hd< zA|fsk9r}H_;=R-8<2aV(H2qRLEa6UXl_c!Oh$!trSHUy8T&JO0ol_V;;V1Bw;`KXm zJQFA)it}~O*GCD~Huv)w+Q{@$TY5;zsoFc-RO=@R=LuBw{+%VZ7?})pJXGA7ICty* z)O`v1lll%t;qO62bK>?H*Uh24(V*YCiI(IS%jLuFk5=1Sd)WtLa%f_S0vM-tAmIFu zoE8S36Qi($Z&ND-z22^x1*iO4w`ib?@2`%f8BEMF8_j=W4CVT9<|y$i^?NoVmhU{F zQ#t$GalOd5K|nL7nhxf}30LbG->ke2B-;eRZ+piJCbPo}$mN`N$BM_sXqCCjO37r0 zF(vC}G$Mj7_dL?XE{#~XCLCi@G7|@bbGkJ&a{ZAy*Gl#+N+w*Zk_z5WJOq11y}l>H zcoz>%sSA{bl{PkT?(80m>@(nZA2oVY2CxtG@_a2Fib7l^W9FQg9-g8zqbvq65Jct= z!z8%nnI~O9T~=3nRWYZ@s>tp}*wx3pOZ=Pqed(oRi0;&cPK8p6aDOqPyLRszWekw9 z?c)Vs!ybWuDgIP*z}8`n8Q8C1e+d{!&ePfXo`YH8Z{G!w>7!xVJnsO3;%hN*Rkj%HIDf!Q{ogB-$?L zynk>X=EHqz@zHg!%wQ?9E=tKF^B3tl0>)eBL_}<`y{2AQ1P4ackOQ&OkFE_jafbuL zbuBl%tRzgv+#;i-%|{Dd;m=`CteeehNnW>BnZurhsIt@~KLPa5ViXWKUlqQ1u{5R|da;j@ zIAGsGP@WQq5u4<{BL()I&WbA8P2vY@j&wh8hzB;RrsTP6}C-X5A@r|81N2ozqLPfdjoAxaSuHbx?xs~s3i%A8*ug;{gBg@G%#=* zqiLJQX!&eQ8~9XO%(uzEr+xEt6Xx#Vz*d~#6G!Q9uycl5v8ljkrnJC|&M zEQ+zTtln%V>moR4b}UX^J^^X$V^!J%E21}!Ho299c2*5N)Vtl!br^l{*+g`{Z-a=y z8PxcHTM20XMf@7Nak{g8xQk`qndC+aJA{je|3VdS<9ED-pAMdlxdGhfH+~Y)J@&L z^(~Fx*KkZmPl!n!a`b;>C-$(LOc@zpT;94E5f(0It7)%YNI0JB&{bGCGcSifD)Kq! z2;mW3I~RA3=qE3r|KF%E|12p>Nyv?yty3%7Lm^S{xXd`~^^j1JDO6qYPI+V#rFJrHRmUOuVcfo=xqW-0i|USPgSl%lzg}${!MOtf6_!XJwd$FM5=u;MB{mLl-ulv}FMFXPCFJ;}e$+Q7K;5MjK9U^=oN}xbXU>phd5|BjeVAURr$4D+SS?{JYhygLk|mYWq&+zW zy(hXe#*Vu5{*ro3f|Jv#vQlG$VvDC^vRU{wO|!kPz-9}*8GZ9e>QI4M|9gwd%-oq+OA2qpD?syZ88RbGb^8E~Jso7wtmF$blGGd)17Hrb9d-yR@ogl)@c}1X zu2(C4XiD=qimv6I38gY_nGrg0e`w1haAf2`7K3&)zaBLBQ%(g5;{A>*LzoKS z)BW_WV9@k_f_I_U3_=T_b$Q4eVmn0#BF9G|0f>nz=JH5Lf@~6C6j9SlQo1MN zh4KF6Kx&Y)b8j7}GsOXte3vYK+rW8c=-(2Wibyvohdru$NMf`-m|Tt0w%x{f-6Njc z+2PGt7wAPpf7Q^M~!n*miyk(D_^rPuA$BGwAd)#^)9XJl`N zm-OqmMsnW?F=Q(?smk+gMwNC^K^W$>FCuYWggT-Qo{D#Z$rwQ+fhbN7xTQM5Az;a0y#nppMQl)9I;+CbW2pt3>{NX z`UL%7A=`}um!QX<% -#include -#include -#include "pair_lj_sf.h" -#include "atom.h" -#include "comm.h" -#include "force.h" -#include "neigh_list.h" -#include "memory.h" -#include "error.h" - -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -PairLJShiftedForce::PairLJShiftedForce(LAMMPS *lmp) : Pair(lmp) -{ - respa_enable = 0; -} - -/* ---------------------------------------------------------------------- */ - -PairLJShiftedForce::~PairLJShiftedForce() -{ - if (allocated) { - memory->destroy(setflag); - memory->destroy(cutsq); - - memory->destroy(cut); - memory->destroy(epsilon); - memory->destroy(sigma); - memory->destroy(lj1); - memory->destroy(lj2); - memory->destroy(lj3); - memory->destroy(lj4); - memory->destroy(foffset); - memory->destroy(offset); - } -} - -/* ---------------------------------------------------------------------- */ - -void PairLJShiftedForce::compute(int eflag, int vflag) -{ - int i,j,ii,jj,inum,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double rsq,r2inv,r6inv,forcelj,factor_lj; - double r,t; - int *ilist,*jlist,*numneigh,**firstneigh; - - evdwl = 0.0; - if (eflag || vflag) ev_setup(eflag,vflag); - else evflag = vflag_fdotr = 0; - - double **x = atom->x; - double **f = atom->f; - int *type = atom->type; - int nlocal = atom->nlocal; - double *special_lj = force->special_lj; - int newton_pair = force->newton_pair; - - inum = list->inum; - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // loop over neighbors of my atoms - - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - xtmp = x[i][0]; - ytmp = x[i][1]; - ztmp = x[i][2]; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_lj = special_lj[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j][0]; - dely = ytmp - x[j][1]; - delz = ztmp - x[j][2]; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; - r6inv = r2inv*r2inv*r2inv; - r = sqrt(rsq); - t = r/cut[itype][jtype]; - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) - - t*foffset[itype][jtype]; - fpair = factor_lj*forcelj*r2inv; - - f[i][0] += delx*fpair; - f[i][1] += dely*fpair; - f[i][2] += delz*fpair; - if (newton_pair || j < nlocal) { - f[j][0] -= delx*fpair; - f[j][1] -= dely*fpair; - f[j][2] -= delz*fpair; - } - - if (eflag) { - evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) + - (t-1.0)*foffset[itype][jtype] - offset[itype][jtype]; - evdwl *= factor_lj; - } - - if (evflag) ev_tally(i,j,nlocal,newton_pair, - evdwl,0.0,fpair,delx,dely,delz); - } - } - } - - if (vflag_fdotr) virial_fdotr_compute(); -} - -/* ---------------------------------------------------------------------- - allocate all arrays -------------------------------------------------------------------------- */ - -void PairLJShiftedForce::allocate() -{ - allocated = 1; - int n = atom->ntypes; - - memory->create(setflag,n+1,n+1,"pair:setflag"); - for (int i = 1; i <= n; i++) - for (int j = i; j <= n; j++) - setflag[i][j] = 0; - - memory->create(cutsq,n+1,n+1,"pair:cutsq"); - - memory->create(cut,n+1,n+1,"pair:cut"); - memory->create(epsilon,n+1,n+1,"pair:epsilon"); - memory->create(sigma,n+1,n+1,"pair:sigma"); - memory->create(lj1,n+1,n+1,"pair:lj1"); - memory->create(lj2,n+1,n+1,"pair:lj2"); - memory->create(lj3,n+1,n+1,"pair:lj3"); - memory->create(lj4,n+1,n+1,"pair:lj4"); - memory->create(foffset,n+1,n+1,"pair:foffset"); - memory->create(offset,n+1,n+1,"pair:offset"); -} - -/* ---------------------------------------------------------------------- - global settings -------------------------------------------------------------------------- */ - -void PairLJShiftedForce::settings(int narg, char **arg) -{ - if (narg != 1) error->all(FLERR,"Illegal pair_style command"); - - cut_global = force->numeric(FLERR,arg[0]); - - if (cut_global <= 0.0) - error->all(FLERR,"Illegal pair_style command"); - - // reset cutoffs that have been explicitly set - - if (allocated) { - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) - if (setflag[i][j]) cut[i][j] = cut_global; - } -} - -/* ---------------------------------------------------------------------- - set coeffs for one or more type pairs -------------------------------------------------------------------------- */ - -void PairLJShiftedForce::coeff(int narg, char **arg) -{ - if (narg < 4 || narg > 5) - error->all(FLERR,"Incorrect args for pair coefficients"); - if (!allocated) allocate(); - - int ilo,ihi,jlo,jhi; - force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); - force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); - - double epsilon_one = force->numeric(FLERR,arg[2]); - double sigma_one = force->numeric(FLERR,arg[3]); - - double cut_one = cut_global; - if (narg == 5) cut_one = force->numeric(FLERR,arg[4]); - - if (cut_one <= 0.0) - error->all(FLERR,"Incorrect args for pair coefficients"); - - int count = 0; - for (int i = ilo; i <= ihi; i++) { - for (int j = MAX(jlo,i); j <= jhi; j++) { - epsilon[i][j] = epsilon_one; - sigma[i][j] = sigma_one; - cut[i][j] = cut_one; - setflag[i][j] = 1; - count++; - } - } - - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); -} - -/* ---------------------------------------------------------------------- - init for one type pair i,j and corresponding j,i -------------------------------------------------------------------------- */ - -double PairLJShiftedForce::init_one(int i, int j) -{ - if (setflag[i][j] == 0) { - epsilon[i][j] = mix_energy(epsilon[i][i],epsilon[j][j], - sigma[i][i],sigma[j][j]); - sigma[i][j] = mix_distance(sigma[i][i],sigma[j][j]); - cut[i][j] = mix_distance(cut[i][i],cut[j][j]); - } - - lj1[i][j] = 48.0 * epsilon[i][j] * pow(sigma[i][j],12.0); - lj2[i][j] = 24.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - lj3[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],12.0); - lj4[i][j] = 4.0 * epsilon[i][j] * pow(sigma[i][j],6.0); - - double ratio = sigma[i][j] / cut[i][j]; - foffset[i][j] = 4.0 * epsilon[i][j] * (12.0 * pow(ratio,12.0) - - 6.0 * pow(ratio,6.0)); - offset[i][j] = 4.0 * epsilon[i][j] * (pow(ratio,12.0) - pow(ratio,6.0)); - - cut[j][i] = cut[i][j]; - lj1[j][i] = lj1[i][j]; - lj2[j][i] = lj2[i][j]; - lj3[j][i] = lj3[i][j]; - lj4[j][i] = lj4[i][j]; - foffset[j][i] = foffset[i][j]; - offset[j][i] = offset[i][j]; - - return cut[i][j]; -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairLJShiftedForce::write_restart(FILE *fp) -{ - write_restart_settings(fp); - - int i,j; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) { - fwrite(&setflag[i][j],sizeof(int),1,fp); - if (setflag[i][j]) { - fwrite(&epsilon[i][j],sizeof(double),1,fp); - fwrite(&sigma[i][j],sizeof(double),1,fp); - fwrite(&cut[i][j],sizeof(double),1,fp); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairLJShiftedForce::read_restart(FILE *fp) -{ - read_restart_settings(fp); - allocate(); - - int i,j; - int me = comm->me; - for (i = 1; i <= atom->ntypes; i++) - for (j = i; j <= atom->ntypes; j++) { - if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); - MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); - if (setflag[i][j]) { - if (me == 0) { - fread(&epsilon[i][j],sizeof(double),1,fp); - fread(&sigma[i][j],sizeof(double),1,fp); - fread(&cut[i][j],sizeof(double),1,fp); - } - MPI_Bcast(&epsilon[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&sigma[i][j],1,MPI_DOUBLE,0,world); - MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); - } - } -} - -/* ---------------------------------------------------------------------- - proc 0 writes to restart file -------------------------------------------------------------------------- */ - -void PairLJShiftedForce::write_restart_settings(FILE *fp) -{ - fwrite(&cut_global,sizeof(double),1,fp); - fwrite(&offset_flag,sizeof(int),1,fp); - fwrite(&mix_flag,sizeof(int),1,fp); -} - -/* ---------------------------------------------------------------------- - proc 0 reads from restart file, bcasts -------------------------------------------------------------------------- */ - -void PairLJShiftedForce::read_restart_settings(FILE *fp) -{ - int me = comm->me; - if (me == 0) { - fread(&cut_global,sizeof(double),1,fp); - fread(&offset_flag,sizeof(int),1,fp); - fread(&mix_flag,sizeof(int),1,fp); - } - MPI_Bcast(&cut_global,1,MPI_DOUBLE,0,world); - MPI_Bcast(&offset_flag,1,MPI_INT,0,world); - MPI_Bcast(&mix_flag,1,MPI_INT,0,world); -} - -/* ---------------------------------------------------------------------- */ - -double PairLJShiftedForce::single(int i, int j, int itype, int jtype, double rsq, - double factor_coul, double factor_lj, - double &fforce) -{ - double r2inv,r6inv,forcelj,philj,r,t; - - r2inv = 1.0/rsq; - r6inv = r2inv*r2inv*r2inv; - r = sqrt(rsq); - t = r/cut[itype][jtype]; - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) - - t*foffset[itype][jtype]; - fforce = factor_lj*forcelj*r2inv; - - philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) + - (t-1.0)*foffset[itype][jtype] - offset[itype][jtype]; - return factor_lj*philj; -} diff --git a/src/USER-MISC/pair_lj_sf.h b/src/USER-MISC/pair_lj_sf.h deleted file mode 100644 index 1a4106b782..0000000000 --- a/src/USER-MISC/pair_lj_sf.h +++ /dev/null @@ -1,53 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -#ifdef PAIR_CLASS - -PairStyle(lj/sf,PairLJShiftedForce) - -#else - -#ifndef LMP_PAIR_LJ_SF_H -#define LMP_PAIR_LJ_SF_H - -#include "pair.h" - -namespace LAMMPS_NS { - -class PairLJShiftedForce : public Pair { - public: - PairLJShiftedForce(class LAMMPS *); - virtual ~PairLJShiftedForce(); - virtual void compute(int, int); - void settings(int, char **); - void coeff(int, char **); - double init_one(int, int); - void write_restart(FILE *); - void read_restart(FILE *); - void write_restart_settings(FILE *); - void read_restart_settings(FILE *); - double single(int, int, int, int, double, double, double, double &); - - protected: - double cut_global; - double **cut; - double **epsilon,**sigma; - double **lj1,**lj2,**lj3,**lj4,**foffset,**offset; - - void allocate(); -}; - -} - -#endif -#endif diff --git a/src/USER-OMP/pair_lj_sf_omp.cpp b/src/USER-OMP/pair_lj_sf_omp.cpp deleted file mode 100644 index bd9d5220b8..0000000000 --- a/src/USER-OMP/pair_lj_sf_omp.cpp +++ /dev/null @@ -1,164 +0,0 @@ -/* ---------------------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - This software is distributed under the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Axel Kohlmeyer (Temple U) -------------------------------------------------------------------------- */ - -#include -#include "pair_lj_sf_omp.h" -#include "atom.h" -#include "comm.h" -#include "force.h" -#include "neighbor.h" -#include "neigh_list.h" - -#include "suffix.h" -using namespace LAMMPS_NS; - -/* ---------------------------------------------------------------------- */ - -PairLJShiftedForceOMP::PairLJShiftedForceOMP(LAMMPS *lmp) : - PairLJShiftedForce(lmp), ThrOMP(lmp, THR_PAIR) -{ - suffix_flag |= Suffix::OMP; - respa_enable = 0; -} - -/* ---------------------------------------------------------------------- */ - -void PairLJShiftedForceOMP::compute(int eflag, int vflag) -{ - if (eflag || vflag) { - ev_setup(eflag,vflag); - } else evflag = vflag_fdotr = 0; - - const int nall = atom->nlocal + atom->nghost; - const int nthreads = comm->nthreads; - const int inum = list->inum; - -#if defined(_OPENMP) -#pragma omp parallel default(none) shared(eflag,vflag) -#endif - { - int ifrom, ito, tid; - - loop_setup_thr(ifrom, ito, tid, inum, nthreads); - ThrData *thr = fix->get_thr(tid); - thr->timer(Timer::START); - ev_setup_thr(eflag, vflag, nall, eatom, vatom, thr); - - if (evflag) { - if (eflag) { - if (force->newton_pair) eval<1,1,1>(ifrom, ito, thr); - else eval<1,1,0>(ifrom, ito, thr); - } else { - if (force->newton_pair) eval<1,0,1>(ifrom, ito, thr); - else eval<1,0,0>(ifrom, ito, thr); - } - } else { - if (force->newton_pair) eval<0,0,1>(ifrom, ito, thr); - else eval<0,0,0>(ifrom, ito, thr); - } - - thr->timer(Timer::PAIR); - reduce_thr(this, eflag, vflag, thr); - } // end of omp parallel region -} - -template -void PairLJShiftedForceOMP::eval(int iifrom, int iito, ThrData * const thr) -{ - int i,j,ii,jj,jnum,itype,jtype; - double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair; - double t,rsq,r2inv,r6inv,forcelj,factor_lj; - int *ilist,*jlist,*numneigh,**firstneigh; - - evdwl = 0.0; - - const dbl3_t * _noalias const x = (dbl3_t *) atom->x[0]; - dbl3_t * _noalias const f = (dbl3_t *) thr->get_f()[0]; - const int * _noalias const type = atom->type; - const int nlocal = atom->nlocal; - const double * _noalias const special_lj = force->special_lj; - double fxtmp,fytmp,fztmp; - - ilist = list->ilist; - numneigh = list->numneigh; - firstneigh = list->firstneigh; - - // loop over neighbors of my atoms - - for (ii = iifrom; ii < iito; ++ii) { - - i = ilist[ii]; - xtmp = x[i].x; - ytmp = x[i].y; - ztmp = x[i].z; - itype = type[i]; - jlist = firstneigh[i]; - jnum = numneigh[i]; - fxtmp=fytmp=fztmp=0.0; - - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - factor_lj = special_lj[sbmask(j)]; - j &= NEIGHMASK; - - delx = xtmp - x[j].x; - dely = ytmp - x[j].y; - delz = ztmp - x[j].z; - rsq = delx*delx + dely*dely + delz*delz; - jtype = type[j]; - - if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; - r6inv = r2inv*r2inv*r2inv; - t = sqrt(rsq)/cut[itype][jtype]; - - forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]) - - t*foffset[itype][jtype]; - - fpair = factor_lj*forcelj*r2inv; - - fxtmp += delx*fpair; - fytmp += dely*fpair; - fztmp += delz*fpair; - if (NEWTON_PAIR || j < nlocal) { - f[j].x -= delx*fpair; - f[j].y -= dely*fpair; - f[j].z -= delz*fpair; - } - - if (EFLAG) { - evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]) + - (t-1.0)*foffset[itype][jtype] - offset[itype][jtype]; - evdwl *= factor_lj; - } - - if (EVFLAG) ev_tally_thr(this,i,j,nlocal,NEWTON_PAIR, - evdwl,0.0,fpair,delx,dely,delz,thr); - } - } - f[i].x += fxtmp; - f[i].y += fytmp; - f[i].z += fztmp; - } -} - -/* ---------------------------------------------------------------------- */ - -double PairLJShiftedForceOMP::memory_usage() -{ - double bytes = memory_usage_thr(); - bytes += PairLJShiftedForce::memory_usage(); - - return bytes; -} diff --git a/src/USER-OMP/pair_lj_sf_omp.h b/src/USER-OMP/pair_lj_sf_omp.h deleted file mode 100644 index 92db973b3d..0000000000 --- a/src/USER-OMP/pair_lj_sf_omp.h +++ /dev/null @@ -1,48 +0,0 @@ -/* -*- c++ -*- ---------------------------------------------------------- - LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator - http://lammps.sandia.gov, Sandia National Laboratories - Steve Plimpton, sjplimp@sandia.gov - - Copyright (2003) Sandia Corporation. Under the terms of Contract - DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains - certain rights in this software. This software is distributed under - the GNU General Public License. - - See the README file in the top-level LAMMPS directory. -------------------------------------------------------------------------- */ - -/* ---------------------------------------------------------------------- - Contributing author: Axel Kohlmeyer (Temple U) -------------------------------------------------------------------------- */ - -#ifdef PAIR_CLASS - -PairStyle(lj/sf/omp,PairLJShiftedForceOMP) - -#else - -#ifndef LMP_PAIR_LJ_SF_OMP_H -#define LMP_PAIR_LJ_SF_OMP_H - -#include "pair_lj_sf.h" -#include "thr_omp.h" - -namespace LAMMPS_NS { - -class PairLJShiftedForceOMP : public PairLJShiftedForce, public ThrOMP { - - public: - PairLJShiftedForceOMP(class LAMMPS *); - - virtual void compute(int, int); - virtual double memory_usage(); - - private: - template - void eval(int ifrom, int ito, ThrData * const thr); -}; - -} - -#endif -#endif diff --git a/src/USER-OMP/pair_lj_smooth_linear_omp.h b/src/USER-OMP/pair_lj_smooth_linear_omp.h index 940c0ea707..874e42eb9f 100644 --- a/src/USER-OMP/pair_lj_smooth_linear_omp.h +++ b/src/USER-OMP/pair_lj_smooth_linear_omp.h @@ -18,6 +18,7 @@ #ifdef PAIR_CLASS PairStyle(lj/smooth/linear/omp,PairLJSmoothLinearOMP) +PairStyle(lj/sf/omp,PairLJSmoothLinearOMP) #else diff --git a/src/pair_lj_smooth_linear.h b/src/pair_lj_smooth_linear.h index 0e3376b789..c18c442a18 100644 --- a/src/pair_lj_smooth_linear.h +++ b/src/pair_lj_smooth_linear.h @@ -14,6 +14,7 @@ #ifdef PAIR_CLASS PairStyle(lj/smooth/linear,PairLJSmoothLinear) +PairStyle(lj/sf,PairLJSmoothLinear) #else From 04ebd81ac5d01c4ed615856e52f1156bd4820699 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 6 Jun 2017 17:26:18 -0400 Subject: [PATCH 13/22] minor whitespace cleanup --- src/write_restart.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/write_restart.cpp b/src/write_restart.cpp index 77e2cb05d9..ad6c756558 100644 --- a/src/write_restart.cpp +++ b/src/write_restart.cpp @@ -185,7 +185,7 @@ void WriteRestart::multiproc_options(int multiproc_caller, int mpiioflag_caller, if (strcmp(arg[iarg],"fileper") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal write_restart command"); if (!multiproc) - error->all(FLERR,"Cannot use write_restart fileper " + error->all(FLERR,"Cannot use write_restart fileper " "without % in restart file name"); int nper = force->inumeric(FLERR,arg[iarg+1]); if (nper <= 0) error->all(FLERR,"Illegal write_restart command"); @@ -203,7 +203,7 @@ void WriteRestart::multiproc_options(int multiproc_caller, int mpiioflag_caller, } else if (strcmp(arg[iarg],"nfile") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal write_restart command"); if (!multiproc) - error->all(FLERR,"Cannot use write_restart nfile " + error->all(FLERR,"Cannot use write_restart nfile " "without % in restart file name"); int nfile = force->inumeric(FLERR,arg[iarg+1]); if (nfile <= 0) error->all(FLERR,"Illegal write_restart command"); From 1f9504c54622f4b1c5eb2a8bd6526a1c7a475583 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 6 Jun 2017 17:31:45 -0400 Subject: [PATCH 14/22] some more bookkeeping updates triggered by the lj/sf style removal --- doc/src/lammps.book | 1 - doc/src/pairs.txt | 1 - src/.gitignore | 2 -- src/USER-MISC/README | 1 - 4 files changed, 5 deletions(-) diff --git a/doc/src/lammps.book b/doc/src/lammps.book index 1769f29825..69c215a2b9 100644 --- a/doc/src/lammps.book +++ b/doc/src/lammps.book @@ -446,7 +446,6 @@ pair_lj96.html pair_lj_cubic.html pair_lj_expand.html pair_lj_long.html -pair_lj_sf.html pair_lj_smooth.html pair_lj_smooth_linear.html pair_lj_soft.html diff --git a/doc/src/pairs.txt b/doc/src/pairs.txt index 538e2a7268..2c1b20f4d3 100644 --- a/doc/src/pairs.txt +++ b/doc/src/pairs.txt @@ -49,7 +49,6 @@ Pair Styles :h1 pair_lj_cubic pair_lj_expand pair_lj_long - pair_lj_sf pair_lj_smooth pair_lj_smooth_linear pair_lj_soft diff --git a/src/.gitignore b/src/.gitignore index 0cddfa6951..e8649c59b0 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -740,8 +740,6 @@ /pair_lj_sdk_coul_long.h /pair_lj_sdk_coul_msm.cpp /pair_lj_sdk_coul_msm.h -/pair_lj_sf.cpp -/pair_lj_sf.h /pair_lj_sf_dipole_sf.cpp /pair_lj_sf_dipole_sf.h /pair_lubricateU.cpp diff --git a/src/USER-MISC/README b/src/USER-MISC/README index cacee41e0c..3b69fd1dec 100644 --- a/src/USER-MISC/README +++ b/src/USER-MISC/README @@ -63,7 +63,6 @@ pair_style gauss/cut, Axel Kohlmeyer, akohlmey at gmail.com, 1 Dec 11 pair_style lennard/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15 pair_style list, Axel Kohlmeyer (Temple U), akohlmey at gmail.com, 1 Jun 13 pair_style lj/mdf, Paolo Raiteri, p.raiteri at curtin.edu.au, 2 Dec 15 -pair_style lj/sf, Laurent Joly (U Lyon), ljoly.ulyon at gmail.com, 8 Aug 11 pair_style kolmogorov/crespi/z, Jaap Kroes (Radboud U), jaapkroes at gmail dot com, 28 Feb 17 pair_style meam/spline, Alexander Stukowski (LLNL), alex at stukowski.com, 1 Feb 12 pair_style meam/sw/spline, Robert Rudd (LLNL), robert.rudd at llnl.gov, 1 Oct 12 From a2edef7c9c3e8b10f72f0d0031129b15392fef39 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 7 Jun 2017 00:23:53 -0400 Subject: [PATCH 15/22] local variable fp in pair style eam/cd was shadowing class member. renamed local variable to fptr --- src/USER-MISC/pair_cdeam.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/USER-MISC/pair_cdeam.cpp b/src/USER-MISC/pair_cdeam.cpp index 91ef598de8..b5607012ce 100644 --- a/src/USER-MISC/pair_cdeam.cpp +++ b/src/USER-MISC/pair_cdeam.cpp @@ -456,11 +456,11 @@ void PairCDEAM::read_h_coeff(char *filename) { if(comm->me == 0) { // Open potential file - FILE *fp; + FILE *fptr; char line[MAXLINE]; char nextline[MAXLINE]; - fp = force->open_potential(filename); - if (fp == NULL) { + fptr = force->open_potential(filename); + if (fptr == NULL) { char str[128]; sprintf(str,"Cannot open EAM potential file %s", filename); error->one(FLERR,str); @@ -468,7 +468,7 @@ void PairCDEAM::read_h_coeff(char *filename) // h coefficients are stored at the end of the file. // Skip to last line of file. - while(fgets(nextline, MAXLINE, fp) != NULL) { + while(fgets(nextline, MAXLINE, fptr) != NULL) { strcpy(line, nextline); } char* ptr = strtok(line, " \t\n\r\f"); @@ -483,7 +483,7 @@ void PairCDEAM::read_h_coeff(char *filename) error->one(FLERR,"Failed to read h(x) function coefficients from EAM file."); // Close the potential file. - fclose(fp); + fclose(fptr); } MPI_Bcast(&nhcoeff, 1, MPI_INT, 0, world); From 31a734b03d620c9aef0ecacb53f2d251ca13e55e Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 7 Jun 2017 17:10:33 -0400 Subject: [PATCH 16/22] sbmask function should be flagged as const indicating no side effects --- src/compute.h | 2 +- src/create_bonds.h | 2 +- src/delete_atoms.h | 2 +- src/pair.h | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/compute.h b/src/compute.h index 7f12cd97e2..279f71829b 100644 --- a/src/compute.h +++ b/src/compute.h @@ -152,7 +152,7 @@ class Compute : protected Pointers { double **vbiasall; // stored velocity bias for all atoms int maxbias; // size of vbiasall array - inline int sbmask(int j) { + inline int sbmask(int j) const { return j >> SBBITS & 3; } diff --git a/src/create_bonds.h b/src/create_bonds.h index 2936506b3f..24b1596e37 100644 --- a/src/create_bonds.h +++ b/src/create_bonds.h @@ -30,7 +30,7 @@ class CreateBonds : protected Pointers { void command(int, char **); private: - inline int sbmask(int j) { + inline int sbmask(int j) const { return j >> SBBITS & 3; } }; diff --git a/src/delete_atoms.h b/src/delete_atoms.h index 62ba47d715..72cf44285f 100644 --- a/src/delete_atoms.h +++ b/src/delete_atoms.h @@ -45,7 +45,7 @@ class DeleteAtoms : protected Pointers { void recount_topology(); void options(int, char **); - inline int sbmask(int j) { + inline int sbmask(int j) const { return j >> SBBITS & 3; } diff --git a/src/pair.h b/src/pair.h index dd859e5f2a..b57004d965 100644 --- a/src/pair.h +++ b/src/pair.h @@ -245,7 +245,7 @@ class Pair : protected Pointers { ubuf(int arg) : i(arg) {} }; - inline int sbmask(int j) { + inline int sbmask(int j) const { return j >> SBBITS & 3; } }; From 0ecdb998856b2ca3d565e3fb0e973f85b52240f0 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 8 Jun 2017 13:50:17 -0400 Subject: [PATCH 17/22] fix uninitialized data access as reported by @martok in #174 --- lib/meam/meam_dens_final.F | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/lib/meam/meam_dens_final.F b/lib/meam/meam_dens_final.F index 92195dcaf4..cdc47067e2 100644 --- a/lib/meam/meam_dens_final.F +++ b/lib/meam/meam_dens_final.F @@ -98,13 +98,14 @@ c Complete the calculation of density call G_gam(Gamma(i),ibar_meam(elti), $ gsmooth_factor,G,errorflag) if (errorflag.ne.0) return + call get_shpfcn(shp,lattce_meam(elti,elti)) if (ibar_meam(elti).le.0) then Gbar = 1.d0 + dGbar = 0.d0 else - call get_shpfcn(shp,lattce_meam(elti,elti)) if (mix_ref_t.eq.1) then - gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2) - $ +t_ave(3,i)*shpi(3))/(Z*Z) + gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2) + $ +t_ave(3,i)*shp(3))/(Z*Z) else gam = (t1_meam(elti)*shp(1)+t2_meam(elti)*shp(2) $ +t3_meam(elti)*shp(3))/(Z*Z) @@ -119,9 +120,8 @@ c Complete the calculation of density Gbar = 1.d0 dGbar = 0.d0 else - call get_shpfcn(shpi,lattce_meam(elti,elti)) - gam = (t_ave(1,i)*shpi(1)+t_ave(2,i)*shpi(2) - $ +t_ave(3,i)*shpi(3))/(Z*Z) + gam = (t_ave(1,i)*shp(1)+t_ave(2,i)*shp(2) + $ +t_ave(3,i)*shp(3))/(Z*Z) call dG_gam(gam,ibar_meam(elti),gsmooth_factor, $ Gbar,dGbar) endif From dd44189d1ff3947f67e006b80dcd39b9b629ff04 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 10 Jun 2017 04:35:11 -0400 Subject: [PATCH 18/22] fix bug in compute orientorder/atom argument parsing --- src/compute_orientorder_atom.cpp | 33 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 17 deletions(-) diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp index 43f13ecc32..90e2830e39 100644 --- a/src/compute_orientorder_atom.cpp +++ b/src/compute_orientorder_atom.cpp @@ -102,23 +102,22 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg if (qlist[iw] > qmax) qmax = qlist[iw]; } iarg += nqlist; - if (strcmp(arg[iarg],"components") == 0) { - qlcompflag = 1; - if (iarg+2 > narg) - error->all(FLERR,"Illegal compute orientorder/atom command"); - qlcomp = force->numeric(FLERR,arg[iarg+1]); - if (qlcomp <= 0) - error->all(FLERR,"Illegal compute orientorder/atom command"); - iqlcomp = -1; - for (int iw = 0; iw < nqlist; iw++) - if (qlcomp == qlist[iw]) { - iqlcomp = iw; - break; - } - if (iqlcomp < 0) - error->all(FLERR,"Illegal compute orientorder/atom command"); - iarg += 2; - } + } else if (strcmp(arg[iarg],"components") == 0) { + qlcompflag = 1; + if (iarg+2 > narg) + error->all(FLERR,"Illegal compute orientorder/atom command"); + qlcomp = force->numeric(FLERR,arg[iarg+1]); + if (qlcomp <= 0) + error->all(FLERR,"Illegal compute orientorder/atom command"); + iqlcomp = -1; + for (int iw = 0; iw < nqlist; iw++) + if (qlcomp == qlist[iw]) { + iqlcomp = iw; + break; + } + if (iqlcomp < 0) + error->all(FLERR,"Illegal compute orientorder/atom command"); + iarg += 2; } else if (strcmp(arg[iarg],"cutoff") == 0) { if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); From daa77176add6f9edb3058a54584c398b0515898b Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Sat, 10 Jun 2017 17:28:17 -0400 Subject: [PATCH 19/22] add restart support to fix deform. only "initial" data is restored and some consistency check performed --- doc/src/fix_deform.txt | 6 +++-- src/fix_deform.cpp | 55 +++++++++++++++++++++++++++++++++--------- src/fix_deform.h | 2 ++ 3 files changed, 49 insertions(+), 14 deletions(-) diff --git a/doc/src/fix_deform.txt b/doc/src/fix_deform.txt index 8c3a9fa499..d3254eece6 100644 --- a/doc/src/fix_deform.txt +++ b/doc/src/fix_deform.txt @@ -565,8 +565,10 @@ more instructions on how to use the accelerated styles effectively. [Restart, fix_modify, output, run start/stop, minimize info:] -No information about this fix is written to "binary restart -files"_restart.html. None of the "fix_modify"_fix_modify.html options +This fix will restore the initial box settings from "binary restart +files"_restart.html, which allows the fix to be properly continue +deformation, when using the start/stop options of the "run"_run.html +command. None of the "fix_modify"_fix_modify.html options are relevant to this fix. No global or per-atom quantities are stored by this fix for access by various "output commands"_Section_howto.html#howto_15. diff --git a/src/fix_deform.cpp b/src/fix_deform.cpp index 171a90ba3a..705f1970af 100644 --- a/src/fix_deform.cpp +++ b/src/fix_deform.cpp @@ -37,7 +37,7 @@ using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -enum{NONE,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME,WIGGLE,VARIABLE}; +enum{NONE=0,FINAL,DELTA,SCALE,VEL,ERATE,TRATE,VOLUME,WIGGLE,VARIABLE}; enum{ONE_FROM_ONE,ONE_FROM_TWO,TWO_FROM_ONE}; // same as domain.cpp, fix_nvt_sllod.cpp, compute_temp_deform.cpp @@ -52,6 +52,7 @@ rfix(NULL), irregular(NULL), set(NULL) if (narg < 4) error->all(FLERR,"Illegal fix deform command"); no_change_box = 1; + restart_global = 1; nevery = force->inumeric(FLERR,arg[3]); if (nevery <= 0) error->all(FLERR,"Illegal fix deform command"); @@ -59,12 +60,7 @@ rfix(NULL), irregular(NULL), set(NULL) // set defaults set = new Set[6]; - set[0].style = set[1].style = set[2].style = - set[3].style = set[4].style = set[5].style = NONE; - set[0].hstr = set[1].hstr = set[2].hstr = - set[3].hstr = set[4].hstr = set[5].hstr = NULL; - set[0].hratestr = set[1].hratestr = set[2].hratestr = - set[3].hratestr = set[4].hratestr = set[5].hratestr = NULL; + memset(set,0,6*sizeof(Set)); // parse arguments @@ -343,11 +339,9 @@ rfix(NULL), irregular(NULL), set(NULL) set[i].hi_initial = domain->boxhi[i]; set[i].vol_initial = domain->xprd * domain->yprd * domain->zprd; } - for (int i = 3; i < 6; i++) { - if (i == 5) set[i].tilt_initial = domain->xy; - else if (i == 4) set[i].tilt_initial = domain->xz; - else if (i == 3) set[i].tilt_initial = domain->yz; - } + set[3].tilt_initial = domain->yz; + set[4].tilt_initial = domain->xz; + set[5].tilt_initial = domain->xy; // reneighboring only forced if flips can occur due to shape changes @@ -955,6 +949,43 @@ void FixDeform::end_of_step() if (kspace_flag) force->kspace->setup(); } +/* ---------------------------------------------------------------------- + write Set data to restart file +------------------------------------------------------------------------- */ + +void FixDeform::write_restart(FILE *fp) +{ + if (comm->me == 0) { + int size = 6*sizeof(Set); + fwrite(&size,sizeof(int),1,fp); + fwrite(set,sizeof(Set),6,fp); + } +} + +/* ---------------------------------------------------------------------- + use selected state info from restart file to restart the Fix +------------------------------------------------------------------------- */ + +void FixDeform::restart(char *buf) +{ + int samestyle = 1; + Set *set_restart = (Set *) buf; + for (int i=0; i<6; ++i) { + // restore data from initial state + set[i].lo_initial = set_restart[i].lo_initial; + set[i].hi_initial = set_restart[i].hi_initial; + set[i].vol_initial = set_restart[i].vol_initial; + set[i].tilt_initial = set_restart[i].tilt_initial; + // check if style settings are consitent (should do the whole set?) + if (set[i].style != set_restart[i].style) + samestyle = 0; + if (set[i].substyle != set_restart[i].substyle) + samestyle = 0; + } + if (!samestyle) + error->all(FLERR,"Fix deform settings not consistent with restart"); +} + /* ---------------------------------------------------------------------- */ void FixDeform::options(int narg, char **arg) diff --git a/src/fix_deform.h b/src/fix_deform.h index cdda1b8547..4d440eb3c7 100644 --- a/src/fix_deform.h +++ b/src/fix_deform.h @@ -35,6 +35,8 @@ class FixDeform : public Fix { void init(); virtual void pre_exchange(); virtual void end_of_step(); + virtual void write_restart(FILE *); + virtual void restart(char *buf); double memory_usage(); protected: From 6b289b07942dbf9d0f86a8047e8205eccf7cbd48 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 12 Jun 2017 14:07:40 -0400 Subject: [PATCH 20/22] change incorrect EINERTIA constant in rigid body integrators from 4.0 to 2.0 (same as in other integrators) --- src/RIGID/fix_rigid.cpp | 2 +- src/RIGID/fix_rigid_small.cpp | 2 +- src/USER-OMP/fix_rigid_omp.cpp | 2 +- src/USER-OMP/fix_rigid_small_omp.cpp | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index e9c06dee5c..708b81d7f0 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -50,7 +50,7 @@ enum{ISO,ANISO,TRICLINIC}; #define EPSILON 1.0e-7 #define SINERTIA 0.4 // moment of inertia prefactor for sphere -#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid +#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid #define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment /* ---------------------------------------------------------------------- */ diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 2d8e736a1e..07b6cc09b3 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -54,7 +54,7 @@ FixRigidSmall *FixRigidSmall::frsptr; #define BIG 1.0e20 #define SINERTIA 0.4 // moment of inertia prefactor for sphere -#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid +#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid #define LINERTIA (1.0/12.0) // moment of inertia prefactor for line segment #define DELTA_BODY 10000 diff --git a/src/USER-OMP/fix_rigid_omp.cpp b/src/USER-OMP/fix_rigid_omp.cpp index d5bacc0320..f3ef1d2985 100644 --- a/src/USER-OMP/fix_rigid_omp.cpp +++ b/src/USER-OMP/fix_rigid_omp.cpp @@ -39,7 +39,7 @@ using namespace MathConst; enum{SINGLE,MOLECULE,GROUP}; // same as in FixRigid -#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid +#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid typedef struct { double x,y,z; } dbl3_t; diff --git a/src/USER-OMP/fix_rigid_small_omp.cpp b/src/USER-OMP/fix_rigid_small_omp.cpp index a260899aef..e3939a829d 100644 --- a/src/USER-OMP/fix_rigid_small_omp.cpp +++ b/src/USER-OMP/fix_rigid_small_omp.cpp @@ -37,7 +37,7 @@ using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; -#define EINERTIA 0.4 // moment of inertia prefactor for ellipsoid +#define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid enum{FULL_BODY,INITIAL,FINAL,FORCE_TORQUE,VCM_ANGMOM,XCM_MASS,ITENSOR,DOF}; From 0c805d0b70ed81b24cfcc340ccef18649d3ffdc6 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 12 Jun 2017 14:20:38 -0400 Subject: [PATCH 21/22] correctly skip over point particles and point dipoles when counting extendend particles in fix rigid/small --- src/RIGID/fix_rigid_small.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 07b6cc09b3..25e528e40b 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -1031,7 +1031,7 @@ int FixRigidSmall::dof(int tgroup) j = atom2body[i]; counts[j][2]++; if (mask[i] & tgroupbit) { - if (extended && eflags[i]) counts[j][1]++; + if (extended && (eflags[i] & ~(POINT | DIPOLE))) counts[j][1]++; else counts[j][0]++; } } From b62d526cc98d53c5c81d9839977dd8fc2dd1ffda Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Thu, 15 Jun 2017 11:01:36 -0400 Subject: [PATCH 22/22] Revert "avoid undesired negative forces for high particle velocities in granular models" This reverts commit 066123007cd9277d4c03669707cf34ab6c05bc06. --- src/GRANULAR/pair_gran_hertz_history.cpp | 2 -- src/GRANULAR/pair_gran_hooke.cpp | 2 -- src/GRANULAR/pair_gran_hooke_history.cpp | 2 -- 3 files changed, 6 deletions(-) diff --git a/src/GRANULAR/pair_gran_hertz_history.cpp b/src/GRANULAR/pair_gran_hertz_history.cpp index 14cbc8eb0a..e52aac10db 100644 --- a/src/GRANULAR/pair_gran_hertz_history.cpp +++ b/src/GRANULAR/pair_gran_hertz_history.cpp @@ -183,7 +183,6 @@ void PairGranHertzHistory::compute(int eflag, int vflag) ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; - if (ccel < 0.0) ccel = 0.0; // relative velocities @@ -391,7 +390,6 @@ double PairGranHertzHistory::single(int i, int j, int itype, int jtype, ccel = kn*(radsum-r)*rinv - damp; polyhertz = sqrt((radsum-r)*radi*radj / radsum); ccel *= polyhertz; - if (ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke.cpp b/src/GRANULAR/pair_gran_hooke.cpp index 2d6d113ab7..66e1c3dd7a 100644 --- a/src/GRANULAR/pair_gran_hooke.cpp +++ b/src/GRANULAR/pair_gran_hooke.cpp @@ -161,7 +161,6 @@ void PairGranHooke::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if (ccel < 0.0) ccel = 0.0; // relative velocities @@ -303,7 +302,6 @@ double PairGranHooke::single(int i, int j, int itype, int jtype, double rsq, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if (ccel < 0.0) ccel = 0.0; // relative velocities diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 1378709cb5..e9662c9e73 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -223,7 +223,6 @@ void PairGranHookeHistory::compute(int eflag, int vflag) damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if (ccel < 0.0) ccel = 0.0; // relative velocities @@ -688,7 +687,6 @@ double PairGranHookeHistory::single(int i, int j, int itype, int jtype, damp = meff*gamman*vnnr*rsqinv; ccel = kn*(radsum-r)*rinv - damp; - if (ccel < 0.0) ccel = 0.0; // relative velocities