diff --git a/doc/src/Howto_bioFF.rst b/doc/src/Howto_bioFF.rst index 72dcec2d31..92dd45b9b6 100644 --- a/doc/src/Howto_bioFF.rst +++ b/doc/src/Howto_bioFF.rst @@ -1,5 +1,5 @@ -CHARMM, AMBER, COMPASS, and DREIDING force fields -================================================= +CHARMM, AMBER, COMPASS, DREIDING, and OPLS force fields +======================================================= A compact summary of the concepts, definitions, and properties of force fields with explicit bonded interactions (like the ones discussed @@ -236,6 +236,40 @@ documentation for the formula it computes. * :doc:`special_bonds ` dreiding +OPLS +---- + +OPLS (Optimized Potentials for Liquid Simulations) is a general force +field for atomistic simulation of organic molecules in solvent. It was +developed by the `Jorgensen group +`_ at Purdue University and +later at Yale University. Multiple versions of the OPLS parameters +exist for united atom representations (OPLS-UA) and for all-atom +representations (OPLS-AA). + +This force field is based on atom types mapped to specific functional +groups in organic and biological molecules. Each atom includes a +static, partial atomic charge reflecting the oxidation state of the +element derived from its bonded neighbors :ref:`(Jorgensen) +` and computed based on increments determined by the +atom type of the atoms bond to it. + +The interaction styles listed below compute force field formulas that +are fully or in part consistent with the OPLS style force fields. See +each command's documentation for the formula it computes. Some are only +compatible with a subset of OPLS interactions. + +* :doc:`bond_style ` harmonic +* :doc:`angle_style ` harmonic +* :doc:`dihedral_style ` opls +* :doc:`improper_style ` cvff +* :doc:`improper_style ` fourier +* :doc:`improper_style ` harmonic +* :doc:`pair_style ` lj/cut/coul/cut +* :doc:`pair_style ` lj/cut/coul/long +* :doc:`pair_modify ` geometric +* :doc:`special_bonds ` lj/coul 0.0 0.0 0.5 + ---------- .. _Typelabel2: @@ -266,3 +300,6 @@ documentation for the formula it computes. **(Mayo)** Mayo, Olfason, Goddard III (1990). J Phys Chem, 94, 8897-8909. https://doi.org/10.1021/j100389a010 +.. _howto-Jorgensen: + +**(Jorgensen)** Jorgensen, Tirado-Rives (1988). J Am Chem Soc, 110, 1657-1666. https://doi.org/10.1021/ja00214a001 diff --git a/doc/src/Howto_chunk.rst b/doc/src/Howto_chunk.rst index f8655b745d..ea000eb22f 100644 --- a/doc/src/Howto_chunk.rst +++ b/doc/src/Howto_chunk.rst @@ -58,28 +58,30 @@ chunk ID for an individual atom can also be static (e.g. a molecule ID), or dynamic (e.g. what spatial bin an atom is in as it moves). Note that this compute allows the per-atom output of other -:doc:`computes `, :doc:`fixes `, and -:doc:`variables ` to be used to define chunk IDs for each -atom. This means you can write your own compute or fix to output a -per-atom quantity to use as chunk ID. See the :doc:`Modify ` -doc pages for info on how to do this. You can also define a :doc:`per-atom variable ` in the input script that uses a formula to -generate a chunk ID for each atom. +:doc:`computes `, :doc:`fixes `, and :doc:`variables +` to be used to define chunk IDs for each atom. This means +you can write your own compute or fix to output a per-atom quantity to +use as chunk ID. See the :doc:`Modify ` doc pages for info on +how to do this. You can also define a :doc:`per-atom variable +` in the input script that uses a formula to generate a chunk +ID for each atom. Fix ave/chunk command: ---------------------- -This fix takes the ID of a :doc:`compute chunk/atom ` command as input. For each chunk, -it then sums one or more specified per-atom values over the atoms in -each chunk. The per-atom values can be any atom property, such as -velocity, force, charge, potential energy, kinetic energy, stress, -etc. Additional keywords are defined for per-chunk properties like -density and temperature. More generally any per-atom value generated -by other :doc:`computes `, :doc:`fixes `, and :doc:`per-atom variables `, can be summed over atoms in each chunk. +This fix takes the ID of a :doc:`compute chunk/atom +` command as input. For each chunk, it then sums +one or more specified per-atom values over the atoms in each chunk. The +per-atom values can be any atom property, such as velocity, force, +charge, potential energy, kinetic energy, stress, etc. Additional +keywords are defined for per-chunk properties like density and +temperature. More generally any per-atom value generated by other +:doc:`computes `, :doc:`fixes `, and :doc:`per-atom +variables `, can be summed over atoms in each chunk. Similar to other averaging fixes, this fix allows the summed per-chunk values to be time-averaged in various ways, and output to a file. The -fix produces a global array as output with one row of values per -chunk. +fix produces a global array as output with one row of values per chunk. Compute \*/chunk commands: -------------------------- @@ -97,17 +99,20 @@ category: * :doc:`compute torque/chunk ` * :doc:`compute vcm/chunk ` -They each take the ID of a :doc:`compute chunk/atom ` command as input. As their names -indicate, they calculate the center-of-mass, radius of gyration, -moments of inertia, mean-squared displacement, temperature, torque, -and velocity of center-of-mass for each chunk of atoms. The :doc:`compute property/chunk ` command can tally the -count of atoms in each chunk and extract other per-chunk properties. +They each take the ID of a :doc:`compute chunk/atom +` command as input. As their names indicate, they +calculate the center-of-mass, radius of gyration, moments of inertia, +mean-squared displacement, temperature, torque, and velocity of +center-of-mass for each chunk of atoms. The :doc:`compute +property/chunk ` command can tally the count of +atoms in each chunk and extract other per-chunk properties. -The reason these various calculations are not part of the :doc:`fix ave/chunk command `, is that each requires a more +The reason these various calculations are not part of the :doc:`fix +ave/chunk command `, is that each requires a more complicated operation than simply summing and averaging over per-atom -values in each chunk. For example, many of them require calculation -of a center of mass, which requires summing mass\*position over the -atoms and then dividing by summed mass. +values in each chunk. For example, many of them require calculation of +a center of mass, which requires summing mass\*position over the atoms +and then dividing by summed mass. All of these computes produce a global vector or global array as output, with one or more values per chunk. The output can be used in @@ -118,9 +123,10 @@ various ways: * As input to the :doc:`fix ave/histo ` command to histogram values across chunks. E.g. a histogram of cluster sizes or molecule diffusion rates. -* As input to special functions of :doc:`equal-style variables `, like sum() and max() and ave(). E.g. to - find the largest cluster or fastest diffusing molecule or average - radius-of-gyration of a set of molecules (chunks). +* As input to special functions of :doc:`equal-style variables + `, like sum() and max() and ave(). E.g. to find the largest + cluster or fastest diffusing molecule or average radius-of-gyration of + a set of molecules (chunks). Other chunk commands: --------------------- @@ -138,9 +144,10 @@ spatially average per-chunk values calculated by a per-chunk compute. The :doc:`compute reduce/chunk ` command reduces a peratom value across the atoms in each chunk to produce a value per -chunk. When used with the :doc:`compute chunk/spread/atom ` command it can -create peratom values that induce a new set of chunks with a second -:doc:`compute chunk/atom ` command. +chunk. When used with the :doc:`compute chunk/spread/atom +` command it can create peratom values that +induce a new set of chunks with a second :doc:`compute chunk/atom +` command. Example calculations with chunks -------------------------------- diff --git a/doc/src/compute_angle_local.rst b/doc/src/compute_angle_local.rst index 126dd309af..d4491c6945 100644 --- a/doc/src/compute_angle_local.rst +++ b/doc/src/compute_angle_local.rst @@ -78,7 +78,7 @@ system and output the statistics in various ways: compute 2 all angle/local eng theta v_cos v_cossq set theta t dump 1 all local 100 tmp.dump c_1[*] c_2[*] - compute 3 all reduce ave c_2[*] + compute 3 all reduce ave c_2[*] inputs local thermo_style custom step temp press c_3[*] fix 10 all ave/histo 10 10 100 -1 1 20 c_2[3] mode vector file tmp.histo diff --git a/doc/src/compute_bond_local.rst b/doc/src/compute_bond_local.rst index 5036358c8c..e070d507b1 100644 --- a/doc/src/compute_bond_local.rst +++ b/doc/src/compute_bond_local.rst @@ -139,7 +139,7 @@ output the statistics in various ways: compute 2 all bond/local engpot dist v_dsq set dist d dump 1 all local 100 tmp.dump c_1[*] c_2[*] - compute 3 all reduce ave c_2[*] + compute 3 all reduce ave c_2[*] inputs local thermo_style custom step temp press c_3[*] fix 10 all ave/histo 10 10 100 0 6 20 c_2[3] mode vector file tmp.histo diff --git a/doc/src/compute_dihedral_local.rst b/doc/src/compute_dihedral_local.rst index 291b870373..d809cd39ce 100644 --- a/doc/src/compute_dihedral_local.rst +++ b/doc/src/compute_dihedral_local.rst @@ -76,7 +76,7 @@ angle in the system and output the statistics in various ways: compute 2 all dihedral/local phi v_cos v_cossq set phi p dump 1 all local 100 tmp.dump c_1[*] c_2[*] - compute 3 all reduce ave c_2[*] + compute 3 all reduce ave c_2[*] inputs local thermo_style custom step temp press c_3[*] fix 10 all ave/histo 10 10 100 -1 1 20 c_2[2] mode vector file tmp.histo diff --git a/doc/src/compute_reduce.rst b/doc/src/compute_reduce.rst index 9eafe7cd5a..e5c99a478f 100644 --- a/doc/src/compute_reduce.rst +++ b/doc/src/compute_reduce.rst @@ -206,11 +206,13 @@ IDs and the bond stretch will be printed with thermodynamic output. The *inputs* keyword allows selection of whether all the inputs are per-atom or local quantities. As noted above, all the inputs must be -the same kind (per-atom or local). Per-atom is the default setting. -If a compute or fix is specified as an input, it must produce per-atom -or local data to match this setting. If it produces both, e.g. for +the same kind (per-atom or local). Per-atom is the default setting. If +a compute or fix is specified as an input, it must produce per-atom or +local data to match this setting. If it produces both, like for example the :doc:`compute voronoi/atom ` command, then -this keyword selects between them. +this keyword selects between them. If a compute *only* produces local +data, like for example the :doc:`compute bond/local command +`, the setting "inputs local" is *required*. ---------- diff --git a/doc/src/compute_reduce_chunk.rst b/doc/src/compute_reduce_chunk.rst index 8ec19ade66..eeadd50a97 100644 --- a/doc/src/compute_reduce_chunk.rst +++ b/doc/src/compute_reduce_chunk.rst @@ -37,55 +37,57 @@ Description Define a calculation that reduces one or more per-atom vectors into per-chunk values. This can be useful for diagnostic output. Or when -used in conjunction with the :doc:`compute chunk/spread/atom ` command it can be -used to create per-atom values that induce a new set of chunks with a -second :doc:`compute chunk/atom ` command. An -example is given below. +used in conjunction with the :doc:`compute chunk/spread/atom +` command it can be used to create per-atom +values that induce a new set of chunks with a second :doc:`compute +chunk/atom ` command. An example is given below. -In LAMMPS, chunks are collections of atoms defined by a :doc:`compute chunk/atom ` command, which assigns each atom -to a single chunk (or no chunk). The ID for this command is specified -as chunkID. For example, a single chunk could be the atoms in a -molecule or atoms in a spatial bin. See the :doc:`compute chunk/atom ` and :doc:`Howto chunk ` -doc pages for details of how chunks can be defined and examples of how -they can be used to measure properties of a system. +In LAMMPS, chunks are collections of atoms defined by a :doc:`compute +chunk/atom ` command, which assigns each atom to a +single chunk (or no chunk). The ID for this command is specified as +chunkID. For example, a single chunk could be the atoms in a molecule +or atoms in a spatial bin. See the :doc:`compute chunk/atom +` and :doc:`Howto chunk ` doc pages for +details of how chunks can be defined and examples of how they can be +used to measure properties of a system. For each atom, this compute accesses its chunk ID from the specified -*chunkID* compute. The per-atom value from an input contributes -to a per-chunk value corresponding the the chunk ID. +*chunkID* compute. The per-atom value from an input contributes to a +per-chunk value corresponding the chunk ID. The reduction operation is specified by the *mode* setting and is performed over all the per-atom values from the atoms in each chunk. -The *sum* option adds the pre-atom values to a per-chunk total. The -*min* or *max* options find the minimum or maximum value of the -per-atom values for each chunk. +The *sum* option adds the per-atom values to a per-chunk total. The +*min* or *max* options find the minimum or maximum value of the per-atom +values for each chunk. -Note that only atoms in the specified group contribute to the -reduction operation. If the *chunkID* compute returns a 0 for the -chunk ID of an atom (i.e., the atom is not in a chunk defined by the -:doc:`compute chunk/atom ` command), that atom will -also not contribute to the reduction operation. An input that is a -compute or fix may define its own group which affects the quantities -it returns. For example, a compute with return a zero value for atoms -that are not in the group specified for that compute. +Note that only atoms in the specified group contribute to the reduction +operation. If the *chunkID* compute returns a 0 for the chunk ID of an +atom (i.e., the atom is not in a chunk defined by the :doc:`compute +chunk/atom ` command), that atom will also not +contribute to the reduction operation. An input that is a compute or +fix may define its own group which affects the quantities it returns. +For example, a compute will return a zero value for atoms that are not +in the group specified for that compute. Each listed input is operated on independently. Each input can be the -result of a :doc:`compute ` or :doc:`fix ` or the evaluation -of an atom-style :doc:`variable `. +result of a :doc:`compute ` or :doc:`fix ` or the +evaluation of an atom-style :doc:`variable `. -Note that for values from a compute or fix, the bracketed index I can -be specified using a wildcard asterisk with the index to effectively +Note that for values from a compute or fix, the bracketed index I can be +specified using a wildcard asterisk with the index to effectively specify multiple values. This takes the form "\*" or "\*n" or "m\*" or -"m\*n". If :math:`N` is the size of the vector (for *mode* = scalar) or the -number of columns in the array (for *mode* = vector), then an asterisk -with no numeric values means all indices from 1 to :math:`N`. A leading -asterisk means all indices from 1 to n (inclusive). A trailing -asterisk means all indices from n to :math:`N` (inclusive). A middle asterisk -means all indices from m to n (inclusive). +"m\*n". If :math:`N` is the size of the vector (for *mode* = scalar) or +the number of columns in the array (for *mode* = vector), then an +asterisk with no numeric values means all indices from 1 to :math:`N`. +A leading asterisk means all indices from 1 to n (inclusive). A +trailing asterisk means all indices from n to :math:`N` (inclusive). A +middle asterisk means all indices from m to n (inclusive). Using a wildcard is the same as if the individual columns of the array -had been listed one by one. For example, the following two compute reduce/chunk -commands are equivalent, since the -:doc:`compute property/chunk ` command creates a per-atom +had been listed one by one. For example, the following two compute +reduce/chunk commands are equivalent, since the :doc:`compute +property/chunk ` command creates a per-atom array with 3 columns: .. code-block:: LAMMPS @@ -164,13 +166,14 @@ Output info """"""""""" This compute calculates a global vector if a single input value is -specified, otherwise a global array is output. The number of columns -in the array is the number of inputs provided. The length of the -vector or the number of vector elements or array rows = the number of -chunks *Nchunk* as calculated by the specified :doc:`compute chunk/atom ` command. The vector or array can -be accessed by any command that uses global values from a compute as -input. See the :doc:`Howto output ` page for an -overview of LAMMPS output options. +specified, otherwise a global array is output. The number of columns in +the array is the number of inputs provided. The length of the vector or +the number of vector elements or array rows = the number of chunks +*Nchunk* as calculated by the specified :doc:`compute chunk/atom +` command. The vector or array can be accessed by +any command that uses global values from a compute as input. See the +:doc:`Howto output ` page for an overview of LAMMPS output +options. The per-atom values for the vector or each column of the array will be in whatever :doc:`units ` the corresponding input value is in. @@ -183,7 +186,9 @@ Restrictions Related commands """""""""""""""" -:doc:`compute chunk/atom `, :doc:`compute reduce `, :doc:`compute chunk/spread/atom ` +:doc:`compute chunk/atom `, +:doc:`compute reduce `, +:doc:`compute chunk/spread/atom ` Default """"""" diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index b4779b8bf3..31c81b5df7 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -129,6 +129,9 @@ package ` doc page on for more info. The method is implemented for orthogonal simulation boxes whose size does not change in time, and axis-aligned planes. +Contributions from bonds, angles, and dihedrals are not compatible +with MPI parallel runs. + The method only works with two-body pair interactions, because it requires the class method ``Pair::single()`` to be implemented, which is not possible for manybody potentials. In particular, compute diff --git a/doc/src/fix_atom_swap.rst b/doc/src/fix_atom_swap.rst index 48f5b33a74..fd4a2f7245 100644 --- a/doc/src/fix_atom_swap.rst +++ b/doc/src/fix_atom_swap.rst @@ -119,15 +119,14 @@ groups of atoms that have different charges, these charges will not be changed when the atom types change. Since this fix computes total potential energies before and after -proposed swaps, so even complicated potential energy calculations are -OK, including the following: +proposed swaps, even complicated potential energy calculations are +acceptable, including the following: * long-range electrostatics (:math:`k`-space) * many body pair styles -* hybrid pair styles -* eam pair styles +* hybrid pair styles (with restrictions) +* EAM pair styles * triclinic systems -* need to include potential energy contributions from other fixes Some fixes have an associated potential energy. Examples of such fixes include: :doc:`efield `, :doc:`gravity `, @@ -181,6 +180,10 @@ This fix is part of the MC package. It is only enabled if LAMMPS was built with that package. See the :doc:`Build package ` doc page for more info. +When this fix is used with a :doc:`hybrid pair style ` +system, only swaps between atom types of the same sub-style (or +combination of sub-styles) are permitted. + This fix cannot be used with systems that do not have per-type masses (e.g. atom style sphere) since the implemented algorithm pre-computes velocity rescaling factors from per-type masses and ignores any per-atom diff --git a/doc/src/fix_halt.rst b/doc/src/fix_halt.rst index 25331804aa..0bcf2fb5ea 100644 --- a/doc/src/fix_halt.rst +++ b/doc/src/fix_halt.rst @@ -101,7 +101,7 @@ hstyle = bondmax option. .. code-block:: LAMMPS compute bdist all bond/local dist - compute bmax all reduce max c_bdist + compute bmax all reduce max c_bdist inputs local variable bondmax equal c_bmax Thus these two versions of a fix halt command will do the same thing: diff --git a/doc/src/fix_langevin.rst b/doc/src/fix_langevin.rst index e04805427e..30e4c48270 100644 --- a/doc/src/fix_langevin.rst +++ b/doc/src/fix_langevin.rst @@ -231,12 +231,6 @@ the particles. As described below, this energy can then be printed out or added to the potential energy of the system to monitor energy conservation. -.. note:: - - This accumulated energy does NOT include kinetic energy removed - by the *zero* flag. LAMMPS will print a warning when both options are - active. - The keyword *zero* can be used to eliminate drift due to the thermostat. Because the random forces on different atoms are independent, they do not sum exactly to zero. As a result, this fix diff --git a/doc/src/group.rst b/doc/src/group.rst index 15ab0c9dc8..a7a29467ff 100644 --- a/doc/src/group.rst +++ b/doc/src/group.rst @@ -162,7 +162,7 @@ potential energy is above the threshold value :math:`-3.0`. compute 1 all pe/atom compute 2 all reduce sum c_1 thermo_style custom step temp pe c_2 - run 0 + run 0 post no variable eatom atom "c_1 > -3.0" group hienergy variable eatom @@ -173,7 +173,7 @@ Note that these lines compute 2 all reduce sum c_1 thermo_style custom step temp pe c_2 - run 0 + run 0 post no are necessary to ensure that the "eatom" variable is current when the group command invokes it. Because the eatom variable computes the diff --git a/doc/utils/requirements.txt b/doc/utils/requirements.txt index 4b9e1a4e48..acf575fe58 100644 --- a/doc/utils/requirements.txt +++ b/doc/utils/requirements.txt @@ -1,4 +1,4 @@ -Sphinx >= 5.3.0, <9.0 +Sphinx >= 5.3.0, <8.2.0 sphinxcontrib-spelling sphinxcontrib-jquery sphinx-design diff --git a/doc/utils/sphinx-config/false_positives.txt b/doc/utils/sphinx-config/false_positives.txt index 2754e4bb1c..0708a14882 100644 --- a/doc/utils/sphinx-config/false_positives.txt +++ b/doc/utils/sphinx-config/false_positives.txt @@ -1286,6 +1286,7 @@ gcc gcmc gdot GeC +Geesthacht Geier gencode Geocomputing @@ -4219,6 +4220,7 @@ zeeman Zeeman Zemer zenodo +Zentrum Zepeda zflag Zhang diff --git a/src/DIELECTRIC/pppm_disp_dielectric.cpp b/src/DIELECTRIC/pppm_disp_dielectric.cpp index 346df8bf93..aba8ccf486 100644 --- a/src/DIELECTRIC/pppm_disp_dielectric.cpp +++ b/src/DIELECTRIC/pppm_disp_dielectric.cpp @@ -715,7 +715,7 @@ void PPPMDispDielectric::fieldforce_c_ad() { int i,l,m,n,nx,ny,nz,mx,my,mz; FFT_SCALAR dx,dy,dz; - FFT_SCALAR ekx,eky,ekz,u; + FFT_SCALAR ekx,eky,ekz; double s1,s2,s3; double sf = 0.0; diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index b8d21d2a4f..7c66c1fe50 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -39,7 +39,7 @@ using namespace LAMMPS_NS; -static constexpr double SMALL = 0.001; +static constexpr double SMALL = 0.001; enum { X, Y, Z }; enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE, DIHEDRAL }; @@ -63,7 +63,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( } else if (strcmp(arg[3], "z") == 0) { dir = Z; } else - error->all(FLERR, "Illegal compute stress/mop command"); + error->all(FLERR, "Illegal compute stress/mop plane direction {}", arg[3]); // Position of the plane @@ -89,7 +89,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( error->all(FLERR, "Plane for compute stress/mop is out of bounds"); } - if (pos < (domain->boxlo[dir] + domain->prd_half[dir])) { + if (pos < (domain->boxlo[dir] + domain->prd_half[dir])) { pos1 = pos + domain->prd[dir]; } else { pos1 = pos - domain->prd[dir]; @@ -133,27 +133,17 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( which[nvalues] = ANGLE; nvalues++; } - } else if (strcmp(arg[iarg],"dihedral") == 0) { - for (i=0; i<3; i++) { + } else if (strcmp(arg[iarg], "dihedral") == 0) { + for (i = 0; i < 3; i++) { which[nvalues] = DIHEDRAL; nvalues++; } } else - error->all(FLERR, "Illegal compute stress/mop command"); //break; + error->all(FLERR, "Unknown compute stress/mop keyword {}", arg[iarg]); iarg++; } - // Error checks: - - // orthogonal simulation box - if (domain->triclinic != 0) - error->all(FLERR, "Compute stress/mop is incompatible with triclinic simulation box"); - - // 2D and pressure calculation in the Z coordinate - if (domain->dimension == 2 && dir == Z) - error->all(FLERR, "Compute stress/mop is incompatible with Z in 2d system"); - // Initialize some variables values_local = values_global = vector = nullptr; @@ -173,8 +163,8 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( memory->create(bond_global, nvalues, "stress/mop:bond_global"); memory->create(angle_local, nvalues, "stress/mop:angle_local"); memory->create(angle_global, nvalues, "stress/mop:angle_global"); - memory->create(dihedral_local,nvalues,"stress/mop:dihedral_local"); - memory->create(dihedral_global,nvalues,"stress/mop:dihedral_global"); + memory->create(dihedral_local, nvalues, "stress/mop:dihedral_local"); + memory->create(dihedral_global, nvalues, "stress/mop:dihedral_global"); size_vector = nvalues; vector_flag = 1; @@ -202,6 +192,16 @@ ComputeStressMop::~ComputeStressMop() void ComputeStressMop::init() { + // Error checks: + + // 2D and pressure calculation in the Z coordinate + if (domain->dimension == 2 && dir == Z) + error->all(FLERR, "Compute stress/mop is incompatible with Z in 2d system"); + + // orthogonal simulation box + + if (domain->triclinic != 0) + error->all(FLERR, "Compute stress/mop is incompatible with triclinic simulation box"); // Conversion constants @@ -224,11 +224,12 @@ void ComputeStressMop::init() dt = update->dt; - // Error check + // Error checks: // Compute stress/mop requires fixed simulation box + if (domain->box_change_size || domain->box_change_shape || domain->deform_flag) - error->all(FLERR, "Compute stress/mop requires a fixed size simulation box"); + error->all(FLERR, "Compute stress/mop requires a fixed simulation box geometry"); // This compute requires a pair style with pair_single method implemented @@ -242,34 +243,46 @@ void ComputeStressMop::init() // issue an error for unimplemented intramolecular potentials or Kspace. - if (force->bond) bondflag = 1; + if (force->bond) { + bondflag = 1; + if (comm->nprocs > 1) + error->one(FLERR, "compute stress/mop with bonds does not (yet) support MPI parallel runs"); + } + if (force->angle) { if (force->angle->born_matrix_enable == 0) { if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0)) - error->all(FLERR, "compute stress/mop does not account for angle potentials"); + error->one(FLERR, "compute stress/mop does not account for angle potentials"); } else { angleflag = 1; + if (comm->nprocs > 1) + error->one(FLERR, + "compute stress/mop with angles does not (yet) support MPI parallel runs"); } } if (force->dihedral) { if (force->dihedral->born_matrix_enable == 0) { if ((strcmp(force->dihedral_style, "zero") != 0) && (strcmp(force->dihedral_style, "none") != 0)) - error->all(FLERR, "compute stress/mop does not account for dihedral potentials"); + error->one(FLERR, "compute stress/mop does not account for dihedral potentials"); } else { dihedralflag = 1; + if (comm->nprocs > 1) + error->one(FLERR, + "compute stress/mop with dihedrals does not (yet) support MPI parallel runs"); } } if (force->improper) { if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) - error->all(FLERR, "compute stress/mop does not account for improper potentials"); + error->one(FLERR, "compute stress/mop does not account for improper potentials"); } if (force->kspace) error->warning(FLERR, "compute stress/mop does not account for kspace contributions"); } // need an occasional half neighbor list + neighbor->add_request(this, NeighConst::REQ_OCCASIONAL); } @@ -324,11 +337,11 @@ void ComputeStressMop::compute_vector() //Compute dihedral contribution on separate procs compute_dihedrals(); } else { - for (int i=0; i pos)) || - ((xi[dir] < pos1) && (xj[dir] > pos1))) { + ((xi[dir] < pos1) && (xj[dir] > pos1))) { pair->single(i, j, itype, jtype, rsq, factor_coul, factor_lj, fpair); values_local[m] -= fpair * (xi[0] - xj[0]) / area * nktv2p; values_local[m + 1] -= fpair * (xi[1] - xj[1]) / area * nktv2p; @@ -869,9 +882,7 @@ void ComputeStressMop::compute_dihedrals() double x_atom_4[3] = {0.0, 0.0, 0.0}; // initialization - for (int i = 0; i < nvalues; i++) { - dihedral_local[i] = 0.0; - } + for (int i = 0; i < nvalues; i++) { dihedral_local[i] = 0.0; } double local_contribution[3] = {0.0, 0.0, 0.0}; for (atom2 = 0; atom2 < nlocal; atom2++) { @@ -889,9 +900,9 @@ void ComputeStressMop::compute_dihedrals() for (i = 0; i < nd; i++) { if (molecular == 1) { if (tag[atom2] != dihedral_atom2[atom2][i]) continue; - atom1 = atom->map(dihedral_atom1[atom2][i]); - atom3 = atom->map(dihedral_atom3[atom2][i]); - atom4 = atom->map(dihedral_atom4[atom2][i]); + atom1 = atom->map(dihedral_atom1[atom2][i]); + atom3 = atom->map(dihedral_atom3[atom2][i]); + atom4 = atom->map(dihedral_atom4[atom2][i]); } else { if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; tagprev = tag[atom2] - iatom - 1; @@ -943,9 +954,9 @@ void ComputeStressMop::compute_dihedrals() double tau_right = (x_atom_2[dir] - pos) / (x_atom_2[dir] - x_atom_1[dir]); double tau_middle = (x_atom_3[dir] - pos) / (x_atom_3[dir] - x_atom_2[dir]); double tau_left = (x_atom_4[dir] - pos) / (x_atom_4[dir] - x_atom_3[dir]); - bool right_cross = ((tau_right >=0) && (tau_right <= 1)); + bool right_cross = ((tau_right >= 0) && (tau_right <= 1)); bool middle_cross = ((tau_middle >= 0) && (tau_middle <= 1)); - bool left_cross = ((tau_left >=0) && (tau_left <= 1)); + bool left_cross = ((tau_left >= 0) && (tau_left <= 1)); // no bonds crossing the plane if (!right_cross && !middle_cross && !left_cross) continue; @@ -972,45 +983,45 @@ void ComputeStressMop::compute_dihedrals() vb3z = x_atom_4[2] - x_atom_3[2]; // c0 calculation - sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); - sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); - sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z); + sb2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z); + sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z); rb1 = sqrt(sb1); rb3 = sqrt(sb3); - c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3; // 1st and 2nd angle - b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z; b1mag = sqrt(b1mag2); - b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; b2mag = sqrt(b2mag2); - b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z; b3mag = sqrt(b3mag2); - ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; - r12c1 = 1.0 / (b1mag*b2mag); + ctmp = vb1x * vb2x + vb1y * vb2y + vb1z * vb2z; + r12c1 = 1.0 / (b1mag * b2mag); c1mag = ctmp * r12c1; - ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; - r12c2 = 1.0 / (b2mag*b3mag); + ctmp = vb2xm * vb3x + vb2ym * vb3y + vb2zm * vb3z; + r12c2 = 1.0 / (b2mag * b3mag); c2mag = ctmp * r12c2; // cos and sin of 2 angles and final c - sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sin2 = MAX(1.0 - c1mag * c1mag, 0.0); sc1 = sqrt(sin2); if (sc1 < SMALL) sc1 = SMALL; - sc1 = 1.0/sc1; + sc1 = 1.0 / sc1; - sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sin2 = MAX(1.0 - c2mag * c2mag, 0.0); sc2 = sqrt(sin2); if (sc2 < SMALL) sc2 = SMALL; - sc2 = 1.0/sc2; + sc2 = 1.0 / sc2; s1 = sc1 * sc1; s2 = sc2 * sc2; s12 = sc1 * sc2; - c = (c0 + c1mag*c2mag) * s12; + c = (c0 + c1mag * c2mag) * s12; // error check if (c > 1.0) c = 1.0; @@ -1020,28 +1031,28 @@ void ComputeStressMop::compute_dihedrals() double a = dudih; c = c * a; s12 = s12 * a; - a11 = c*sb1*s1; - a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); - a33 = c*sb3*s2; - a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12); - a13 = -rb1*rb3*s12; - a23 = r12c2 * (c2mag*c*s2 + c1mag*s12); + a11 = c * sb1 * s1; + a22 = -sb2 * (2.0 * c0 * s12 - c * (s1 + s2)); + a33 = c * sb3 * s2; + a12 = -r12c1 * (c1mag * c * s1 + c2mag * s12); + a13 = -rb1 * rb3 * s12; + a23 = r12c2 * (c2mag * c * s2 + c1mag * s12); - sx2 = a12*vb1x + a22*vb2x + a23*vb3x; - sy2 = a12*vb1y + a22*vb2y + a23*vb3y; - sz2 = a12*vb1z + a22*vb2z + a23*vb3z; + sx2 = a12 * vb1x + a22 * vb2x + a23 * vb3x; + sy2 = a12 * vb1y + a22 * vb2y + a23 * vb3y; + sz2 = a12 * vb1z + a22 * vb2z + a23 * vb3z; - f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; - f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; - f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; + f1[0] = a11 * vb1x + a12 * vb2x + a13 * vb3x; + f1[1] = a11 * vb1y + a12 * vb2y + a13 * vb3y; + f1[2] = a11 * vb1z + a12 * vb2z + a13 * vb3z; f2[0] = -sx2 - f1[0]; f2[1] = -sy2 - f1[1]; f2[2] = -sz2 - f1[2]; - f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; - f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; - f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; + f4[0] = a13 * vb1x + a23 * vb2x + a33 * vb3x; + f4[1] = a13 * vb1y + a23 * vb2y + a33 * vb3y; + f4[2] = a13 * vb1z + a23 * vb2z + a33 * vb3z; f3[0] = sx2 - f4[0]; f3[1] = sy2 - f4[1]; @@ -1106,9 +1117,9 @@ void ComputeStressMop::compute_dihedrals() df[1] = sgn * (f1[1] + f3[1]); df[2] = sgn * (f1[2] + f3[2]); } - local_contribution[0] += df[0]/area*nktv2p; - local_contribution[1] += df[1]/area*nktv2p; - local_contribution[2] += df[2]/area*nktv2p; + local_contribution[0] += df[0] / area * nktv2p; + local_contribution[1] += df[1] / area * nktv2p; + local_contribution[2] += df[2] / area * nktv2p; } } @@ -1116,11 +1127,10 @@ void ComputeStressMop::compute_dihedrals() int m = 0; while (m < nvalues) { if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == DIHEDRAL)) { - dihedral_local[m] = local_contribution[0]; - dihedral_local[m+1] = local_contribution[1]; - dihedral_local[m+2] = local_contribution[2]; + dihedral_local[m] = local_contribution[0]; + dihedral_local[m + 1] = local_contribution[1]; + dihedral_local[m + 2] = local_contribution[2]; } m += 3; } - } diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index ca2d095fd9..c8087b60a9 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -39,7 +39,7 @@ using namespace LAMMPS_NS; -static constexpr double SMALL = 0.001; +static constexpr double SMALL = 0.001; enum { X, Y, Z }; enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE, DIHEDRAL }; @@ -64,7 +64,7 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a } else if (strcmp(arg[3], "z") == 0) { dir = Z; } else - error->all(FLERR, "Illegal compute stress/mop/profile command"); + error->all(FLERR, "Illegal compute stress/mop/profile plane direction {}", arg[3]); // bin parameters @@ -118,30 +118,18 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a which[nvalues] = ANGLE; nvalues++; } - } else if (strcmp(arg[iarg],"dihedral") == 0) { - for (i=0; i<3; i++) { + } else if (strcmp(arg[iarg], "dihedral") == 0) { + for (i = 0; i < 3; i++) { which[nvalues] = DIHEDRAL; nvalues++; } } else - error->all(FLERR, "Illegal compute stress/mop/profile command"); //break; + error->all(FLERR, "Unknown compute stress/mop/profile keyword {}", arg[iarg]); iarg++; } - // check domain related errors - - // 3D only - - if (domain->dimension == 2 && dir == Z) - error->all(FLERR, "Compute stress/mop/profile incompatible with Z in 2d system"); - - // orthogonal simulation box - - if (domain->triclinic != 0) - error->all(FLERR, "Compute stress/mop/profile incompatible with triclinic simulation box"); - - // initialize some variables + // Initialize some variables nbins = 0; coord = coordp = nullptr; @@ -192,12 +180,25 @@ ComputeStressMopProfile::~ComputeStressMopProfile() void ComputeStressMopProfile::init() { - // conversion constants + // check domain related errors + + // 3D only + + if (domain->dimension == 2 && dir == Z) + error->all(FLERR, "Compute stress/mop/profile is incompatible with Z in 2d system"); + + // check for orthogonal simulation box + + if (domain->triclinic != 0) + error->all(FLERR, "Compute stress/mop/profile is incompatible with triclinic simulation box"); + + // Conversion constants nktv2p = force->nktv2p; ftm2v = force->ftm2v; - // plane area + // Plane area + if (domain->dimension == 3) { area = 1; int i; @@ -208,7 +209,7 @@ void ComputeStressMopProfile::init() area = (dir == X) ? domain->prd[1] : domain->prd[0]; } - // timestep Value + // Timestep Value dt = update->dt; @@ -217,46 +218,59 @@ void ComputeStressMopProfile::init() // Compute stress/mop/profile requires fixed simulation box if (domain->box_change_size || domain->box_change_shape || domain->deform_flag) - error->all(FLERR, "Compute stress/mop/profile requires a fixed simulation box"); + error->all(FLERR, "Compute stress/mop/profile requires a fixed simulation box geometry"); // This compute requires a pair style with pair_single method implemented if (!force->pair) error->all(FLERR, "No pair style is defined for compute stress/mop/profile"); if (force->pair->single_enable == 0) - error->all(FLERR, "Pair style does not support compute stress/mop/profile"); + error->all(FLERR, "Pair style {} does not support compute stress/mop/profile", + force->pair_style); // Errors if (comm->me == 0) { - // Compute stress/mop/profile only accounts for pair interactions. - // issue an error if any intramolecular potential or Kspace is defined. + // issue an error for unimplemented intramolecular potentials or Kspace. - if (force->bond) bondflag = 1; + if (force->bond) { + bondflag = 1; + if (comm->nprocs > 1) + error->one( + FLERR, + "compute stress/mop/profile with bonds does not (yet) support MPI parallel runs"); + } if (force->angle) { if (force->angle->born_matrix_enable == 0) { if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0)) - error->all(FLERR,"compute stress/mop/profile does not account for angle potentials"); + error->one(FLERR, "compute stress/mop/profile does not account for angle potentials"); } else { angleflag = 1; + if (comm->nprocs > 1) + error->one( + FLERR, + "compute stress/mop/profile with angles does not (yet) support MPI parallel runs"); } } - if (force->dihedral) { if (force->dihedral->born_matrix_enable == 0) { if ((strcmp(force->dihedral_style, "zero") != 0) && (strcmp(force->dihedral_style, "none") != 0)) - error->all(FLERR, "compute stress/mop/profile does not account for dihedral potentials"); + error->one(FLERR, "compute stress/mop/profile does not account for dihedral potentials"); } else { dihedralflag = 1; + if (comm->nprocs > 1) + error->one( + FLERR, + "compute stress/mop/profile with dihedrals does not (yet) support MPI parallel runs"); } } - - if (force->improper) + if (force->improper) { if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) - error->all(FLERR, "compute stress/mop/profile does not account for improper potentials"); + error->one(FLERR, "compute stress/mop/profile does not account for improper potentials"); + } if (force->kspace) error->warning(FLERR, "compute stress/mop/profile does not account for kspace contributions"); } @@ -308,28 +322,26 @@ void ComputeStressMopProfile::compute_array() compute_angles(); } else { for (int m = 0; m < nbins; m++) { - for (int i = 0; i < nvalues; i++) { - angle_local[m][i] = 0.0; - } + for (int i = 0; i < nvalues; i++) { angle_local[m][i] = 0.0; } } } // sum angle contribution over all procs - MPI_Allreduce(&angle_local[0][0],&angle_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&angle_local[0][0], &angle_global[0][0], nbins * nvalues, MPI_DOUBLE, MPI_SUM, + world); if (dihedralflag) { //Compute dihedral contribution on separate procs compute_dihedrals(); } else { for (int m = 0; m < nbins; m++) { - for (int i = 0; i < nvalues; i++) { - dihedral_local[m][i] = 0.0; - } + for (int i = 0; i < nvalues; i++) { dihedral_local[m][i] = 0.0; } } } // sum dihedral contribution over all procs - MPI_Allreduce(&dihedral_local[0][0],&dihedral_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&dihedral_local[0][0], &dihedral_global[0][0], nbins * nvalues, MPI_DOUBLE, MPI_SUM, + world); for (int ibin = 0; ibin < nbins; ibin++) { array[ibin][0] = coord[ibin]; @@ -337,7 +349,8 @@ void ComputeStressMopProfile::compute_array() int mo = 1; int m = 0; while (m < nvalues) { - array[ibin][m + mo] = values_global[ibin][m] + bond_global[ibin][m] + angle_global[ibin][m] + dihedral_global[ibin][m]; + array[ibin][m + mo] = values_global[ibin][m] + bond_global[ibin][m] + angle_global[ibin][m] + + dihedral_global[ibin][m]; m++; } } @@ -726,15 +739,12 @@ void ComputeStressMopProfile::compute_angles() // initialization for (int m = 0; m < nbins; m++) { - for (int i = 0; i < nvalues; i++) { - angle_local[m][i] = 0.0; - } + for (int i = 0; i < nvalues; i++) { angle_local[m][i] = 0.0; } local_contribution[m][0] = 0.0; local_contribution[m][1] = 0.0; local_contribution[m][2] = 0.0; } - for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; @@ -765,7 +775,7 @@ void ComputeStressMopProfile::compute_angles() if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; if (atype <= 0) continue; - for (int ibin = 0; ibin=0) && (tau_right <= 1)); - bool left_cross = ((tau_left >=0) && (tau_left <= 1)); + bool right_cross = ((tau_right >= 0) && (tau_right <= 1)); + bool left_cross = ((tau_left >= 0) && (tau_left <= 1)); // no bonds crossing the plane if (!right_cross && !left_cross) continue; // compute the cos(theta) of the angle - r1 = sqrt(dx_left[0]*dx_left[0] + dx_left[1]*dx_left[1] + dx_left[2]*dx_left[2]); - r2 = sqrt(dx_right[0]*dx_right[0] + dx_right[1]*dx_right[1] + dx_right[2]*dx_right[2]); - cos_theta = -(dx_right[0]*dx_left[0] + dx_right[1]*dx_left[1] + dx_right[2]*dx_left[2])/(r1*r2); + r1 = sqrt(dx_left[0] * dx_left[0] + dx_left[1] * dx_left[1] + dx_left[2] * dx_left[2]); + r2 = + sqrt(dx_right[0] * dx_right[0] + dx_right[1] * dx_right[1] + dx_right[2] * dx_right[2]); + cos_theta = + -(dx_right[0] * dx_left[0] + dx_right[1] * dx_left[1] + dx_right[2] * dx_left[2]) / + (r1 * r2); - if (cos_theta > 1.0) cos_theta = 1.0; + if (cos_theta > 1.0) cos_theta = 1.0; if (cos_theta < -1.0) cos_theta = -1.0; // The method returns derivative with regards to cos(theta) angle->born_matrix(atype, atom1, atom2, atom3, duang, du2ang); // only right bond crossing the plane - if (right_cross && !left_cross) - { + if (right_cross && !left_cross) { double sgn = copysign(1.0, x_angle_right[dir] - pos); - dcos_theta[0] = sgn*(dx_right[0]*cos_theta/r2 + dx_left[0]/r1)/r2; - dcos_theta[1] = sgn*(dx_right[1]*cos_theta/r2 + dx_left[1]/r1)/r2; - dcos_theta[2] = sgn*(dx_right[2]*cos_theta/r2 + dx_left[2]/r1)/r2; + dcos_theta[0] = sgn * (dx_right[0] * cos_theta / r2 + dx_left[0] / r1) / r2; + dcos_theta[1] = sgn * (dx_right[1] * cos_theta / r2 + dx_left[1] / r1) / r2; + dcos_theta[2] = sgn * (dx_right[2] * cos_theta / r2 + dx_left[2] / r1) / r2; } // only left bond crossing the plane - if (!right_cross && left_cross) - { + if (!right_cross && left_cross) { double sgn = copysign(1.0, x_angle_left[dir] - pos); - dcos_theta[0] = -sgn*(dx_left[0]*cos_theta/r1 + dx_right[0]/r2)/r1; - dcos_theta[1] = -sgn*(dx_left[1]*cos_theta/r1 + dx_right[1]/r2)/r1; - dcos_theta[2] = -sgn*(dx_left[2]*cos_theta/r1 + dx_right[2]/r2)/r1; + dcos_theta[0] = -sgn * (dx_left[0] * cos_theta / r1 + dx_right[0] / r2) / r1; + dcos_theta[1] = -sgn * (dx_left[1] * cos_theta / r1 + dx_right[1] / r2) / r1; + dcos_theta[2] = -sgn * (dx_left[2] * cos_theta / r1 + dx_right[2] / r2) / r1; } // both bonds crossing the plane - if (right_cross && left_cross) - { + if (right_cross && left_cross) { // due to right bond double sgn = copysign(1.0, x_angle_middle[dir] - pos); - dcos_theta[0] = -sgn*(dx_right[0]*cos_theta/r2 + dx_left[0]/r1)/r2; - dcos_theta[1] = -sgn*(dx_right[1]*cos_theta/r2 + dx_left[1]/r1)/r2; - dcos_theta[2] = -sgn*(dx_right[2]*cos_theta/r2 + dx_left[2]/r1)/r2; + dcos_theta[0] = -sgn * (dx_right[0] * cos_theta / r2 + dx_left[0] / r1) / r2; + dcos_theta[1] = -sgn * (dx_right[1] * cos_theta / r2 + dx_left[1] / r1) / r2; + dcos_theta[2] = -sgn * (dx_right[2] * cos_theta / r2 + dx_left[2] / r1) / r2; // due to left bond - dcos_theta[0] += sgn*(dx_left[0]*cos_theta/r1 + dx_right[0]/r2)/r1; - dcos_theta[1] += sgn*(dx_left[1]*cos_theta/r1 + dx_right[1]/r2)/r1; - dcos_theta[2] += sgn*(dx_left[2]*cos_theta/r1 + dx_right[2]/r2)/r1; + dcos_theta[0] += sgn * (dx_left[0] * cos_theta / r1 + dx_right[0] / r2) / r1; + dcos_theta[1] += sgn * (dx_left[1] * cos_theta / r1 + dx_right[1] / r2) / r1; + dcos_theta[2] += sgn * (dx_left[2] * cos_theta / r1 + dx_right[2] / r2) / r1; } // final contribution of the given angle term - local_contribution[ibin][0] += duang*dcos_theta[0]/area*nktv2p; - local_contribution[ibin][1] += duang*dcos_theta[1]/area*nktv2p; - local_contribution[ibin][2] += duang*dcos_theta[2]/area*nktv2p; + local_contribution[ibin][0] += duang * dcos_theta[0] / area * nktv2p; + local_contribution[ibin][1] += duang * dcos_theta[1] / area * nktv2p; + local_contribution[ibin][2] += duang * dcos_theta[2] / area * nktv2p; } } } @@ -863,8 +873,8 @@ void ComputeStressMopProfile::compute_angles() if (which[m] == CONF || which[m] == TOTAL || which[m] == ANGLE) { for (int ibin = 0; ibin < nbins; ibin++) { angle_local[ibin][m] = local_contribution[ibin][0]; - angle_local[ibin][m+1] = local_contribution[ibin][1]; - angle_local[ibin][m+2] = local_contribution[ibin][2]; + angle_local[ibin][m + 1] = local_contribution[ibin][1]; + angle_local[ibin][m + 2] = local_contribution[ibin][2]; } } m += 3; @@ -917,9 +927,7 @@ void ComputeStressMopProfile::compute_dihedrals() // initialization for (int m = 0; m < nbins; m++) { - for (int i = 0; i < nvalues; i++) { - dihedral_local[m][i] = 0.0; - } + for (int i = 0; i < nvalues; i++) { dihedral_local[m][i] = 0.0; } local_contribution[m][0] = 0.0; local_contribution[m][1] = 0.0; local_contribution[m][2] = 0.0; @@ -940,9 +948,9 @@ void ComputeStressMopProfile::compute_dihedrals() for (i = 0; i < nd; i++) { if (molecular == 1) { if (tag[atom2] != dihedral_atom2[atom2][i]) continue; - atom1 = atom->map(dihedral_atom1[atom2][i]); - atom3 = atom->map(dihedral_atom3[atom2][i]); - atom4 = atom->map(dihedral_atom4[atom2][i]); + atom1 = atom->map(dihedral_atom1[atom2][i]); + atom3 = atom->map(dihedral_atom3[atom2][i]); + atom4 = atom->map(dihedral_atom4[atom2][i]); } else { if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; tagprev = tag[atom2] - iatom - 1; @@ -955,7 +963,7 @@ void ComputeStressMopProfile::compute_dihedrals() if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; if (atom4 < 0 || !(mask[atom4] & groupbit)) continue; - for (int ibin = 0; ibin=0) && (tau_right <= 1)); + bool right_cross = ((tau_right >= 0) && (tau_right <= 1)); bool middle_cross = ((tau_middle >= 0) && (tau_middle <= 1)); - bool left_cross = ((tau_left >=0) && (tau_left <= 1)); + bool left_cross = ((tau_left >= 0) && (tau_left <= 1)); // no bonds crossing the plane if (!right_cross && !middle_cross && !left_cross) continue; @@ -1026,45 +1034,45 @@ void ComputeStressMopProfile::compute_dihedrals() vb3z = x_atom_4[2] - x_atom_3[2]; // c0 calculation - sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); - sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); - sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + sb1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z); + sb2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z); + sb3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z); rb1 = sqrt(sb1); rb3 = sqrt(sb3); - c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + c0 = (vb1x * vb3x + vb1y * vb3y + vb1z * vb3z) * rb1 * rb3; // 1st and 2nd angle - b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag2 = vb1x * vb1x + vb1y * vb1y + vb1z * vb1z; b1mag = sqrt(b1mag2); - b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag2 = vb2x * vb2x + vb2y * vb2y + vb2z * vb2z; b2mag = sqrt(b2mag2); - b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag2 = vb3x * vb3x + vb3y * vb3y + vb3z * vb3z; b3mag = sqrt(b3mag2); - ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; - r12c1 = 1.0 / (b1mag*b2mag); + ctmp = vb1x * vb2x + vb1y * vb2y + vb1z * vb2z; + r12c1 = 1.0 / (b1mag * b2mag); c1mag = ctmp * r12c1; - ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; - r12c2 = 1.0 / (b2mag*b3mag); + ctmp = vb2xm * vb3x + vb2ym * vb3y + vb2zm * vb3z; + r12c2 = 1.0 / (b2mag * b3mag); c2mag = ctmp * r12c2; // cos and sin of 2 angles and final c - sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sin2 = MAX(1.0 - c1mag * c1mag, 0.0); sc1 = sqrt(sin2); if (sc1 < SMALL) sc1 = SMALL; - sc1 = 1.0/sc1; + sc1 = 1.0 / sc1; - sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sin2 = MAX(1.0 - c2mag * c2mag, 0.0); sc2 = sqrt(sin2); if (sc2 < SMALL) sc2 = SMALL; - sc2 = 1.0/sc2; + sc2 = 1.0 / sc2; s1 = sc1 * sc1; s2 = sc2 * sc2; s12 = sc1 * sc2; - c = (c0 + c1mag*c2mag) * s12; + c = (c0 + c1mag * c2mag) * s12; // error check if (c > 1.0) c = 1.0; @@ -1074,36 +1082,35 @@ void ComputeStressMopProfile::compute_dihedrals() double a = dudih; c = c * a; s12 = s12 * a; - a11 = c*sb1*s1; - a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); - a33 = c*sb3*s2; - a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12); - a13 = -rb1*rb3*s12; - a23 = r12c2 * (c2mag*c*s2 + c1mag*s12); + a11 = c * sb1 * s1; + a22 = -sb2 * (2.0 * c0 * s12 - c * (s1 + s2)); + a33 = c * sb3 * s2; + a12 = -r12c1 * (c1mag * c * s1 + c2mag * s12); + a13 = -rb1 * rb3 * s12; + a23 = r12c2 * (c2mag * c * s2 + c1mag * s12); - sx2 = a12*vb1x + a22*vb2x + a23*vb3x; - sy2 = a12*vb1y + a22*vb2y + a23*vb3y; - sz2 = a12*vb1z + a22*vb2z + a23*vb3z; + sx2 = a12 * vb1x + a22 * vb2x + a23 * vb3x; + sy2 = a12 * vb1y + a22 * vb2y + a23 * vb3y; + sz2 = a12 * vb1z + a22 * vb2z + a23 * vb3z; - f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; - f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; - f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; + f1[0] = a11 * vb1x + a12 * vb2x + a13 * vb3x; + f1[1] = a11 * vb1y + a12 * vb2y + a13 * vb3y; + f1[2] = a11 * vb1z + a12 * vb2z + a13 * vb3z; f2[0] = -sx2 - f1[0]; f2[1] = -sy2 - f1[1]; f2[2] = -sz2 - f1[2]; - f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; - f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; - f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; + f4[0] = a13 * vb1x + a23 * vb2x + a33 * vb3x; + f4[1] = a13 * vb1y + a23 * vb2y + a33 * vb3y; + f4[2] = a13 * vb1z + a23 * vb2z + a33 * vb3z; f3[0] = sx2 - f4[0]; f3[1] = sy2 - f4[1]; f3[2] = sz2 - f4[2]; // only right bond crossing the plane - if (right_cross && !middle_cross && !left_cross) - { + if (right_cross && !middle_cross && !left_cross) { double sgn = copysign(1.0, x_atom_1[dir] - pos); df[0] = sgn * f1[0]; df[1] = sgn * f1[1]; @@ -1111,8 +1118,7 @@ void ComputeStressMopProfile::compute_dihedrals() } // only middle bond crossing the plane - if (!right_cross && middle_cross && !left_cross) - { + if (!right_cross && middle_cross && !left_cross) { double sgn = copysign(1.0, x_atom_2[dir] - pos); df[0] = sgn * (f2[0] + f1[0]); df[1] = sgn * (f2[1] + f1[1]); @@ -1120,8 +1126,7 @@ void ComputeStressMopProfile::compute_dihedrals() } // only left bond crossing the plane - if (!right_cross && !middle_cross && left_cross) - { + if (!right_cross && !middle_cross && left_cross) { double sgn = copysign(1.0, x_atom_4[dir] - pos); df[0] = sgn * f4[0]; df[1] = sgn * f4[1]; @@ -1129,8 +1134,7 @@ void ComputeStressMopProfile::compute_dihedrals() } // only right & middle bonds crossing the plane - if (right_cross && middle_cross && !left_cross) - { + if (right_cross && middle_cross && !left_cross) { double sgn = copysign(1.0, x_atom_2[dir] - pos); df[0] = sgn * f2[0]; df[1] = sgn * f2[1]; @@ -1138,8 +1142,7 @@ void ComputeStressMopProfile::compute_dihedrals() } // only right & left bonds crossing the plane - if (right_cross && !middle_cross && left_cross) - { + if (right_cross && !middle_cross && left_cross) { double sgn = copysign(1.0, x_atom_1[dir] - pos); df[0] = sgn * (f1[0] + f4[0]); df[1] = sgn * (f1[1] + f4[1]); @@ -1147,8 +1150,7 @@ void ComputeStressMopProfile::compute_dihedrals() } // only middle & left bonds crossing the plane - if (!right_cross && middle_cross && left_cross) - { + if (!right_cross && middle_cross && left_cross) { double sgn = copysign(1.0, x_atom_3[dir] - pos); df[0] = sgn * f3[0]; df[1] = sgn * f3[1]; @@ -1156,17 +1158,16 @@ void ComputeStressMopProfile::compute_dihedrals() } // all three bonds crossing the plane - if (right_cross && middle_cross && left_cross) - { + if (right_cross && middle_cross && left_cross) { double sgn = copysign(1.0, x_atom_1[dir] - pos); df[0] = sgn * (f1[0] + f3[0]); df[1] = sgn * (f1[1] + f3[1]); df[2] = sgn * (f1[2] + f3[2]); } - local_contribution[ibin][0] += df[0]/area*nktv2p; - local_contribution[ibin][1] += df[1]/area*nktv2p; - local_contribution[ibin][2] += df[2]/area*nktv2p; + local_contribution[ibin][0] += df[0] / area * nktv2p; + local_contribution[ibin][1] += df[1] / area * nktv2p; + local_contribution[ibin][2] += df[2] / area * nktv2p; } } } @@ -1177,13 +1178,12 @@ void ComputeStressMopProfile::compute_dihedrals() if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == DIHEDRAL)) { for (int ibin = 0; ibin < nbins; ibin++) { dihedral_local[ibin][m] = local_contribution[ibin][0]; - dihedral_local[ibin][m+1] = local_contribution[ibin][1]; - dihedral_local[ibin][m+2] = local_contribution[ibin][2]; + dihedral_local[ibin][m + 1] = local_contribution[ibin][1]; + dihedral_local[ibin][m + 2] = local_contribution[ibin][2]; } } m += 3; } - } /* ---------------------------------------------------------------------- @@ -1203,11 +1203,11 @@ void ComputeStressMopProfile::setup_bins() if ((origin > domain->boxhi[dir]) || (origin < domain->boxlo[dir])) error->all(FLERR, "Origin of bins for compute stress/mop/profile is out of bounds"); - n = static_cast ((origin - boxlo[dir]) * invdelta); - lo = origin - n*delta; + n = static_cast((origin - boxlo[dir]) * invdelta); + lo = origin - n * delta; - n = static_cast ((boxhi[dir] - origin) * invdelta); - hi = origin + n*delta; + n = static_cast((boxhi[dir] - origin) * invdelta); + hi = origin + n * delta; offset = lo; nbins = static_cast((hi - lo) * invdelta + 1.5); @@ -1221,8 +1221,8 @@ void ComputeStressMopProfile::setup_bins() memory->create(bond_global, nbins, nvalues, "stress/mop/profile:bond_global"); memory->create(angle_local, nbins, nvalues, "stress/mop/profile:angle_local"); memory->create(angle_global, nbins, nvalues, "stress/mop/profile:angle_global"); - memory->create(dihedral_local,nbins,nvalues,"stress/mop/profile:dihedral_local"); - memory->create(dihedral_global,nbins,nvalues,"stress/mop/profile:dihedral_global"); + memory->create(dihedral_local, nbins, nvalues, "stress/mop/profile:dihedral_local"); + memory->create(dihedral_global, nbins, nvalues, "stress/mop/profile:dihedral_global"); memory->create(local_contribution, nbins, 3, "stress/mop/profile:local_contribution"); // set bin coordinates diff --git a/src/MC/fix_atom_swap.cpp b/src/MC/fix_atom_swap.cpp index 1348f1248e..209c5cf8c4 100644 --- a/src/MC/fix_atom_swap.cpp +++ b/src/MC/fix_atom_swap.cpp @@ -35,6 +35,7 @@ #include "modify.h" #include "neighbor.h" #include "pair.h" +#include "pair_hybrid.h" #include "random_park.h" #include "region.h" #include "update.h" @@ -223,6 +224,28 @@ void FixAtomSwap::init() error->all(FLERR, "Mu not allowed when not using semi-grand in fix atom/swap command"); } + // check if constraints for hybrid pair styles are fulfilled + + if (utils::strmatch(force->pair_style, "^hybrid")) { + auto *hybrid = dynamic_cast(force->pair); + if (hybrid) { + for (int i = 0; i < nswaptypes - 1; ++i) { + int type1 = type_list[i]; + for (int j = i + 1; j < nswaptypes; ++j) { + int type2 = type_list[j]; + if (hybrid->nmap[type1][type1] != hybrid->nmap[type2][type2]) + error->all(FLERR, "Pair {} substyles for atom types {} and {} are not compatible", + force->pair_style, type1, type2); + for (int k = 0; k < hybrid->nmap[type1][type1]; ++k) { + if (hybrid->map[type1][type1][k] != hybrid->map[type2][type2][k]) + error->all(FLERR, "Pair {} substyles for atom types {} and {} are not compatible", + force->pair_style, type1, type2); + } + } + } + } + } + // set index and check validity of region if (idregion) { diff --git a/src/compute_angle_local.cpp b/src/compute_angle_local.cpp index 3e8b15fd64..b4a9334416 100644 --- a/src/compute_angle_local.cpp +++ b/src/compute_angle_local.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -13,38 +12,40 @@ ------------------------------------------------------------------------- */ #include "compute_angle_local.h" -#include -#include + +#include "angle.h" #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "update.h" #include "domain.h" +#include "error.h" #include "force.h" -#include "angle.h" #include "input.h" -#include "variable.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "molecule.h" +#include "update.h" +#include "variable.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; static constexpr int DELTA = 10000; -enum{THETA,ENG,VARIABLE}; +enum { THETA, ENG, VARIABLE }; /* ---------------------------------------------------------------------- */ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - bstyle(nullptr), vvar(nullptr), tstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr) + Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), tstr(nullptr), vstr(nullptr), + vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute angle/local command"); + if (narg < 4) error->all(FLERR, "Illegal compute angle/local command"); if (atom->avec->angles_allow == 0) - error->all(FLERR,"Compute angle/local used when angles are not allowed"); + error->all(FLERR, "Compute angle/local used when angles are not allowed"); local_flag = 1; @@ -52,7 +53,7 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : nvalues = narg - 3; bstyle = new int[nvalues]; - vstr = new char*[nvalues]; + vstr = new char *[nvalues]; vvar = new int[nvalues]; nvalues = 0; @@ -61,16 +62,17 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : int iarg; for (iarg = 3; iarg < narg; iarg++) { - if (strcmp(arg[iarg],"theta") == 0) { + if (strcmp(arg[iarg], "theta") == 0) { bstyle[nvalues++] = THETA; tflag = 1; - } else if (strcmp(arg[iarg],"eng") == 0) { + } else if (strcmp(arg[iarg], "eng") == 0) { bstyle[nvalues++] = ENG; - } else if (strncmp(arg[iarg],"v_",2) == 0) { + } else if (strncmp(arg[iarg], "v_", 2) == 0) { bstyle[nvalues++] = VARIABLE; vstr[nvar] = utils::strdup(&arg[iarg][2]); nvar++; - } else break; + } else + break; } // optional args @@ -79,45 +81,46 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : tstr = nullptr; while (iarg < narg) { - if (strcmp(arg[iarg],"set") == 0) { + if (strcmp(arg[iarg], "set") == 0) { setflag = 1; - if (iarg+3 > narg) error->all(FLERR,"Illegal compute angle/local command"); - if (strcmp(arg[iarg+1],"theta") == 0) { - delete [] tstr; - tstr = utils::strdup(arg[iarg+2]); + if (iarg + 3 > narg) error->all(FLERR, "Illegal compute angle/local command"); + if (strcmp(arg[iarg + 1], "theta") == 0) { + delete[] tstr; + tstr = utils::strdup(arg[iarg + 2]); tflag = 1; - } else error->all(FLERR,"Illegal compute angle/local command"); + } else + error->all(FLERR, "Illegal compute angle/local command"); iarg += 3; - } else error->all(FLERR,"Illegal compute angle/local command"); + } else + error->all(FLERR, "Illegal compute angle/local command"); } // error check if (nvar) { - if (!setflag) - error->all(FLERR,"Compute angle/local variable requires a set variable"); + if (!setflag) error->all(FLERR, "Compute angle/local variable requires a set variable"); for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) - error->all(FLERR,"Variable name for copute angle/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for compute angle/local does not exist"); if (!input->variable->equalstyle(vvar[i])) - error->all(FLERR,"Variable for compute angle/local is invalid style"); + error->all(FLERR, "Variable for compute angle/local is invalid style"); } if (tstr) { tvar = input->variable->find(tstr); - if (tvar < 0) - error->all(FLERR,"Variable name for compute angle/local does not exist"); + if (tvar < 0) error->all(FLERR, "Variable name for compute angle/local does not exist"); if (!input->variable->internalstyle(tvar)) - error->all(FLERR,"Variable for compute angle/local is invalid style"); + error->all(FLERR, "Variable for compute angle/local is invalid style"); } } else if (setflag) - error->all(FLERR,"Compute angle/local set with no variable"); + error->all(FLERR, "Compute angle/local set with no variable"); // initialize output - if (nvalues == 1) size_local_cols = 0; - else size_local_cols = nvalues; + if (nvalues == 1) + size_local_cols = 0; + else + size_local_cols = nvalues; nmax = 0; vlocal = nullptr; @@ -128,12 +131,12 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : ComputeAngleLocal::~ComputeAngleLocal() { - delete [] bstyle; - for (int i = 0; i < nvar; i++) delete [] vstr[i]; - delete [] vstr; - delete [] vvar; + delete[] bstyle; + for (int i = 0; i < nvar; i++) delete[] vstr[i]; + delete[] vstr; + delete[] vvar; - delete [] tstr; + delete[] tstr; memory->destroy(vlocal); memory->destroy(alocal); @@ -144,19 +147,17 @@ ComputeAngleLocal::~ComputeAngleLocal() void ComputeAngleLocal::init() { if (force->angle == nullptr) - error->all(FLERR,"No angle style is defined for compute angle/local"); + error->all(FLERR, "No angle style is defined for compute angle/local"); if (nvar) { for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) - error->all(FLERR,"Variable name for compute angle/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for compute angle/local does not exist"); } if (tstr) { tvar = input->variable->find(tstr); - if (tvar < 0) - error->all(FLERR,"Variable name for compute angle/local does not exist"); + if (tvar < 0) error->all(FLERR, "Variable name for compute angle/local does not exist"); } } @@ -194,10 +195,10 @@ void ComputeAngleLocal::compute_local() int ComputeAngleLocal::compute_angles(int flag) { - int i,m,na,atom1,atom2,atom3,imol,iatom,atype,ivar; + int i, m, na, atom1, atom2, atom3, imol, iatom, atype, ivar; tagint tagprev; - double delx1,dely1,delz1,delx2,dely2,delz2; - double rsq1,rsq2,r1,r2,c,theta; + double delx1, dely1, delz1, delx2, dely2, delz2; + double rsq1, rsq2, r1, r2, c, theta; double *ptr; double **x = atom->x; @@ -224,7 +225,8 @@ int ComputeAngleLocal::compute_angles(int flag) for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; - if (molecular == Atom::MOLECULAR) na = num_angle[atom2]; + if (molecular == Atom::MOLECULAR) + na = num_angle[atom2]; else { if (molindex[atom2] < 0) continue; imol = molindex[atom2]; @@ -242,8 +244,8 @@ int ComputeAngleLocal::compute_angles(int flag) if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue; atype = onemols[imol]->angle_type[atom2][i]; tagprev = tag[atom2] - iatom - 1; - atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i]+tagprev); - atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i]+tagprev); + atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i] + tagprev); } if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; @@ -261,50 +263,54 @@ int ComputeAngleLocal::compute_angles(int flag) delx1 = x[atom1][0] - x[atom2][0]; dely1 = x[atom1][1] - x[atom2][1]; delz1 = x[atom1][2] - x[atom2][2]; - domain->minimum_image(delx1,dely1,delz1); + domain->minimum_image(delx1, dely1, delz1); - rsq1 = delx1*delx1 + dely1*dely1 + delz1*delz1; + rsq1 = delx1 * delx1 + dely1 * dely1 + delz1 * delz1; r1 = sqrt(rsq1); delx2 = x[atom3][0] - x[atom2][0]; dely2 = x[atom3][1] - x[atom2][1]; delz2 = x[atom3][2] - x[atom2][2]; - domain->minimum_image(delx2,dely2,delz2); + domain->minimum_image(delx2, dely2, delz2); - rsq2 = delx2*delx2 + dely2*dely2 + delz2*delz2; + rsq2 = delx2 * delx2 + dely2 * dely2 + delz2 * delz2; r2 = sqrt(rsq2); // c = cosine of angle // theta = angle in radians - c = delx1*delx2 + dely1*dely2 + delz1*delz2; - c /= r1*r2; + c = delx1 * delx2 + dely1 * dely2 + delz1 * delz2; + c /= r1 * r2; if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; theta = acos(c); } - if (nvalues == 1) ptr = &vlocal[m]; - else ptr = alocal[m]; + if (nvalues == 1) + ptr = &vlocal[m]; + else + ptr = alocal[m]; if (nvar) { ivar = 0; - if (tstr) input->variable->internal_set(tvar,theta); + if (tstr) input->variable->internal_set(tvar, theta); } for (int n = 0; n < nvalues; n++) { switch (bstyle[n]) { - case THETA: - ptr[n] = 180.0*theta/MY_PI; - break; - case ENG: - if (atype > 0) ptr[n] = angle->single(atype,atom1,atom2,atom3); - else ptr[n] = 0.0; - break; - case VARIABLE: - ptr[n] = input->variable->compute_equal(vvar[ivar]); - ivar++; - break; + case THETA: + ptr[n] = 180.0 * theta / MY_PI; + break; + case ENG: + if (atype > 0) + ptr[n] = angle->single(atype, atom1, atom2, atom3); + else + ptr[n] = 0.0; + break; + case VARIABLE: + ptr[n] = input->variable->compute_equal(vvar[ivar]); + ivar++; + break; } } @@ -325,11 +331,11 @@ void ComputeAngleLocal::reallocate(int n) if (nvalues == 1) { memory->destroy(vlocal); - memory->create(vlocal,nmax,"angle/local:vector_local"); + memory->create(vlocal, nmax, "angle/local:vector_local"); vector_local = vlocal; } else { memory->destroy(alocal); - memory->create(alocal,nmax,nvalues,"angle/local:array_local"); + memory->create(alocal, nmax, nvalues, "angle/local:array_local"); array_local = alocal; } } @@ -340,6 +346,6 @@ void ComputeAngleLocal::reallocate(int n) double ComputeAngleLocal::memory_usage() { - double bytes = (double)nmax*nvalues * sizeof(double); + double bytes = (double) nmax * nvalues * sizeof(double); return bytes; } diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index e3c3e26a85..e9632d254f 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -35,18 +34,35 @@ using namespace LAMMPS_NS; static constexpr int DELTA = 10000; -enum{DIST,DX,DY,DZ,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,FX,FY,FZ,VARIABLE,BN}; +enum { + DIST, + DX, + DY, + DZ, + VELVIB, + OMEGA, + ENGTRANS, + ENGVIB, + ENGROT, + ENGPOT, + FORCE, + FX, + FY, + FZ, + VARIABLE, + BN +}; /* ---------------------------------------------------------------------- */ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - bstyle(nullptr), vvar(nullptr), dstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr) + Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), dstr(nullptr), vstr(nullptr), + vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute bond/local command"); + if (narg < 4) error->all(FLERR, "Illegal compute bond/local command"); if (atom->avec->bonds_allow == 0) - error->all(FLERR,"Compute bond/local used when bonds are not allowed"); + error->all(FLERR, "Compute bond/local used when bonds are not allowed"); local_flag = 1; comm_forward = 3; @@ -56,7 +72,7 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : nvalues = narg - 3; bstyle = new int[nvalues]; bindex = new int[nvalues]; - vstr = new char*[nvalues]; + vstr = new char *[nvalues]; vvar = new int[nvalues]; nvalues = 0; @@ -64,30 +80,45 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : int iarg; for (iarg = 3; iarg < narg; iarg++) { - if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST; - else if (strcmp(arg[iarg],"dx") == 0) bstyle[nvalues++] = DX; - else if (strcmp(arg[iarg],"dy") == 0) bstyle[nvalues++] = DY; - else if (strcmp(arg[iarg],"dz") == 0) bstyle[nvalues++] = DZ; - else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT; - else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE; - else if (strcmp(arg[iarg],"fx") == 0) bstyle[nvalues++] = FX; - else if (strcmp(arg[iarg],"fy") == 0) bstyle[nvalues++] = FY; - else if (strcmp(arg[iarg],"fz") == 0) bstyle[nvalues++] = FZ; - else if (strcmp(arg[iarg],"engvib") == 0) bstyle[nvalues++] = ENGVIB; - else if (strcmp(arg[iarg],"engrot") == 0) bstyle[nvalues++] = ENGROT; - else if (strcmp(arg[iarg],"engtrans") == 0) bstyle[nvalues++] = ENGTRANS; - else if (strcmp(arg[iarg],"omega") == 0) bstyle[nvalues++] = OMEGA; - else if (strcmp(arg[iarg],"velvib") == 0) bstyle[nvalues++] = VELVIB; - else if (strncmp(arg[iarg],"v_",2) == 0) { + if (strcmp(arg[iarg], "dist") == 0) + bstyle[nvalues++] = DIST; + else if (strcmp(arg[iarg], "dx") == 0) + bstyle[nvalues++] = DX; + else if (strcmp(arg[iarg], "dy") == 0) + bstyle[nvalues++] = DY; + else if (strcmp(arg[iarg], "dz") == 0) + bstyle[nvalues++] = DZ; + else if (strcmp(arg[iarg], "engpot") == 0) + bstyle[nvalues++] = ENGPOT; + else if (strcmp(arg[iarg], "force") == 0) + bstyle[nvalues++] = FORCE; + else if (strcmp(arg[iarg], "fx") == 0) + bstyle[nvalues++] = FX; + else if (strcmp(arg[iarg], "fy") == 0) + bstyle[nvalues++] = FY; + else if (strcmp(arg[iarg], "fz") == 0) + bstyle[nvalues++] = FZ; + else if (strcmp(arg[iarg], "engvib") == 0) + bstyle[nvalues++] = ENGVIB; + else if (strcmp(arg[iarg], "engrot") == 0) + bstyle[nvalues++] = ENGROT; + else if (strcmp(arg[iarg], "engtrans") == 0) + bstyle[nvalues++] = ENGTRANS; + else if (strcmp(arg[iarg], "omega") == 0) + bstyle[nvalues++] = OMEGA; + else if (strcmp(arg[iarg], "velvib") == 0) + bstyle[nvalues++] = VELVIB; + else if (strncmp(arg[iarg], "v_", 2) == 0) { bstyle[nvalues++] = VARIABLE; vstr[nvar] = utils::strdup(&arg[iarg][2]); nvar++; - } else if (utils::strmatch(arg[iarg], "^b\\d+$")) { // b1, b2, b3, ... bN + } else if (utils::strmatch(arg[iarg], "^b\\d+$")) { // b1, b2, b3, ... bN int n = std::stoi(&arg[iarg][1]); if (n <= 0) error->all(FLERR, "Invalid keyword {} in compute bond/local command", arg[iarg]); bstyle[nvalues] = BN; bindex[nvalues++] = n - 1; - } else break; + } else + break; } // optional args @@ -96,40 +127,39 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : dstr = nullptr; while (iarg < narg) { - if (strcmp(arg[iarg],"set") == 0) { + if (strcmp(arg[iarg], "set") == 0) { setflag = 1; - if (iarg+3 > narg) utils::missing_cmd_args(FLERR,"compute bond/local set", error); - if (strcmp(arg[iarg+1],"dist") == 0) { - delete [] dstr; - dstr = utils::strdup(arg[iarg+2]); - } else error->all(FLERR,"Unknown compute bond/local set keyword: {}", arg[iarg+2]); + if (iarg + 3 > narg) utils::missing_cmd_args(FLERR, "compute bond/local set", error); + if (strcmp(arg[iarg + 1], "dist") == 0) { + delete[] dstr; + dstr = utils::strdup(arg[iarg + 2]); + } else + error->all(FLERR, "Unknown compute bond/local set keyword: {}", arg[iarg + 2]); iarg += 3; - } else error->all(FLERR,"Unknown compute bond/local keyword: {}", arg[iarg]); + } else + error->all(FLERR, "Unknown compute bond/local keyword: {}", arg[iarg]); } // error check if (nvar) { - if (!setflag) - error->all(FLERR,"Compute bond/local variable requires a set variable"); + if (!setflag) error->all(FLERR, "Compute bond/local variable requires a set variable"); for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); if (vvar[i] < 0) - error->all(FLERR,"Variable name {} for copute bond/local does not exist", vstr[i]); + error->all(FLERR, "Variable name {} for compute bond/local does not exist", vstr[i]); if (!input->variable->equalstyle(vvar[i])) - error->all(FLERR,"Variable {} for compute bond/local is invalid style", vstr[i]); + error->all(FLERR, "Variable {} for compute bond/local is invalid style", vstr[i]); } if (dstr) { dvar = input->variable->find(dstr); - if (dvar < 0) - error->all(FLERR,"Variable name for compute bond/local does not exist"); + if (dvar < 0) error->all(FLERR, "Variable name for compute bond/local does not exist"); if (!input->variable->internalstyle(dvar)) - error->all(FLERR,"Variable for compute bond/local is invalid style"); + error->all(FLERR, "Variable for compute bond/local is invalid style"); } } else if (setflag) - error->all(FLERR,"Compute bond/local set used with without a variable"); - + error->all(FLERR, "Compute bond/local set used with without a variable"); // set singleflag if need to call bond->single() // set velflag if compute any quantities based on velocities @@ -137,16 +167,20 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : singleflag = 0; velflag = 0; for (int i = 0; i < nvalues; i++) { - if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX || - bstyle[i] == FY || bstyle[i] == FZ || bstyle[i] == BN) singleflag = 1; - if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS || - bstyle[i] == ENGVIB || bstyle[i] == ENGROT) velflag = 1; + if (bstyle[i] == ENGPOT || bstyle[i] == FORCE || bstyle[i] == FX || bstyle[i] == FY || + bstyle[i] == FZ || bstyle[i] == BN) + singleflag = 1; + if (bstyle[i] == VELVIB || bstyle[i] == OMEGA || bstyle[i] == ENGTRANS || bstyle[i] == ENGVIB || + bstyle[i] == ENGROT) + velflag = 1; } // initialize output - if (nvalues == 1) size_local_cols = 0; - else size_local_cols = nvalues; + if (nvalues == 1) + size_local_cols = 0; + else + size_local_cols = nvalues; nmax = 0; vlocal = nullptr; @@ -157,13 +191,13 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : ComputeBondLocal::~ComputeBondLocal() { - delete [] bstyle; - delete [] bindex; - for (int i = 0; i < nvar; i++) delete [] vstr[i]; - delete [] vstr; - delete [] vvar; + delete[] bstyle; + delete[] bindex; + for (int i = 0; i < nvar; i++) delete[] vstr[i]; + delete[] vstr; + delete[] vvar; - delete [] dstr; + delete[] dstr; memory->destroy(vlocal); memory->destroy(alocal); @@ -173,8 +207,7 @@ ComputeBondLocal::~ComputeBondLocal() void ComputeBondLocal::init() { - if (force->bond == nullptr) - error->all(FLERR,"No bond style is defined for compute bond/local"); + if (force->bond == nullptr) error->all(FLERR, "No bond style is defined for compute bond/local"); for (int i = 0; i < nvalues; i++) if (bstyle[i] == BN && bindex[i] >= force->bond->single_extra) @@ -183,21 +216,21 @@ void ComputeBondLocal::init() if (nvar) { for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) - error->all(FLERR,"Variable name for compute bond/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for compute bond/local does not exist"); } if (dstr) { dvar = input->variable->find(dstr); - if (dvar < 0) - error->all(FLERR,"Variable name for compute bond/local does not exist"); + if (dvar < 0) error->all(FLERR, "Variable name for compute bond/local does not exist"); } } // set ghostvelflag if need to acquire ghost atom velocities - if (velflag && !comm->ghost_velocity) ghostvelflag = 1; - else ghostvelflag = 0; + if (velflag && !comm->ghost_velocity) + ghostvelflag = 1; + else + ghostvelflag = 0; // do initial memory allocation so that memory_usage() is correct @@ -236,17 +269,17 @@ void ComputeBondLocal::compute_local() int ComputeBondLocal::compute_bonds(int flag) { - int i,m,nb,atom1,atom2,imol,iatom,btype,ivar; + int i, m, nb, atom1, atom2, imol, iatom, btype, ivar; tagint tagprev; - double dx,dy,dz,rsq; - double mass1,mass2,masstotal,invmasstotal; - double xcm[3],vcm[3]; - double delr1[3],delr2[3],delv1[3],delv2[3]; - double r12[3],vpar1,vpar2; - double vvib,vrotsq; - double inertia,omegasq; + double dx, dy, dz, rsq; + double mass1, mass2, masstotal, invmasstotal; + double xcm[3], vcm[3]; + double delr1[3], delr2[3], delv1[3], delv2[3]; + double r12[3], vpar1, vpar2; + double vvib, vrotsq; + double inertia, omegasq; double mvv2e; - double engpot,engtrans,engvib,engrot,fbond; + double engpot, engtrans, engvib, engrot, fbond; double *ptr; double **x = atom->x; @@ -280,7 +313,8 @@ int ComputeBondLocal::compute_bonds(int flag) for (atom1 = 0; atom1 < nlocal; atom1++) { if (!(mask[atom1] & groupbit)) continue; - if (molecular == Atom::MOLECULAR) nb = num_bond[atom1]; + if (molecular == Atom::MOLECULAR) + nb = num_bond[atom1]; else { if (molindex[atom1] < 0) continue; imol = molindex[atom1]; @@ -295,7 +329,7 @@ int ComputeBondLocal::compute_bonds(int flag) } else { tagprev = tag[atom1] - iatom - 1; btype = onemols[imol]->bond_type[iatom][i]; - atom2 = atom->map(onemols[imol]->bond_atom[iatom][i]+tagprev); + atom2 = atom->map(onemols[imol]->bond_atom[iatom][i] + tagprev); } if (atom2 < 0 || !(mask[atom2] & groupbit)) continue; @@ -310,55 +344,56 @@ int ComputeBondLocal::compute_bonds(int flag) dx = x[atom1][0] - x[atom2][0]; dy = x[atom1][1] - x[atom2][1]; dz = x[atom1][2] - x[atom2][2]; - domain->minimum_image(dx,dy,dz); - rsq = dx*dx + dy*dy + dz*dz; + domain->minimum_image(dx, dy, dz); + rsq = dx * dx + dy * dy + dz * dz; if (btype == 0) { fbond = 0.0; } else { - if (singleflag) engpot = bond->single(btype,rsq,atom1,atom2,fbond); - else fbond = engpot = 0.0; + if (singleflag) + engpot = bond->single(btype, rsq, atom1, atom2, fbond); + else + fbond = engpot = 0.0; engvib = engrot = engtrans = omegasq = vvib = 0.0; if (velflag) { if (rmass) { mass1 = rmass[atom1]; mass2 = rmass[atom2]; - } - else { + } else { mass1 = mass[type[atom1]]; mass2 = mass[type[atom2]]; } - masstotal = mass1+mass2; + masstotal = mass1 + mass2; invmasstotal = 1.0 / (masstotal); - xcm[0] = (mass1*x[atom1][0] + mass2*x[atom2][0]) * invmasstotal; - xcm[1] = (mass1*x[atom1][1] + mass2*x[atom2][1]) * invmasstotal; - xcm[2] = (mass1*x[atom1][2] + mass2*x[atom2][2]) * invmasstotal; - vcm[0] = (mass1*v[atom1][0] + mass2*v[atom2][0]) * invmasstotal; - vcm[1] = (mass1*v[atom1][1] + mass2*v[atom2][1]) * invmasstotal; - vcm[2] = (mass1*v[atom1][2] + mass2*v[atom2][2]) * invmasstotal; + xcm[0] = (mass1 * x[atom1][0] + mass2 * x[atom2][0]) * invmasstotal; + xcm[1] = (mass1 * x[atom1][1] + mass2 * x[atom2][1]) * invmasstotal; + xcm[2] = (mass1 * x[atom1][2] + mass2 * x[atom2][2]) * invmasstotal; + vcm[0] = (mass1 * v[atom1][0] + mass2 * v[atom2][0]) * invmasstotal; + vcm[1] = (mass1 * v[atom1][1] + mass2 * v[atom2][1]) * invmasstotal; + vcm[2] = (mass1 * v[atom1][2] + mass2 * v[atom2][2]) * invmasstotal; engtrans = 0.5 * masstotal * MathExtra::lensq3(vcm); // r12 = unit bond vector from atom1 to atom2 - MathExtra::sub3(x[atom2],x[atom1],r12); + MathExtra::sub3(x[atom2], x[atom1], r12); MathExtra::norm3(r12); // delr = vector from COM to each atom // delv = velocity of each atom relative to COM - MathExtra::sub3(x[atom1],xcm,delr1); - MathExtra::sub3(x[atom2],xcm,delr2); - MathExtra::sub3(v[atom1],vcm,delv1); - MathExtra::sub3(v[atom2],vcm,delv2); + MathExtra::sub3(x[atom1], xcm, delr1); + MathExtra::sub3(x[atom2], xcm, delr2); + MathExtra::sub3(v[atom1], vcm, delv1); + MathExtra::sub3(v[atom2], vcm, delv2); // vpar = component of delv parallel to bond vector - vpar1 = MathExtra::dot3(delv1,r12); - vpar2 = MathExtra::dot3(delv2,r12); - engvib = 0.5 * (mass1*vpar1*vpar1 + mass2*vpar2*vpar2); + vpar1 = MathExtra::dot3(delv1, r12); + vpar2 = MathExtra::dot3(delv2, r12); + engvib = 0.5 * (mass1 * vpar1 * vpar1 + mass2 * vpar2 * vpar2); // vvib = relative velocity of 2 atoms along bond direction // vvib < 0 for 2 atoms moving towards each other @@ -369,9 +404,8 @@ int ComputeBondLocal::compute_bonds(int flag) // vrotsq = tangential speed squared of atom1 only // omegasq = omega squared, and is the same for atom1 and atom2 - inertia = mass1*MathExtra::lensq3(delr1) + - mass2*MathExtra::lensq3(delr2); - vrotsq = MathExtra::lensq3(delv1) - vpar1*vpar1; + inertia = mass1 * MathExtra::lensq3(delr1) + mass2 * MathExtra::lensq3(delr2); + vrotsq = MathExtra::lensq3(delv1) - vpar1 * vpar1; omegasq = vrotsq / MathExtra::lensq3(delr1); engrot = 0.5 * inertia * omegasq; @@ -384,12 +418,14 @@ int ComputeBondLocal::compute_bonds(int flag) engrot *= mvv2e; } - if (nvalues == 1) ptr = &vlocal[m]; - else ptr = alocal[m]; + if (nvalues == 1) + ptr = &vlocal[m]; + else + ptr = alocal[m]; if (nvar) { ivar = 0; - if (dstr) input->variable->internal_set(dvar,sqrt(rsq)); + if (dstr) input->variable->internal_set(dvar, sqrt(rsq)); } // to make sure dx, dy and dz are always from the lower to the higher id @@ -397,55 +433,55 @@ int ComputeBondLocal::compute_bonds(int flag) for (int n = 0; n < nvalues; n++) { switch (bstyle[n]) { - case DIST: - ptr[n] = sqrt(rsq); - break; - case DX: - ptr[n] = dx*directionCorrection; - break; - case DY: - ptr[n] = dy*directionCorrection; - break; - case DZ: - ptr[n] = dz*directionCorrection; - break; - case ENGPOT: - ptr[n] = engpot; - break; - case FORCE: - ptr[n] = sqrt(rsq)*fbond; - break; - case FX: - ptr[n] = dx*fbond; - break; - case FY: - ptr[n] = dy*fbond; - break; - case FZ: - ptr[n] = dz*fbond; - break; - case ENGVIB: - ptr[n] = engvib; - break; - case ENGROT: - ptr[n] = engrot; - break; - case ENGTRANS: - ptr[n] = engtrans; - break; - case OMEGA: - ptr[n] = sqrt(omegasq); - break; - case VELVIB: - ptr[n] = vvib; - break; - case VARIABLE: - ptr[n] = input->variable->compute_equal(vvar[ivar]); - ivar++; - break; - case BN: - ptr[n] = bond->svector[bindex[n]]; - break; + case DIST: + ptr[n] = sqrt(rsq); + break; + case DX: + ptr[n] = dx * directionCorrection; + break; + case DY: + ptr[n] = dy * directionCorrection; + break; + case DZ: + ptr[n] = dz * directionCorrection; + break; + case ENGPOT: + ptr[n] = engpot; + break; + case FORCE: + ptr[n] = sqrt(rsq) * fbond; + break; + case FX: + ptr[n] = dx * fbond; + break; + case FY: + ptr[n] = dy * fbond; + break; + case FZ: + ptr[n] = dz * fbond; + break; + case ENGVIB: + ptr[n] = engvib; + break; + case ENGROT: + ptr[n] = engrot; + break; + case ENGTRANS: + ptr[n] = engtrans; + break; + case OMEGA: + ptr[n] = sqrt(omegasq); + break; + case VELVIB: + ptr[n] = vvib; + break; + case VARIABLE: + ptr[n] = input->variable->compute_equal(vvar[ivar]); + ivar++; + break; + case BN: + ptr[n] = bond->svector[bindex[n]]; + break; } } } @@ -459,10 +495,10 @@ int ComputeBondLocal::compute_bonds(int flag) /* ---------------------------------------------------------------------- */ -int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) +int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf, int /*pbc_flag*/, + int * /*pbc*/) { - int i,j,m; + int i, j, m; double **v = atom->v; @@ -481,7 +517,7 @@ int ComputeBondLocal::pack_forward_comm(int n, int *list, double *buf, void ComputeBondLocal::unpack_forward_comm(int n, int first, double *buf) { - int i,m,last; + int i, m, last; double **v = atom->v; @@ -504,11 +540,11 @@ void ComputeBondLocal::reallocate(int n) if (nvalues == 1) { memory->destroy(vlocal); - memory->create(vlocal,nmax,"bond/local:vector_local"); + memory->create(vlocal, nmax, "bond/local:vector_local"); vector_local = vlocal; } else { memory->destroy(alocal); - memory->create(alocal,nmax,nvalues,"bond/local:array_local"); + memory->create(alocal, nmax, nvalues, "bond/local:array_local"); array_local = alocal; } } @@ -519,6 +555,6 @@ void ComputeBondLocal::reallocate(int n) double ComputeBondLocal::memory_usage() { - double bytes = (double)nmax*nvalues * sizeof(double); + double bytes = (double) nmax * nvalues * sizeof(double); return bytes; } diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp index 894d0e33e4..bd1e389e85 100644 --- a/src/compute_dihedral_local.cpp +++ b/src/compute_dihedral_local.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -13,19 +12,21 @@ ------------------------------------------------------------------------- */ #include "compute_dihedral_local.h" -#include -#include + #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "update.h" #include "domain.h" +#include "error.h" #include "force.h" #include "input.h" -#include "variable.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "molecule.h" +#include "update.h" +#include "variable.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; @@ -37,14 +38,13 @@ enum { PHI, VARIABLE }; /* ---------------------------------------------------------------------- */ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - bstyle(nullptr), vvar(nullptr), pstr(nullptr), vstr(nullptr), vlocal(nullptr), alocal(nullptr) + Compute(lmp, narg, arg), bstyle(nullptr), vvar(nullptr), pstr(nullptr), vstr(nullptr), + vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute dihedral/local command"); + if (narg < 4) error->all(FLERR, "Illegal compute dihedral/local command"); if (atom->avec->dihedrals_allow == 0) - error->all(FLERR, - "Compute dihedral/local used when dihedrals are not allowed"); + error->all(FLERR, "Compute dihedral/local used when dihedrals are not allowed"); local_flag = 1; @@ -52,7 +52,7 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : nvalues = narg - 3; bstyle = new int[nvalues]; - vstr = new char*[nvalues]; + vstr = new char *[nvalues]; vvar = new int[nvalues]; nvalues = 0; @@ -60,13 +60,14 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : int iarg; for (iarg = 3; iarg < narg; iarg++) { - if (strcmp(arg[iarg],"phi") == 0) { + if (strcmp(arg[iarg], "phi") == 0) { bstyle[nvalues++] = PHI; - } else if (strncmp(arg[iarg],"v_",2) == 0) { + } else if (strncmp(arg[iarg], "v_", 2) == 0) { bstyle[nvalues++] = VARIABLE; vstr[nvar] = utils::strdup(&arg[iarg][2]); nvar++; - } else break; + } else + break; } // optional args @@ -75,47 +76,45 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : pstr = nullptr; while (iarg < narg) { - if (strcmp(arg[iarg],"set") == 0) { + if (strcmp(arg[iarg], "set") == 0) { setflag = 1; - if (iarg+3 > narg) - error->all(FLERR,"Illegal compute dihedral/local command"); - if (strcmp(arg[iarg+1],"phi") == 0) { - delete [] pstr; - pstr = utils::strdup(arg[iarg+2]); - } else error->all(FLERR,"Illegal compute dihedral/local command"); + if (iarg + 3 > narg) error->all(FLERR, "Illegal compute dihedral/local command"); + if (strcmp(arg[iarg + 1], "phi") == 0) { + delete[] pstr; + pstr = utils::strdup(arg[iarg + 2]); + } else + error->all(FLERR, "Illegal compute dihedral/local command"); iarg += 3; - } else error->all(FLERR,"Illegal compute dihedral/local command"); + } else + error->all(FLERR, "Illegal compute dihedral/local command"); } // error check if (nvar) { - if (!setflag) - error->all(FLERR,"Compute dihedral/local variable requires a set variable"); + if (!setflag) error->all(FLERR, "Compute dihedral/local variable requires a set variable"); for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) - error->all(FLERR, - "Variable name for copute dihedral/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist"); if (!input->variable->equalstyle(vvar[i])) - error->all(FLERR,"Variable for compute dihedral/local is invalid style"); + error->all(FLERR, "Variable for compute dihedral/local is invalid style"); } if (pstr) { pvar = input->variable->find(pstr); - if (pvar < 0) - error->all(FLERR, - "Variable name for compute dihedral/local does not exist"); + if (pvar < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist"); if (!input->variable->internalstyle(pvar)) - error->all(FLERR,"Variable for compute dihedral/local is invalid style"); + error->all(FLERR, "Variable for compute dihedral/local is invalid style"); } } else if (setflag) - error->all(FLERR,"Compute dihedral/local set with no variable"); + error->all(FLERR, "Compute dihedral/local set with no variable"); // initialize output - if (nvalues == 1) size_local_cols = 0; - else size_local_cols = nvalues; + if (nvalues == 1) + size_local_cols = 0; + else + size_local_cols = nvalues; nmax = 0; vlocal = nullptr; @@ -126,12 +125,12 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : ComputeDihedralLocal::~ComputeDihedralLocal() { - delete [] bstyle; - for (int i = 0; i < nvar; i++) delete [] vstr[i]; - delete [] vstr; - delete [] vvar; + delete[] bstyle; + for (int i = 0; i < nvar; i++) delete[] vstr[i]; + delete[] vstr; + delete[] vvar; - delete [] pstr; + delete[] pstr; memory->destroy(vlocal); memory->destroy(alocal); @@ -142,21 +141,17 @@ ComputeDihedralLocal::~ComputeDihedralLocal() void ComputeDihedralLocal::init() { if (force->dihedral == nullptr) - error->all(FLERR,"No dihedral style is defined for compute dihedral/local"); + error->all(FLERR, "No dihedral style is defined for compute dihedral/local"); if (nvar) { for (int i = 0; i < nvar; i++) { vvar[i] = input->variable->find(vstr[i]); - if (vvar[i] < 0) - error->all(FLERR, - "Variable name for compute dihedral/local does not exist"); + if (vvar[i] < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist"); } if (pstr) { pvar = input->variable->find(pstr); - if (pvar < 0) - error->all(FLERR, - "Variable name for compute dihedral/local does not exist"); + if (pvar < 0) error->all(FLERR, "Variable name for compute dihedral/local does not exist"); } } @@ -191,11 +186,11 @@ void ComputeDihedralLocal::compute_local() int ComputeDihedralLocal::compute_dihedrals(int flag) { - int i,m,nd,atom1,atom2,atom3,atom4,imol,iatom,ivar; + int i, m, nd, atom1, atom2, atom3, atom4, imol, iatom, ivar; tagint tagprev; - double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; - double ax,ay,az,bx,by,bz,rasq,rbsq,rgsq,rg,ra2inv,rb2inv,rabinv; - double s,c,phi; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z, vb2xm, vb2ym, vb2zm; + double ax, ay, az, bx, by, bz, rasq, rbsq, rgsq, rg, ra2inv, rb2inv, rabinv; + double s, c, phi; double *ptr; double **x = atom->x; @@ -220,7 +215,8 @@ int ComputeDihedralLocal::compute_dihedrals(int flag) for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; - if (molecular == Atom::MOLECULAR) nd = num_dihedral[atom2]; + if (molecular == Atom::MOLECULAR) + nd = num_dihedral[atom2]; else { if (molindex[atom2] < 0) continue; imol = molindex[atom2]; @@ -237,9 +233,9 @@ int ComputeDihedralLocal::compute_dihedrals(int flag) } else { if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; tagprev = tag[atom2] - iatom - 1; - atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i]+tagprev); - atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i]+tagprev); - atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i]+tagprev); + atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i] + tagprev); + atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i] + tagprev); } if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; @@ -256,64 +252,66 @@ int ComputeDihedralLocal::compute_dihedrals(int flag) vb1x = x[atom1][0] - x[atom2][0]; vb1y = x[atom1][1] - x[atom2][1]; vb1z = x[atom1][2] - x[atom2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); + domain->minimum_image(vb1x, vb1y, vb1z); vb2x = x[atom3][0] - x[atom2][0]; vb2y = x[atom3][1] - x[atom2][1]; vb2z = x[atom3][2] - x[atom2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); + domain->minimum_image(vb2x, vb2y, vb2z); vb2xm = -vb2x; vb2ym = -vb2y; vb2zm = -vb2z; - domain->minimum_image(vb2xm,vb2ym,vb2zm); + domain->minimum_image(vb2xm, vb2ym, vb2zm); vb3x = x[atom4][0] - x[atom3][0]; vb3y = x[atom4][1] - x[atom3][1]; vb3z = x[atom4][2] - x[atom3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); + domain->minimum_image(vb3x, vb3y, vb3z); - ax = vb1y*vb2zm - vb1z*vb2ym; - ay = vb1z*vb2xm - vb1x*vb2zm; - az = vb1x*vb2ym - vb1y*vb2xm; - bx = vb3y*vb2zm - vb3z*vb2ym; - by = vb3z*vb2xm - vb3x*vb2zm; - bz = vb3x*vb2ym - vb3y*vb2xm; + ax = vb1y * vb2zm - vb1z * vb2ym; + ay = vb1z * vb2xm - vb1x * vb2zm; + az = vb1x * vb2ym - vb1y * vb2xm; + bx = vb3y * vb2zm - vb3z * vb2ym; + by = vb3z * vb2xm - vb3x * vb2zm; + bz = vb3x * vb2ym - vb3y * vb2xm; - rasq = ax*ax + ay*ay + az*az; - rbsq = bx*bx + by*by + bz*bz; - rgsq = vb2xm*vb2xm + vb2ym*vb2ym + vb2zm*vb2zm; + rasq = ax * ax + ay * ay + az * az; + rbsq = bx * bx + by * by + bz * bz; + rgsq = vb2xm * vb2xm + vb2ym * vb2ym + vb2zm * vb2zm; rg = sqrt(rgsq); ra2inv = rb2inv = 0.0; - if (rasq > 0) ra2inv = 1.0/rasq; - if (rbsq > 0) rb2inv = 1.0/rbsq; - rabinv = sqrt(ra2inv*rb2inv); + if (rasq > 0) ra2inv = 1.0 / rasq; + if (rbsq > 0) rb2inv = 1.0 / rbsq; + rabinv = sqrt(ra2inv * rb2inv); - c = (ax*bx + ay*by + az*bz)*rabinv; - s = rg*rabinv*(ax*vb3x + ay*vb3y + az*vb3z); + c = (ax * bx + ay * by + az * bz) * rabinv; + s = rg * rabinv * (ax * vb3x + ay * vb3y + az * vb3z); if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; - phi = atan2(s,c); + phi = atan2(s, c); - if (nvalues == 1) ptr = &vlocal[m]; - else ptr = alocal[m]; + if (nvalues == 1) + ptr = &vlocal[m]; + else + ptr = alocal[m]; if (nvar) { ivar = 0; - if (pstr) input->variable->internal_set(pvar,phi); + if (pstr) input->variable->internal_set(pvar, phi); } for (int n = 0; n < nvalues; n++) { switch (bstyle[n]) { - case PHI: - ptr[n] = 180.0*phi/MY_PI; - break; - case VARIABLE: - ptr[n] = input->variable->compute_equal(vvar[ivar]); - ivar++; - break; + case PHI: + ptr[n] = 180.0 * phi / MY_PI; + break; + case VARIABLE: + ptr[n] = input->variable->compute_equal(vvar[ivar]); + ivar++; + break; } } @@ -334,11 +332,11 @@ void ComputeDihedralLocal::reallocate(int n) if (nvalues == 1) { memory->destroy(vlocal); - memory->create(vlocal,nmax,"dihedral/local:vector_local"); + memory->create(vlocal, nmax, "dihedral/local:vector_local"); vector_local = vlocal; } else { memory->destroy(alocal); - memory->create(alocal,nmax,nvalues,"dihedral/local:array_local"); + memory->create(alocal, nmax, nvalues, "dihedral/local:array_local"); array_local = alocal; } } @@ -349,6 +347,6 @@ void ComputeDihedralLocal::reallocate(int n) double ComputeDihedralLocal::memory_usage() { - double bytes = (double)nmax*nvalues * sizeof(double); + double bytes = (double) nmax * nvalues * sizeof(double); return bytes; } diff --git a/src/compute_improper_local.cpp b/src/compute_improper_local.cpp index a58f4f4d0d..ff8f069a03 100644 --- a/src/compute_improper_local.cpp +++ b/src/compute_improper_local.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -13,36 +12,35 @@ ------------------------------------------------------------------------- */ #include "compute_improper_local.h" -#include -#include + #include "atom.h" #include "atom_vec.h" -#include "molecule.h" -#include "update.h" #include "domain.h" +#include "error.h" #include "force.h" #include "math_const.h" #include "memory.h" -#include "error.h" +#include "molecule.h" +#include "update.h" + +#include +#include using namespace LAMMPS_NS; using namespace MathConst; static constexpr int DELTA = 10000; - -static constexpr double SMALL = 0.001; +static constexpr double SMALL = 0.001; /* ---------------------------------------------------------------------- */ ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), - vlocal(nullptr), alocal(nullptr) + Compute(lmp, narg, arg), vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR,"Illegal compute improper/local command"); + if (narg < 4) error->all(FLERR, "Illegal compute improper/local command"); if (atom->avec->impropers_allow == 0) - error->all(FLERR, - "Compute improper/local used when impropers are not allowed"); + error->all(FLERR, "Compute improper/local used when impropers are not allowed"); local_flag = 1; nvalues = narg - 3; @@ -50,12 +48,16 @@ ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) : nvalues = 0; for (int iarg = 3; iarg < narg; iarg++) { - if (strcmp(arg[iarg],"chi") == 0) cflag = nvalues++; - else error->all(FLERR,"Invalid keyword in compute improper/local command"); + if (strcmp(arg[iarg], "chi") == 0) + cflag = nvalues++; + else + error->all(FLERR, "Invalid keyword in compute improper/local command"); } - if (nvalues == 1) size_local_cols = 0; - else size_local_cols = nvalues; + if (nvalues == 1) + size_local_cols = 0; + else + size_local_cols = nvalues; nmax = 0; vlocal = nullptr; @@ -75,7 +77,7 @@ ComputeImproperLocal::~ComputeImproperLocal() void ComputeImproperLocal::init() { if (force->improper == nullptr) - error->all(FLERR,"No improper style is defined for compute improper/local"); + error->all(FLERR, "No improper style is defined for compute improper/local"); // do initial memory allocation so that memory_usage() is correct @@ -108,11 +110,11 @@ void ComputeImproperLocal::compute_local() int ComputeImproperLocal::compute_impropers(int flag) { - int i,m,n,ni,atom1,atom2,atom3,atom4,imol,iatom; + int i, m, n, ni, atom1, atom2, atom3, atom4, imol, iatom; tagint tagprev; - double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z; - double ss1,ss2,ss3,r1,r2,r3,c0,c1,c2,s1,s2; - double s12,c; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z; + double ss1, ss2, ss3, r1, r2, r3, c0, c1, c2, s1, s2; + double s12, c; double *cbuf; double **x = atom->x; @@ -135,8 +137,10 @@ int ComputeImproperLocal::compute_impropers(int flag) if (nvalues == 1) { if (cflag >= 0) cbuf = vlocal; } else { - if (cflag >= 0 && alocal) cbuf = &alocal[0][cflag]; - else cbuf = nullptr; + if (cflag >= 0 && alocal) + cbuf = &alocal[0][cflag]; + else + cbuf = nullptr; } } @@ -144,7 +148,8 @@ int ComputeImproperLocal::compute_impropers(int flag) for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; - if (molecular == Atom::MOLECULAR) ni = num_improper[atom2]; + if (molecular == Atom::MOLECULAR) + ni = num_improper[atom2]; else { if (molindex[atom2] < 0) continue; imol = molindex[atom2]; @@ -161,9 +166,9 @@ int ComputeImproperLocal::compute_impropers(int flag) } else { if (tag[atom2] != onemols[imol]->improper_atom2[atom2][i]) continue; tagprev = tag[atom2] - iatom - 1; - atom1 = atom->map(onemols[imol]->improper_atom1[atom2][i]+tagprev); - atom3 = atom->map(onemols[imol]->improper_atom3[atom2][i]+tagprev); - atom4 = atom->map(onemols[imol]->improper_atom4[atom2][i]+tagprev); + atom1 = atom->map(onemols[imol]->improper_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->improper_atom3[atom2][i] + tagprev); + atom4 = atom->map(onemols[imol]->improper_atom4[atom2][i] + tagprev); } if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; @@ -178,21 +183,21 @@ int ComputeImproperLocal::compute_impropers(int flag) vb1x = x[atom1][0] - x[atom2][0]; vb1y = x[atom1][1] - x[atom2][1]; vb1z = x[atom1][2] - x[atom2][2]; - domain->minimum_image(vb1x,vb1y,vb1z); + domain->minimum_image(vb1x, vb1y, vb1z); vb2x = x[atom3][0] - x[atom2][0]; vb2y = x[atom3][1] - x[atom2][1]; vb2z = x[atom3][2] - x[atom2][2]; - domain->minimum_image(vb2x,vb2y,vb2z); + domain->minimum_image(vb2x, vb2y, vb2z); vb3x = x[atom4][0] - x[atom3][0]; vb3y = x[atom4][1] - x[atom3][1]; vb3z = x[atom4][2] - x[atom3][2]; - domain->minimum_image(vb3x,vb3y,vb3z); + domain->minimum_image(vb3x, vb3y, vb3z); - ss1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); - ss2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); - ss3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + ss1 = 1.0 / (vb1x * vb1x + vb1y * vb1y + vb1z * vb1z); + ss2 = 1.0 / (vb2x * vb2x + vb2y * vb2y + vb2z * vb2z); + ss3 = 1.0 / (vb3x * vb3x + vb3y * vb3y + vb3z * vb3z); r1 = sqrt(ss1); r2 = sqrt(ss2); @@ -202,20 +207,20 @@ int ComputeImproperLocal::compute_impropers(int flag) c1 = (vb1x * vb2x + vb1y * vb2y + vb1z * vb2z) * r1 * r2; c2 = -(vb3x * vb2x + vb3y * vb2y + vb3z * vb2z) * r3 * r2; - s1 = 1.0 - c1*c1; + s1 = 1.0 - c1 * c1; if (s1 < SMALL) s1 = SMALL; s1 = 1.0 / s1; - s2 = 1.0 - c2*c2; + s2 = 1.0 - c2 * c2; if (s2 < SMALL) s2 = SMALL; s2 = 1.0 / s2; - s12 = sqrt(s1*s2); - c = (c1*c2 + c0) * s12; + s12 = sqrt(s1 * s2); + c = (c1 * c2 + c0) * s12; if (c > 1.0) c = 1.0; if (c < -1.0) c = -1.0; - cbuf[n] = 180.0*acos(c)/MY_PI; + cbuf[n] = 180.0 * acos(c) / MY_PI; } n += nvalues; } @@ -237,11 +242,11 @@ void ComputeImproperLocal::reallocate(int n) if (nvalues == 1) { memory->destroy(vlocal); - memory->create(vlocal,nmax,"improper/local:vector_local"); + memory->create(vlocal, nmax, "improper/local:vector_local"); vector_local = vlocal; } else { memory->destroy(alocal); - memory->create(alocal,nmax,nvalues,"improper/local:array_local"); + memory->create(alocal, nmax, nvalues, "improper/local:array_local"); array_local = alocal; } } @@ -252,6 +257,6 @@ void ComputeImproperLocal::reallocate(int n) double ComputeImproperLocal::memory_usage() { - double bytes = (double)nmax*nvalues * sizeof(double); + double bytes = (double) nmax * nvalues * sizeof(double); return bytes; } diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index 351b499468..57f15264f0 100644 --- a/src/compute_pair_local.cpp +++ b/src/compute_pair_local.cpp @@ -38,7 +38,7 @@ enum { TYPE, RADIUS }; ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), pstyle(nullptr), pindex(nullptr), vlocal(nullptr), alocal(nullptr) { - if (narg < 4) error->all(FLERR, "Illegal compute pair/local command"); + if (narg < 4) utils::missing_cmd_args(FLERR, "compute pair/local", error); local_flag = 1; nvalues = narg - 3; @@ -97,7 +97,7 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : // error check if (cutstyle == RADIUS && !atom->radius_flag) - error->all(FLERR, "This compute pair/local requires atom attribute radius"); + error->all(FLERR, "Compute pair/local with ID {} requires atom attribute radius", id); // set singleflag if need to call pair->single() diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index 64f3859117..872f4173ec 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -39,7 +39,7 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), vlocal(nullptr), alocal(nullptr), indices(nullptr), pack_choice(nullptr) { - if (narg < 4) error->all(FLERR, "Illegal compute property/local command"); + if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/local", error); local_flag = 1; nvalues = narg - 3; @@ -202,25 +202,23 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : while (iarg < narg) { if (strcmp(arg[iarg], "cutoff") == 0) { - if (iarg + 2 > narg) error->all(FLERR, "Illegal compute property/local command"); + if (iarg + 2 > narg) utils::missing_cmd_args(FLERR, "compute property/local cutoff", error); if (strcmp(arg[iarg + 1], "type") == 0) cutstyle = TYPE; else if (strcmp(arg[iarg + 1], "radius") == 0) cutstyle = RADIUS; else - error->all(FLERR, "Illegal compute property/local command"); + error->all(FLERR, "Unknown compute property/local cutoff keyword: {}", arg[iarg + 1]); iarg += 2; } else - error->all(FLERR, "Illegal compute property/local command"); + error->all(FLERR, "Unknown compute property/local keyword: {}", arg[iarg]); } // error check if (atom->molecular == 2 && (kindflag == BOND || kindflag == ANGLE || kindflag == DIHEDRAL || kindflag == IMPROPER)) - error->all(FLERR, - "Compute property/local does not (yet) work " - "with atom_style template"); + error->all(FLERR, "Compute property/local does not (yet) work with atom_style template"); if (kindflag == BOND && atom->avec->bonds_allow == 0) error->all(FLERR, "Compute property/local for property that isn't allocated"); @@ -231,7 +229,7 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : if (kindflag == IMPROPER && atom->avec->impropers_allow == 0) error->all(FLERR, "Compute property/local for property that isn't allocated"); if (cutstyle == RADIUS && !atom->radius_flag) - error->all(FLERR, "Compute property/local requires atom attribute radius"); + error->all(FLERR, "Compute property/local with ID {} requires atom attribute radius", id); nmax = 0; vlocal = nullptr; @@ -255,8 +253,6 @@ void ComputePropertyLocal::init() if (kindflag == NEIGH || kindflag == PAIR) { if (force->pair == nullptr) error->all(FLERR, "No pair style is defined for compute property/local"); - if (force->pair->single_enable == 0) - error->all(FLERR, "Pair style does not support compute property/local"); } // for NEIGH/PAIR need an occasional half neighbor list diff --git a/src/fix_nve.cpp b/src/fix_nve.cpp index 30c1f24867..58a5677b92 100644 --- a/src/fix_nve.cpp +++ b/src/fix_nve.cpp @@ -15,6 +15,7 @@ #include "fix_nve.h" #include "atom.h" +#include "error.h" #include "force.h" #include "respa.h" #include "update.h" @@ -28,7 +29,11 @@ FixNVE::FixNVE(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg) { if (!utils::strmatch(style,"^nve/sphere") && narg < 3) - utils::missing_cmd_args(FLERR, "fix nve", error); + utils::missing_cmd_args(FLERR, fmt::format("fix {}", style), error); + + auto plain_style = utils::strip_style_suffix(style, lmp); + if (utils::strmatch(plain_style, "^nve$") && narg > 3) + error->all(FLERR, "Unsupported additional arguments for fix {}", style); dynamic_group_allow = 1; time_integrate = 1; diff --git a/src/pair_hybrid.h b/src/pair_hybrid.h index b59cc82cb4..995de37e59 100644 --- a/src/pair_hybrid.h +++ b/src/pair_hybrid.h @@ -28,6 +28,7 @@ namespace LAMMPS_NS { class PairHybrid : public Pair { friend class AtomVecDielectric; friend class ComputeSpin; + friend class FixAtomSwap; friend class FixGPU; friend class FixIntel; friend class FixNVESpin; diff --git a/tools/lammps-gui/chartviewer.cpp b/tools/lammps-gui/chartviewer.cpp index eb5444de8e..8543da166a 100644 --- a/tools/lammps-gui/chartviewer.cpp +++ b/tools/lammps-gui/chartviewer.cpp @@ -423,27 +423,6 @@ void ChartViewer::add_data(int step, double data) { if (last_step < step) { last_step = step; - - // do not add data that deviates by more than 4 sigma from the average - // over the last 5 to 20 data items. this is a hack to work around - // getting corrupted data from lammps_get_last_thermo() - const auto &points = series->points(); - const auto count = points.count(); - if (count > 4) { - double ysum = 0.0; - double ysumsq = 0.0; - int first = count - 20; - if (first < 0) first = 0; - for (int i = first; i < count; ++i) { - double val = points[i].y(); - ysum += val; - ysumsq += val * val; - } - const double num = count - first; - const double avg = ysum / num; - const double avgsq = ysumsq / num; - if (fabs(data - avg) > (5.0 * sqrt(avgsq - (avg * avg)))) return; - } series->append(step, data); QSettings settings; diff --git a/tools/lammps-gui/lammpsgui.cpp b/tools/lammps-gui/lammpsgui.cpp index fe6b8c5391..fb38b6fc92 100644 --- a/tools/lammps-gui/lammpsgui.cpp +++ b/tools/lammps-gui/lammpsgui.cpp @@ -988,6 +988,7 @@ void LammpsGui::logupdate() if (logwindow) { const auto text = capturer->GetChunk(); if (text.size() > 0) { + logwindow->moveCursor(QTextCursor::End); logwindow->insertPlainText(text.c_str()); logwindow->moveCursor(QTextCursor::End); logwindow->textCursor().deleteChar(); @@ -2016,8 +2017,8 @@ bool LammpsGui::eventFilter(QObject *watched, QEvent *event) } // LAMMPS geturl command with current location of the input and solution files on the web -static const QString geturl = "geturl https://raw.githubusercontent.com/akohlmey/" - "lammps-tutorials-inputs/main/tutorial%1/%2 output %2 verify no"; +static const QString geturl = "geturl https://raw.githubusercontent.com/lammpstutorials/" + "lammpstutorials-article/refs/heads/main/files/tutorial%1/%2 output %2 verify no"; void LammpsGui::setup_tutorial(int tutno, const QString &dir, bool purgedir, bool getsolution) {