diff --git a/doc/src/Section_commands.txt b/doc/src/Section_commands.txt index 010e45a34e..7284b4b0f6 100644 --- a/doc/src/Section_commands.txt +++ b/doc/src/Section_commands.txt @@ -702,6 +702,8 @@ package"_Section_start.html#start_3. "meso"_fix_meso.html, "manifoldforce"_fix_manifoldforce.html, "meso/stationary"_fix_meso_stationary.html, +"nve/dot"_fix_nve_dot.html, +"nve/dotc/langevin"_fix_nve_dotc_langevin.html, "nve/manifold/rattle"_fix_nve_manifold_rattle.html, "nvk"_fix_nvk.html, "nvt/manifold/rattle"_fix_nvt_manifold_rattle.html, @@ -1035,6 +1037,11 @@ package"_Section_start.html#start_3. "morse/soft"_pair_morse.html, "multi/lucy"_pair_multi_lucy.html, "multi/lucy/rx"_pair_multi_lucy_rx.html, +"oxdna/coaxstk"_pair_oxdna.html, +"oxdna/excv"_pair_oxdna.html, +"oxdna/hbond"_pair_oxdna.html, +"oxdna/stk"_pair_oxdna.html, +"oxdna/xstk"_pair_oxdna.html, "quip"_pair_quip.html, "reax/c (k)"_pair_reax_c.html, "smd/hertz"_pair_smd_hertz.html, @@ -1083,7 +1090,8 @@ if "LAMMPS is built with the appropriate package"_Section_start.html#start_3. "harmonic/shift (o)"_bond_harmonic_shift.html, -"harmonic/shift/cut (o)"_bond_harmonic_shift_cut.html :tb(c=4,ea=c) +"harmonic/shift/cut (o)"_bond_harmonic_shift_cut.html, +"oxdna/fene"_bond_oxdna_fene.html :tb(c=4,ea=c) :line diff --git a/doc/src/Section_packages.txt b/doc/src/Section_packages.txt index 61e667e933..b9b96d66ad 100644 --- a/doc/src/Section_packages.txt +++ b/doc/src/Section_packages.txt @@ -84,6 +84,7 @@ Package, Description, Author(s), Doc page, Example, Library "PERI"_#PERI, Peridynamics models, Mike Parks (Sandia), "pair_style peri"_pair_peri.html, peri, - "POEMS"_#POEMS, coupled rigid body motion, Rudra Mukherjee (JPL), "fix poems"_fix_poems.html, rigid, lib/poems "PYTHON"_#PYTHON, embed Python code in an input script, -, "python"_python.html, python, lib/python +"REAX"_#REAX, ReaxFF potential, Aidan Thompson (Sandia), "pair_style reax"_pair_reax.html, reax, lib/reax "REPLICA"_#REPLICA, multi-replica methods, -, "Section 6.6.5"_Section_howto.html#howto_5, tad, - "RIGID"_#RIGID, rigid bodies, -, "fix rigid"_fix_rigid.html, rigid, - "SHOCK"_#SHOCK, shock loading methods, -, "fix msst"_fix_msst.html, -, - @@ -1286,27 +1287,26 @@ him directly if you have questions. USER-CGDNA package :link(USER-CGDNA),h5 -Contents: The CGDNA package implements coarse-grained force fields -for single- and double-stranded DNA. This is at the moment mainly -the oxDNA model, developed by Doye, Louis and Ouldridge at the -University of Oxford. -The package also contains Langevin-type rigid-body integrators -with improved stability. +Contents: The CGDNA package implements coarse-grained force fields for +single- and double-stranded DNA. This is at the moment mainly the +oxDNA model, developed by Doye, Louis and Ouldridge at the University +of Oxford. The package also contains Langevin-type rigid-body +integrators with improved stability. See these doc pages to get started: "bond_style oxdna_fene"_bond_oxdna_fene.html - "pair_style oxdna_excv"_pair_oxdna_excv.html +"fix nve/dotc/langevin"_fix_nve_dotc_langevin.html :ul -"fix nve/dotc/langevin"_fix_nve_dotc_langevin.html - -Supporting info: /src/USER-CGDNA/README, "bond_style oxdna_fene"_bond_oxdna_fene.html, -"pair_style oxdna_excv"_pair_oxdna_excv.html, "fix nve/dotc/langevin"_fix_nve_dotc_langevin.html +Supporting info: /src/USER-CGDNA/README, "bond_style +oxdna_fene"_bond_oxdna_fene.html, "pair_style +oxdna_excv"_pair_oxdna_excv.html, "fix +nve/dotc/langevin"_fix_nve_dotc_langevin.html Author: Oliver Henrich at the University of Edinburgh, UK (o.henrich -at epcc.ed.ac.uk or ohenrich at ph.ed.ac.uk). Contact him directly -if you have any questions. +at epcc.ed.ac.uk or ohenrich at ph.ed.ac.uk). Contact him directly if +you have any questions. :line diff --git a/doc/src/compute_coord_atom.txt b/doc/src/compute_coord_atom.txt index f1a6bf7ffb..0367b2fe95 100644 --- a/doc/src/compute_coord_atom.txt +++ b/doc/src/compute_coord_atom.txt @@ -12,19 +12,16 @@ compute coord/atom command :h3 compute ID group-ID coord/atom cstyle args ... :pre -ID, group-ID are documented in "compute"_compute.html command -coord/atom = style name of this compute command -one cstyle must be appended :ul - -cstyle = {cutoff} or {orientorder} - -{cutoff} args = cutoff typeN - cutoff = distance within which to count coordination neighbors (distance units) - typeN = atom type for Nth coordination count (see asterisk form below) :pre - -{orientorder} args = orientorderID threshold - orientorderID = ID of a previously defined orientorder/atom compute - threshold = minimum value of the scalar product between two 'connected' atoms (see text for explanation) :pre +ID, group-ID are documented in "compute"_compute.html command :ulb,l +coord/atom = style name of this compute command :l +cstyle = {cutoff} or {orientorder} :l + {cutoff} args = cutoff typeN + cutoff = distance within which to count coordination neighbors (distance units) + typeN = atom type for Nth coordination count (see asterisk form below) + {orientorder} args = orientorderID threshold + orientorderID = ID of an orientorder/atom compute + threshold = minimum value of the product of two "connected" atoms :pre +:ule [Examples:] @@ -35,21 +32,21 @@ compute 1 all coord/atom orientorder 2 0.5 :pre [Description:] -This compute performs generic calculations between neighboring atoms. So far, -there are two cstyles implemented: {cutoff} and {orientorder}. -The {cutoff} cstyle calculates one or more coordination numbers -for each atom in a group. +This compute performs calculations between neighboring atoms to +determine a coordination value. The specific calculation and the +meaning of the resulting value depend on the {cstyle} keyword used. -A coordination number is defined as the number of neighbor atoms with -specified atom type(s) that are within the specified cutoff distance -from the central atom. Atoms not in the group are included in a -coordination number of atoms in the group. +The {cutoff} cstyle calculates one or more traditional coordination +numbers for each atom. A coordination number is defined as the number +of neighbor atoms with specified atom type(s) that are within the +specified cutoff distance from the central atom. Atoms not in the +specified group are included in the coordination number tally. -The {typeN} keywords allow you to specify which atom types contribute -to each coordination number. One coordination number is computed for -each of the {typeN} keywords listed. If no {typeN} keywords are -listed, a single coordination number is calculated, which includes -atoms of all types (same as the "*" format, see below). +The {typeN} keywords allow specification of which atom types +contribute to each coordination number. One coordination number is +computed for each of the {typeN} keywords listed. If no {typeN} +keywords are listed, a single coordination number is calculated, which +includes atoms of all types (same as the "*" format, see below). The {typeN} keywords can be specified in one of two ways. An explicit numeric value can be used, as in the 2nd example above. Or a @@ -61,16 +58,27 @@ from 1 to N. A leading asterisk means all types from 1 to n (inclusive). A middle asterisk means all types from m to n (inclusive). -The {orientorder} cstyle calculates the number of 'connected' atoms j -around each atom i. The atom j is connected to i if the scalar product -({Ybar_lm(i)},{Ybar_lm(j)}) is larger than {threshold}. Thus, this cstyle -will work only if a "compute orientorder/atom"_compute_orientorder_atom.html -has been previously defined. This cstyle allows one to apply the -ten Wolde's criterion to identify cristal-like atoms in a system -(see "ten Wolde et al."_#tenWolde). +The {orientorder} cstyle calculates the number of "connected" neighbor +atoms J around each central atom I. For this {cstyle}, connected is +defined by the orientational order parameter calculated by the +"compute orientorder/atom"_compute_orientorder_atom.html command. +This {cstyle} thus allows one to apply the ten Wolde's criterion to +identify crystal-like atoms in a system, as discussed in "ten +Wolde"_#tenWolde. -The value of all coordination numbers will be 0.0 for atoms not in the -specified compute group. +The ID of the previously specified "compute +orientorder/atom"_compute_orientorder/atom command is specified as +{orientorderID}. The compute must invoke its {components} option to +calculate components of the {Ybar_lm} vector for each atoms, as +described in its documenation. Note that orientorder/atom compute +defines its own criteria for identifying neighboring atoms. If the +scalar product ({Ybar_lm(i)},{Ybar_lm(j)}), calculated by the +orientorder/atom compute is larger than the specified {threshold}, +then I and J are connected, and the coordination value of I is +incremented by one. + +For all {cstyle} settings, all coordination values will be 0.0 for +atoms not in the specified compute group. The neighbor list needed to compute this quantity is constructed each time the calculation is performed (i.e. each time a snapshot of atoms @@ -92,21 +100,23 @@ the neighbor list. [Output info:] -If single {type1} keyword is specified (or if none are specified), -this compute calculates a per-atom vector. If multiple {typeN} -keywords are specified, this compute calculates a per-atom array, with -N columns. These values can be accessed by any command that uses -per-atom values from a compute as input. See "Section +For {cstyle} cutoff, this compute can calculate a per-atom vector or +array. If single {type1} keyword is specified (or if none are +specified), this compute calculates a per-atom vector. If multiple +{typeN} keywords are specified, this compute calculates a per-atom +array, with N columns. + +For {cstyle} orientorder, this compute calculates a per-atom vector. + +These values can be accessed by any command that uses per-atom values +from a compute as input. See "Section 6.15"_Section_howto.html#howto_15 for an overview of LAMMPS output options. The per-atom vector or array values will be a number >= 0.0, as explained above. -[Restrictions:] -The cstyle {orientorder} can only be used if a -"compute orientorder/atom"_compute_orientorder_atom.html command -was previously defined. Otherwise, an error message will be issued. +[Restrictions:] none [Related commands:] @@ -118,4 +128,5 @@ was previously defined. Otherwise, an error message will be issued. :line :link(tenWolde) -[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, J. Chem. Phys. 104, 9932 (1996). +[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, +J. Chem. Phys. 104, 9932 (1996). diff --git a/doc/src/compute_orientorder_atom.txt b/doc/src/compute_orientorder_atom.txt index 74426dd33b..2aa76ab5ac 100644 --- a/doc/src/compute_orientorder_atom.txt +++ b/doc/src/compute_orientorder_atom.txt @@ -19,7 +19,7 @@ keyword = {cutoff} or {nnn} or {degrees} or {components} {cutoff} value = distance cutoff {nnn} value = number of nearest neighbors {degrees} values = nlvalues, l1, l2,... - {components} value = l :pre + {components} value = ldegree :pre :ule @@ -64,21 +64,21 @@ specified distance cutoff are used. The optional keyword {degrees} defines the list of order parameters to be computed. The first argument {nlvalues} is the number of order parameters. This is followed by that number of integers giving the -degree of each order parameter. Because {Q}2 and all odd-degree -order parameters are zero for atoms in cubic crystals -(see "Steinhardt"_#Steinhardt), the default order parameters -are {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12. For the -FCC crystal with {nnn}=12, {Q}4 = sqrt(7/3)/8 = 0.19094.... -The numerical values of all order parameters up to {Q}12 -for a range of commonly encountered high-symmetry structures are given -in Table I of "Mickel et al."_#Mickel. +degree of each order parameter. Because {Q}2 and all odd-degree order +parameters are zero for atoms in cubic crystals (see +"Steinhardt"_#Steinhardt), the default order parameters are {Q}4, +{Q}6, {Q}8, {Q}10, and {Q}12. For the FCC crystal with {nnn}=12, {Q}4 += sqrt(7/3)/8 = 0.19094.... The numerical values of all order +parameters up to {Q}12 for a range of commonly encountered +high-symmetry structures are given in Table I of "Mickel et +al."_#Mickel. -The optional keyword {components} will output the components of -the normalized complex vector {Ybar_lm} of degree {l}, which must be +The optional keyword {components} will output the components of the +normalized complex vector {Ybar_lm} of degree {ldegree}, which must be explicitly included in the keyword {degrees}. This option can be used in conjunction with "compute coord_atom"_compute_coord_atom.html to -calculate the ten Wolde's criterion to identify crystal-like particles -(see "ten Wolde et al."_#tenWolde96). +calculate the ten Wolde's criterion to identify crystal-like +particles, as discussed in "ten Wolde"_#tenWolde. The value of {Ql} is set to zero for atoms not in the specified compute group, as well as for atoms that have less than @@ -104,14 +104,16 @@ the neighbor list. [Output info:] -This compute calculates a per-atom array with {nlvalues} columns, giving the -{Ql} values for each atom, which are real numbers on the range 0 <= {Ql} <= 1. +This compute calculates a per-atom array with {nlvalues} columns, +giving the {Ql} values for each atom, which are real numbers on the +range 0 <= {Ql} <= 1. -If the keyword {components} is set, then the real and imaginary parts of each -component of (normalized) {Ybar_lm} will be added to the output array in the -following order: -Re({Ybar_-m}) Im({Ybar_-m}) Re({Ybar_-m+1}) Im({Ybar_-m+1}) ... Re({Ybar_m}) Im({Ybar_m}). -This way, the per-atom array will have a total of {nlvalues}+2*(2{l}+1) columns. +If the keyword {components} is set, then the real and imaginary parts +of each component of (normalized) {Ybar_lm} will be added to the +output array in the following order: Re({Ybar_-m}) Im({Ybar_-m}) +Re({Ybar_-m+1}) Im({Ybar_-m+1}) ... Re({Ybar_m}) Im({Ybar_m}). This +way, the per-atom array will have a total of {nlvalues}+2*(2{l}+1) +columns. These values can be accessed by any command that uses per-atom values from a compute as input. See "Section @@ -122,19 +124,25 @@ options. [Related commands:] -"compute coord/atom"_compute_coord_atom.html, "compute centro/atom"_compute_centro_atom.html, "compute hexorder/atom"_compute_hexorder_atom.html +"compute coord/atom"_compute_coord_atom.html, "compute +centro/atom"_compute_centro_atom.html, "compute +hexorder/atom"_compute_hexorder_atom.html [Default:] -The option defaults are {cutoff} = pair style cutoff, {nnn} = 12, {degrees} = 5 4 6 8 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12. +The option defaults are {cutoff} = pair style cutoff, {nnn} = 12, +{degrees} = 5 4 6 8 10 12 i.e. {Q}4, {Q}6, {Q}8, {Q}10, and {Q}12. :line :link(Steinhardt) -[(Steinhardt)] P. Steinhardt, D. Nelson, and M. Ronchetti, Phys. Rev. B 28, 784 (1983). +[(Steinhardt)] P. Steinhardt, D. Nelson, and M. Ronchetti, +Phys. Rev. B 28, 784 (1983). :link(Mickel) -[(Mickel)] W. Mickel, S. C. Kapfer, G. E. Schroeder-Turkand, K. Mecke, J. Chem. Phys. 138, 044501 (2013). +[(Mickel)] W. Mickel, S. C. Kapfer, G. E. Schroeder-Turkand, K. Mecke, +J. Chem. Phys. 138, 044501 (2013). -:link(tenWolde96) -[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, J. Chem. Phys. 104, 9932 (1996). +:link(tenWolde) +[(tenWolde)] P. R. ten Wolde, M. J. Ruiz-Montero, D. Frenkel, +J. Chem. Phys. 104, 9932 (1996). diff --git a/src/KOKKOS/nbin_kokkos.cpp b/src/KOKKOS/nbin_kokkos.cpp index 9a73e23b42..5e41787247 100644 --- a/src/KOKKOS/nbin_kokkos.cpp +++ b/src/KOKKOS/nbin_kokkos.cpp @@ -74,10 +74,7 @@ void NBinKokkos::bin_atoms_setup(int nall) k_bincount = DAT::tdual_int_1d("Neighbor::d_bincount",mbins); bincount = k_bincount.view(); - last_bin_memory = update->ntimestep; } - - last_bin = update->ntimestep; } /* ---------------------------------------------------------------------- @@ -87,6 +84,8 @@ void NBinKokkos::bin_atoms_setup(int nall) template void NBinKokkos::bin_atoms() { + last_bin = update->ntimestep; + h_resize() = 1; while(h_resize() > 0) { @@ -116,7 +115,6 @@ void NBinKokkos::bin_atoms() k_bins = DAT::tdual_int_2d("bins", mbins, atoms_per_bin); bins = k_bins.view(); c_bins = bins; - last_bin_memory = update->ntimestep; } } } diff --git a/src/MC/fix_atom_swap.cpp b/src/MC/fix_atom_swap.cpp index 3156e851ea..80bebefec0 100644 --- a/src/MC/fix_atom_swap.cpp +++ b/src/MC/fix_atom_swap.cpp @@ -57,7 +57,10 @@ using namespace MathConst; FixAtomSwap::FixAtomSwap(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - idregion(NULL), type_list(NULL), mu(NULL), qtype(NULL), sqrt_mass_ratio(NULL), local_swap_iatom_list(NULL), local_swap_jatom_list(NULL), local_swap_atom_list(NULL), random_equal(NULL), random_unequal(NULL), c_pe(NULL) + idregion(NULL), type_list(NULL), mu(NULL), qtype(NULL), + sqrt_mass_ratio(NULL), local_swap_iatom_list(NULL), + local_swap_jatom_list(NULL), local_swap_atom_list(NULL), + random_equal(NULL), random_unequal(NULL), c_pe(NULL) { if (narg < 10) error->all(FLERR,"Illegal fix atom/swap command"); diff --git a/src/USER-CGDNA/README b/src/USER-CGDNA/README index 8415e12fc9..8158b83ab8 100644 --- a/src/USER-CGDNA/README +++ b/src/USER-CGDNA/README @@ -17,21 +17,20 @@ self-assembly, DPhil. University of Oxford (2011). gradient thermostats for rigid body dynamics", J. Chem. Phys. 142, 144114 (2015). -Example input and data files can be found in -/examples/USER/cgdna/examples/duplex1/ and /duplex2/. -A simple python setup tool which creates single straight or helical DNA -strands as well as DNA duplexes and arrays of duplexes can be found in -/examples/USER/cgdna/util/. -A technical report with more information on the model, the structure -of the input and data file, the setup tool and the performance of the -LAMMPS-implementation of oxDNA can be found in -/doc/src/PDF/USER-CGDNA-overview.pdf. +Example input and data files can be found in +/examples/USER/cgdna/examples/duplex1/ and /duplex2/. A simple python +setup tool which creates single straight or helical DNA strands as +well as DNA duplexes and arrays of duplexes can be found in +/examples/USER/cgdna/util/. A technical report with more information +on the model, the structure of the input and data file, the setup tool +and the performance of the LAMMPS-implementation of oxDNA can be found +in /doc/src/PDF/USER-CGDNA-overview.pdf. IMPORTANT NOTE: This package can only be used if LAMMPS is compiled -with the MOLECULE and ASPHERE packages. These should be included -in the LAMMPS build by typing "make yes-asphere yes-molecule" prior -to the usual compilation (see the "Including/excluding packages" -section of the LAMMPS manual). +with the MOLECULE and ASPHERE packages. These should be included in +the LAMMPS build by typing "make yes-asphere yes-molecule" prior to +the usual compilation (see the "Including/excluding packages" section +of the LAMMPS manual). The creator of this package is: @@ -39,14 +38,14 @@ Dr Oliver Henrich University of Edinburgh, UK ohenrich@ph.ed.ac.uk o.henrich@epcc.ed.ac.uk + -------------------------------------------------------------------------- -Bond styles provided by this package: +** Bond styles provided by this package: bond_oxdna_fene.cpp: backbone connectivity, a modified FENE potential - -Pair styles provided by this package: +** Pair styles provided by this package: pair_oxdna_excv.cpp: excluded volume interaction between the nucleotides @@ -60,8 +59,7 @@ pair_oxdna_xstk.cpp: cross-stacking interaction between nucleotides pair_oxdna_coaxstk.cpp: coaxial stacking interaction between nucleotides - -Fixes provided by this package: +** Fixes provided by this package: fix_nve_dotc_langevin.cpp: fix for Langevin-type rigid body integrator "C" in above Ref. [3] diff --git a/src/USER-DPD/nbin_ssa.cpp b/src/USER-DPD/nbin_ssa.cpp index d55acb6040..73da5e0df3 100644 --- a/src/USER-DPD/nbin_ssa.cpp +++ b/src/USER-DPD/nbin_ssa.cpp @@ -18,6 +18,7 @@ #include "nbin_ssa.h" #include "atom.h" +#include "update.h" #include "group.h" #include "memory.h" #include "error.h" @@ -59,6 +60,8 @@ void NBinSSA::bin_atoms() int *mask = atom->mask; int *ssaAIR = atom->ssaAIR; + last_bin = update->ntimestep; + for (i = 0; i < mbins; i++) { gbinhead_ssa[i] = -1; binhead_ssa[i] = -1; @@ -126,4 +129,3 @@ bigint NBinSSA::memory_usage() } return bytes; } - diff --git a/src/USER-INTEL/nbin_intel.cpp b/src/USER-INTEL/nbin_intel.cpp index 8dceafd3ab..194b9a5f97 100644 --- a/src/USER-INTEL/nbin_intel.cpp +++ b/src/USER-INTEL/nbin_intel.cpp @@ -86,7 +86,6 @@ void NBinIntel::bin_atoms_setup(int nall) nocopy(binhead:length(maxbin+1) alloc_if(1) free_if(0)) } #endif - last_bin_memory = update->ntimestep; } // bins = per-atom vector @@ -127,11 +126,7 @@ void NBinIntel::bin_atoms_setup(int nall) _fix->get_single_buffers()->set_bininfo(_atombin,_binpacked); else _fix->get_double_buffers()->set_bininfo(_atombin,_binpacked); - - last_bin_memory = update->ntimestep; } - - last_bin = update->ntimestep; } /* ---------------------------------------------------------------------- @@ -140,6 +135,8 @@ void NBinIntel::bin_atoms_setup(int nall) void NBinIntel::bin_atoms() { + last_bin = update->ntimestep; + if (_precision_mode == FixIntel::PREC_MODE_MIXED) bin_atoms(_fix->get_mixed_buffers()); else if (_precision_mode == FixIntel::PREC_MODE_SINGLE) diff --git a/src/USER-INTEL/npair_intel.cpp b/src/USER-INTEL/npair_intel.cpp index 8ec40260f3..bffb31b710 100644 --- a/src/USER-INTEL/npair_intel.cpp +++ b/src/USER-INTEL/npair_intel.cpp @@ -46,9 +46,7 @@ NPairIntel::~NPairIntel() { #endif } -/* ---------------------------------------------------------------------- - copy needed info from NStencil class to this build class -------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- */ #ifdef _LMP_INTEL_OFFLOAD void NPairIntel::grow_stencil() diff --git a/src/compute_coord_atom.cpp b/src/compute_coord_atom.cpp index 36f0b63504..1df5ad64ca 100644 --- a/src/compute_coord_atom.cpp +++ b/src/compute_coord_atom.cpp @@ -81,11 +81,11 @@ ComputeCoordAtom::ComputeCoordAtom(LAMMPS *lmp, int narg, char **arg) : if (iorientorder < 0) error->all(FLERR,"Could not find compute coord/atom compute ID"); if (strcmp(modify->compute[iorientorder]->style,"orientorder/atom") != 0) - error->all(FLERR,"Compute coord/atom compute ID does not compute orientorder/atom"); + error->all(FLERR,"Compute coord/atom compute ID is not orientorder/atom"); threshold = force->numeric(FLERR,arg[5]); if (threshold <= -1.0 || threshold >= 1.0) - error->all(FLERR,"Compute coord/atom threshold value must lie between -1 and 1"); + error->all(FLERR,"Compute coord/atom threshold not between -1 and 1"); ncol = 1; typelo = new int[ncol]; @@ -126,7 +126,7 @@ void ComputeCoordAtom::init() comm_forward = 2*(2*l+1); if (c_orientorder->iqlcomp < 0) error->all(FLERR,"Compute coord/atom requires components " - "option in compute orientorder/atom be defined"); + "option in compute orientorder/atom"); } if (force->pair == NULL) @@ -169,9 +169,6 @@ void ComputeCoordAtom::compute_peratom() invoked_peratom = update->ntimestep; -// printf("Number of degrees %i components degree %i",nqlist,l); -// printf("Particle \t %i \t Norm \t %g \n",0,norm[0][0]); - // grow coordination array if necessary if (atom->nmax > nmax) { @@ -195,8 +192,6 @@ void ComputeCoordAtom::compute_peratom() } nqlist = c_orientorder->nqlist; int ltmp = l; -// l = c_orientorder->qlcomp; - if (ltmp != l) error->all(FLERR,"Debug error, ltmp != l\n"); normv = c_orientorder->array_atom; comm->forward_comm_compute(this); } @@ -317,7 +312,7 @@ void ComputeCoordAtom::compute_peratom() /* ---------------------------------------------------------------------- */ int ComputeCoordAtom::pack_forward_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) + int pbc_flag, int *pbc) { int i,m=0,j; for (i = 0; i < n; ++i) { @@ -340,7 +335,6 @@ void ComputeCoordAtom::unpack_forward_comm(int n, int first, double *buf) normv[i][j] = buf[m++]; } } - } /* ---------------------------------------------------------------------- diff --git a/src/compute_orientorder_atom.cpp b/src/compute_orientorder_atom.cpp index 5f78b33b61..43f13ecc32 100644 --- a/src/compute_orientorder_atom.cpp +++ b/src/compute_orientorder_atom.cpp @@ -73,7 +73,8 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg int iarg = 3; while (iarg < narg) { if (strcmp(arg[iarg],"nnn") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); + if (iarg+2 > narg) + error->all(FLERR,"Illegal compute orientorder/atom command"); if (strcmp(arg[iarg+1],"NULL") == 0) { nnn = 0; } else { @@ -91,7 +92,8 @@ ComputeOrientOrderAtom::ComputeOrientOrderAtom(LAMMPS *lmp, int narg, char **arg memory->destroy(qlist); memory->create(qlist,nqlist,"orientorder/atom:qlist"); iarg += 2; - if (iarg+nqlist > narg) error->all(FLERR,"Illegal compute orientorder/atom command"); + if (iarg+nqlist > narg) + error->all(FLERR,"Illegal compute orientorder/atom command"); qmax = 0; for (int iw = 0; iw < nqlist; iw++) { qlist[iw] = force->numeric(FLERR,arg[iarg+iw]); @@ -157,11 +159,12 @@ ComputeOrientOrderAtom::~ComputeOrientOrderAtom() void ComputeOrientOrderAtom::init() { if (force->pair == NULL) - error->all(FLERR,"Compute orientorder/atom requires a pair style be defined"); + error->all(FLERR,"Compute orientorder/atom requires a " + "pair style be defined"); if (cutsq == 0.0) cutsq = force->pair->cutforce * force->pair->cutforce; else if (sqrt(cutsq) > force->pair->cutforce) - error->all(FLERR, - "Compute orientorder/atom cutoff is longer than pairwise cutoff"); + error->all(FLERR,"Compute orientorder/atom cutoff is " + "longer than pairwise cutoff"); memory->create(qnm_r,qmax,2*qmax+1,"orientorder/atom:qnm_r"); memory->create(qnm_i,qmax,2*qmax+1,"orientorder/atom:qnm_i"); @@ -488,7 +491,8 @@ void ComputeOrientOrderAtom::calc_boop(double **rlist, calculate scalar distance ------------------------------------------------------------------------- */ -double ComputeOrientOrderAtom::dist(const double r[]) { +double ComputeOrientOrderAtom::dist(const double r[]) +{ return sqrt(r[0]*r[0] + r[1]*r[1] + r[2]*r[2]); } @@ -497,8 +501,8 @@ double ComputeOrientOrderAtom::dist(const double r[]) { Y_l^m (theta, phi) = prefactor(l, m, cos(theta)) * exp(i*m*phi) ------------------------------------------------------------------------- */ -double ComputeOrientOrderAtom:: -polar_prefactor(int l, int m, double costheta) { +double ComputeOrientOrderAtom::polar_prefactor(int l, int m, double costheta) +{ const int mabs = abs(m); double prefactor = 1.0; @@ -517,8 +521,8 @@ polar_prefactor(int l, int m, double costheta) { associated legendre polynomial ------------------------------------------------------------------------- */ -double ComputeOrientOrderAtom:: -associated_legendre(int l, int m, double x) { +double ComputeOrientOrderAtom::associated_legendre(int l, int m, double x) +{ if (l < m) return 0.0; double p(1.0), pm1(0.0), pm2(0.0); diff --git a/src/compute_rdf.cpp b/src/compute_rdf.cpp index 7b8b808d17..27fd09f5e0 100644 --- a/src/compute_rdf.cpp +++ b/src/compute_rdf.cpp @@ -40,7 +40,8 @@ using namespace MathConst; ComputeRDF::ComputeRDF(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), rdfpair(NULL), nrdfpair(NULL), ilo(NULL), ihi(NULL), jlo(NULL), jhi(NULL), - hist(NULL), histall(NULL), typecount(NULL), icount(NULL), jcount(NULL), duplicates(NULL) + hist(NULL), histall(NULL), typecount(NULL), icount(NULL), jcount(NULL), + duplicates(NULL) { if (narg < 4 || (narg-4) % 2) error->all(FLERR,"Illegal compute rdf command"); diff --git a/src/nbin.cpp b/src/nbin.cpp index ef8543fdf9..aa7f6954d4 100644 --- a/src/nbin.cpp +++ b/src/nbin.cpp @@ -24,7 +24,7 @@ using namespace LAMMPS_NS; NBin::NBin(LAMMPS *lmp) : Pointers(lmp) { - last_setup = last_bin = last_bin_memory = -1; + last_bin = -1; maxbin = maxatom = 0; binhead = NULL; bins = NULL; @@ -71,7 +71,6 @@ void NBin::bin_atoms_setup(int nall) maxbin = mbins; memory->destroy(binhead); memory->create(binhead,maxbin,"neigh:binhead"); - last_bin_memory = update->ntimestep; } // bins = per-atom vector @@ -80,10 +79,7 @@ void NBin::bin_atoms_setup(int nall) maxatom = nall; memory->destroy(bins); memory->create(bins,maxatom,"neigh:bins"); - last_bin_memory = update->ntimestep; } - - last_bin = update->ntimestep; } /* ---------------------------------------------------------------------- diff --git a/src/nbin.h b/src/nbin.h index f000df75b0..0599bd0a99 100644 --- a/src/nbin.h +++ b/src/nbin.h @@ -21,9 +21,7 @@ namespace LAMMPS_NS { class NBin : protected Pointers { public: int istyle; // 1-N index into binnames - - bigint last_setup,last_bin; // timesteps for last operations performed - bigint last_bin_memory; + bigint last_bin; // last timestep atoms were binned int nbinx,nbiny,nbinz; // # of global bins int mbins; // # of local bins and offset on this proc diff --git a/src/nbin_standard.cpp b/src/nbin_standard.cpp index 97257df717..2a72d996a5 100644 --- a/src/nbin_standard.cpp +++ b/src/nbin_standard.cpp @@ -54,8 +54,6 @@ NBinStandard::NBinStandard(LAMMPS *lmp) : NBin(lmp) {} void NBinStandard::setup_bins(int style) { - last_setup = update->ntimestep; - // bbox = size of bbox of entire domain // bsubbox lo/hi = bounding box of my subdomain extended by comm->cutghost // for triclinic: @@ -197,6 +195,7 @@ void NBinStandard::bin_atoms() { int i,ibin; + last_bin = update->ntimestep; for (i = 0; i < mbins; i++) binhead[i] = -1; // bin in reverse order so linked list will be in forward order diff --git a/src/neighbor.cpp b/src/neighbor.cpp index 834bdffb7a..9d6b93819a 100644 --- a/src/neighbor.cpp +++ b/src/neighbor.cpp @@ -695,27 +695,38 @@ void Neighbor::init_pair() } // (C) rule: - // for fix/compute requests, occasional or not does not matter + // for fix/compute requests to be copied: // 1st check: + // occasional or not does not matter + // Kokkos host/device flags must match + // SSA flag must match + // newton setting must match (else list has different neighbors in it) + // 2nd check: // if request = half and non-skip/copy pair half/respaouter request exists, // or if request = full and non-skip/copy pair full request exists, // or if request = gran and non-skip/copy pair gran request exists, // then morph to copy of the matching parent list - // 2nd check: only if no match to 1st check + // 3rd check: only if no match to 1st check // if request = half and non-skip/copy pair full request exists, // then morph to half_from_full of the matching parent list // for 1st or 2nd check, parent can be copy list or pair or fix + int inewton,jnewton; + for (i = 0; i < nrequest; i++) { if (!requests[i]->fix && !requests[i]->compute) continue; for (j = 0; j < nrequest; j++) { - // Kokkos flags must match if (requests[i]->kokkos_device != requests[j]->kokkos_device) continue; if (requests[i]->kokkos_host != requests[j]->kokkos_host) continue; if (requests[i]->ssa != requests[j]->ssa) continue; - // newton 2 and newton 0 both are newton off - if ((requests[i]->newton & 2) != (requests[j]->newton & 2)) continue; + + // IJ newton = 1 for newton on, 2 for newton off + inewton = requests[i]->newton; + if (inewton == 0) inewton = force->newton_pair ? 1 : 2; + jnewton = requests[i]->newton; + if (jnewton == 0) jnewton = force->newton_pair ? 1 : 2; + if (inewton != jnewton) continue; if (requests[i]->half && requests[j]->pair && !requests[j]->skip && requests[j]->half && !requests[j]->copy) @@ -899,33 +910,37 @@ void Neighbor::init_pair() } // reorder plist vector if necessary - // relevant for lists that copy/skip/half-full from parent + // relevant for lists that are derived from a parent list: + // half-full,copy,skip // the child index must appear in plist after the parent index // swap two indices within plist when dependency is mis-ordered - // done when entire pass thru plist results in no swaps - + // start double loop check again whenever a swap is made + // done when entire double loop test results in no swaps + NeighList *ptr; int done = 0; while (!done) { done = 1; for (i = 0; i < npair_perpetual; i++) { - ptr = NULL; - if (lists[plist[i]]->listfull) ptr = lists[plist[i]]->listfull; - if (lists[plist[i]]->listcopy) ptr = lists[plist[i]]->listcopy; - // listskip check must be after listfull check - if (lists[plist[i]]->listskip) ptr = lists[plist[i]]->listskip; - if (ptr == NULL) continue; - for (m = 0; m < nrequest; m++) - if (ptr == lists[m]) break; - for (j = 0; j < npair_perpetual; j++) - if (m == plist[j]) break; - if (j < i) continue; - int tmp = plist[i]; // swap I,J indices - plist[i] = plist[j]; - plist[j] = tmp; - done = 0; - break; + for (k = 0; k < 3; k++) { + ptr = NULL; + if (k == 0) ptr = lists[plist[i]]->listcopy; + if (k == 1) ptr = lists[plist[i]]->listskip; + if (k == 2) ptr = lists[plist[i]]->listfull; + if (ptr == NULL) continue; + for (m = 0; m < nrequest; m++) + if (ptr == lists[m]) break; + for (j = 0; j < npair_perpetual; j++) + if (m == plist[j]) break; + if (j < i) continue; + int tmp = plist[i]; // swap I,J indices + plist[i] = plist[j]; + plist[j] = tmp; + done = 0; + break; + } + if (!done) break; } } @@ -1841,7 +1856,7 @@ void Neighbor::build_one(class NeighList *mylist, int preflag) // create stencil if hasn't been created since last setup_bins() call NStencil *ns = np->ns; - if (ns && ns->last_create < last_setup_bins) { + if (ns && ns->last_stencil < last_setup_bins) { ns->create_setup(); ns->create(); } @@ -1874,29 +1889,22 @@ void Neighbor::set(int narg, char **arg) /* ---------------------------------------------------------------------- reset timestamps in all NeignBin, NStencil, NPair classes so that neighbor lists will rebuild properly with timestep change + ditto for lastcall and last_setup_bins ------------------------------------------------------------------------- */ void Neighbor::reset_timestep(bigint ntimestep) { - for (int i = 0; i < nbin; i++) { - neigh_bin[i]->last_setup = -1; + for (int i = 0; i < nbin; i++) neigh_bin[i]->last_bin = -1; - neigh_bin[i]->last_bin_memory = -1; - } - - for (int i = 0; i < nstencil; i++) { - neigh_stencil[i]->last_create = -1; - neigh_stencil[i]->last_stencil_memory = -1; - neigh_stencil[i]->last_copy_bin = -1; - } - + for (int i = 0; i < nstencil; i++) + neigh_stencil[i]->last_stencil = -1; for (int i = 0; i < nlist; i++) { if (!neigh_pair[i]) continue; neigh_pair[i]->last_build = -1; - neigh_pair[i]->last_copy_bin_setup = -1; - neigh_pair[i]->last_copy_bin = -1; - neigh_pair[i]->last_copy_stencil = -1; } + + lastcall = -1; + last_setup_bins = -1; } /* ---------------------------------------------------------------------- diff --git a/src/npair.cpp b/src/npair.cpp index 3c38c40f50..bba3a2567f 100644 --- a/src/npair.cpp +++ b/src/npair.cpp @@ -28,8 +28,6 @@ NPair::NPair(LAMMPS *lmp) : Pointers(lmp), nb(NULL), ns(NULL), bins(NULL), stencil(NULL) { last_build = -1; - last_copy_bin_setup = last_copy_bin = last_copy_stencil = -1; - molecular = atom->molecular; } @@ -76,10 +74,10 @@ void NPair::copy_neighbor_info() } /* ---------------------------------------------------------------------- - copy bin geometry info from NBin class to this build class + copy info from NBin class to this build class ------------------------------------------------------------------------- */ -void NPair::copy_bin_setup_info() +void NPair::copy_bin_info() { nbinx = nb->nbinx; nbiny = nb->nbiny; @@ -95,20 +93,13 @@ void NPair::copy_bin_setup_info() bininvx = nb->bininvx; bininvy = nb->bininvy; bininvz = nb->bininvz; -} -/* ---------------------------------------------------------------------- - copy per-atom and per-bin vectors from NBin class to this build class -------------------------------------------------------------------------- */ - -void NPair::copy_bin_info() -{ bins = nb->bins; binhead = nb->binhead; } /* ---------------------------------------------------------------------- - copy needed info from NStencil class to this build class + copy info from NStencil class to this build class ------------------------------------------------------------------------- */ void NPair::copy_stencil_info() @@ -122,23 +113,15 @@ void NPair::copy_stencil_info() } /* ---------------------------------------------------------------------- - copy needed info from NStencil class to this build class + copy info from NBin and NStencil classes to this build class ------------------------------------------------------------------------- */ void NPair::build_setup() { - if (nb && last_copy_bin_setup < nb->last_setup) { - copy_bin_setup_info(); - last_copy_bin_setup = update->ntimestep; - } - if (nb && ((last_copy_bin < nb->last_bin_memory) || (bins != nb->bins))) { - copy_bin_info(); - last_copy_bin = update->ntimestep; - } - if (ns && ((last_copy_stencil < ns->last_create) || (stencil != ns->stencil))) { - copy_stencil_info(); - last_copy_stencil = update->ntimestep; - } + if (nb) copy_bin_info(); + if (ns) copy_stencil_info(); + + // set here, since build_setup() always called before build() last_build = update->ntimestep; } diff --git a/src/npair.h b/src/npair.h index a6440faddf..0c8dee5e2b 100644 --- a/src/npair.h +++ b/src/npair.h @@ -23,11 +23,7 @@ class NPair : protected Pointers { int istyle; // 1-N index into pairnames class NBin *nb; // ptr to NBin instance I depend on class NStencil *ns; // ptr to NStencil instance I depend on - bigint last_build; // last timestep build performed - bigint last_copy_bin_setup; // last timestep I invoked copy_bin_setup_info() - bigint last_copy_bin; // last step I invoked copy_bin_info() - bigint last_copy_stencil; // last step I invoked copy_bin_stencil_info() NPair(class LAMMPS *); virtual ~NPair() {} @@ -93,7 +89,6 @@ class NPair : protected Pointers { // methods for all NPair variants - void copy_bin_setup_info(); virtual void copy_bin_info(); virtual void copy_stencil_info(); diff --git a/src/nstencil.cpp b/src/nstencil.cpp index 68d5f80412..4ad2cbe56c 100644 --- a/src/nstencil.cpp +++ b/src/nstencil.cpp @@ -56,11 +56,9 @@ enum{NSQ,BIN,MULTI}; // also in Neighbor NStencil::NStencil(LAMMPS *lmp) : Pointers(lmp) { - last_create = last_stencil_memory = -1; - last_copy_bin = -1; + last_stencil = -1; xyzflag = 0; - maxstencil = maxstencil_multi = 0; stencil = NULL; stencilxyz = NULL; @@ -126,10 +124,8 @@ void NStencil::copy_bin_info() void NStencil::create_setup() { - if (nb && last_copy_bin < nb->last_setup) { - copy_bin_info(); - last_copy_bin = update->ntimestep; - } + if (nb) copy_bin_info(); + last_stencil = update->ntimestep; // sx,sy,sz = max range of stencil in each dim // smax = max possible size of entire 3d stencil @@ -157,7 +153,6 @@ void NStencil::create_setup() memory->destroy(stencilxyz); memory->create(stencilxyz,maxstencil,3,"neighstencil:stencilxyz"); } - last_stencil_memory = update->ntimestep; } } else { @@ -172,7 +167,6 @@ void NStencil::create_setup() stencil_multi[i] = NULL; distsq_multi[i] = NULL; } - last_stencil_memory = update->ntimestep; } if (smax > maxstencil_multi) { maxstencil_multi = smax; @@ -183,12 +177,9 @@ void NStencil::create_setup() "neighstencil:stencil_multi"); memory->create(distsq_multi[i],maxstencil_multi, "neighstencil:distsq_multi"); - last_stencil_memory = update->ntimestep; } } } - - last_create = update->ntimestep; } /* ---------------------------------------------------------------------- diff --git a/src/nstencil.h b/src/nstencil.h index 8672584a19..e71cf163c5 100644 --- a/src/nstencil.h +++ b/src/nstencil.h @@ -22,10 +22,7 @@ class NStencil : protected Pointers { public: int istyle; // 1-N index into binnames class NBin *nb; // ptr to NBin instance I depend on - - bigint last_create; // timesteps for last operations performed - bigint last_stencil_memory; - bigint last_copy_bin; + bigint last_stencil; // last timestep stencil was created int nstencil; // # of bins in stencil int *stencil; // list of bin offsets