diff --git a/doc/src/compute_stress_atom.rst b/doc/src/compute_stress_atom.rst index 16afd49548..9ebb51a3eb 100644 --- a/doc/src/compute_stress_atom.rst +++ b/doc/src/compute_stress_atom.rst @@ -33,32 +33,31 @@ Examples Description """"""""""" -Define a computation that computes per-atom stress -tensor for each atom in a group. In case of compute *stress/atom*, -the tensor for each atom is symmetric with 6 -components and is stored as a 6-element vector in the following order: -:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`. -In case of compute *centroid/stress/atom*, -the tensor for each atom is asymmetric with 9 components -and is stored as a 9-element vector in the following order: -:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, :math:`yz`, -:math:`yx`, :math:`zx`, :math:`zy`. -See the :doc:`compute pressure ` command if you want the stress tensor +Define a computation that computes per-atom stress tensor for each +atom in a group. In case of compute *stress/atom*, the tensor for +each atom is symmetric with 6 components and is stored as a 6-element +vector in the following order: :math:`xx`, :math:`yy`, :math:`zz`, +:math:`xy`, :math:`xz`, :math:`yz`. In case of compute +*centroid/stress/atom*, the tensor for each atom is asymmetric with 9 +components and is stored as a 9-element vector in the following order: +:math:`xx`, :math:`yy`, :math:`zz`, :math:`xy`, :math:`xz`, +:math:`yz`, :math:`yx`, :math:`zx`, :math:`zy`. See the :doc:`compute +pressure ` command if you want the stress tensor (pressure) of the entire system. -The stress tensor for atom :math:`I` is given by the following formula, -where :math:`a` and :math:`b` take on values :math:`x`, :math:`y`, :math:`z` -to generate the components of the tensor: +The stress tensor for atom :math:`I` is given by the following +formula, where :math:`a` and :math:`b` take on values :math:`x`, +:math:`y`, :math:`z` to generate the components of the tensor: .. math:: S_{ab} = - m v_a v_b - W_{ab} -The first term is a kinetic energy contribution for atom :math:`I`. See -details below on how the specified *temp-ID* can affect the velocities -used in this calculation. The second term is the virial -contribution due to intra and intermolecular interactions, -where the exact computation details are determined by the compute style. +The first term is a kinetic energy contribution for atom :math:`I`. +See details below on how the specified *temp-ID* can affect the +velocities used in this calculation. The second term is the virial +contribution due to intra and intermolecular interactions, where the +exact computation details are determined by the compute style. In case of compute *stress/atom*, the virial contribution is: @@ -68,29 +67,26 @@ In case of compute *stress/atom*, the virial contribution is: & + \frac{1}{3} \sum_{n = 1}^{N_a} (r_{1_a} F_{1_b} + r_{2_a} F_{2_b} + r_{3_a} F_{3_b}) + \frac{1}{4} \sum_{n = 1}^{N_d} (r_{1_a} F_{1_b} + r_{2_a} F_{2_b} + r_{3_a} F_{3_b} + r_{4_a} F_{4_b}) \\ & + \frac{1}{4} \sum_{n = 1}^{N_i} (r_{1_a} F_{1_b} + r_{2_a} F_{2_b} + r_{3_a} F_{3_b} + r_{4_a} F_{4_b}) + {\rm Kspace}(r_{i_a},F_{i_b}) + \sum_{n = 1}^{N_f} r_{i_a} F_{i_b} -The first term is a pairwise energy -contribution where :math:`n` loops over the :math:`N_p` -neighbors of atom :math:`I`, :math:`\mathbf{r}_1` and :math:`\mathbf{r}_2` -are the positions of the 2 atoms in the pairwise interaction, -and :math:`\mathbf{F}_1` and :math:`\mathbf{F}_2` are the forces -on the 2 atoms resulting from the pairwise interaction. -The second term is a bond contribution of -similar form for the :math:`N_b` bonds which atom :math:`I` is part of. -There are similar terms for the :math:`N_a` angle, :math:`N_d` dihedral, -and :math:`N_i` improper interactions atom :math:`I` is part of. -There is also a term for the KSpace -contribution from long-range Coulombic interactions, if defined. -Finally, there is a term for the :math:`N_f` :doc:`fixes ` that apply -internal constraint forces to atom :math:`I`. Currently, only the -:doc:`fix shake ` and :doc:`fix rigid ` commands -contribute to this term. -As the coefficients in the formula imply, a virial contribution -produced by a small set of atoms (e.g. 4 atoms in a dihedral or 3 -atoms in a Tersoff 3-body interaction) is assigned in equal portions -to each atom in the set. E.g. 1/4 of the dihedral virial to each of -the 4 atoms, or 1/3 of the fix virial due to SHAKE constraints applied -to atoms in a water molecule via the :doc:`fix shake ` -command. +The first term is a pairwise energy contribution where :math:`n` loops +over the :math:`N_p` neighbors of atom :math:`I`, :math:`\mathbf{r}_1` +and :math:`\mathbf{r}_2` are the positions of the 2 atoms in the +pairwise interaction, and :math:`\mathbf{F}_1` and +:math:`\mathbf{F}_2` are the forces on the 2 atoms resulting from the +pairwise interaction. The second term is a bond contribution of +similar form for the :math:`N_b` bonds which atom :math:`I` is part +of. There are similar terms for the :math:`N_a` angle, :math:`N_d` +dihedral, and :math:`N_i` improper interactions atom :math:`I` is part +of. There is also a term for the KSpace contribution from long-range +Coulombic interactions, if defined. Finally, there is a term for the +:math:`N_f` :doc:`fixes ` that apply internal constraint forces +to atom :math:`I`. Currently, only the :doc:`fix shake ` +and :doc:`fix rigid ` commands contribute to this term. As +the coefficients in the formula imply, a virial contribution produced +by a small set of atoms (e.g. 4 atoms in a dihedral or 3 atoms in a +Tersoff 3-body interaction) is assigned in equal portions to each atom +in the set. E.g. 1/4 of the dihedral virial to each of the 4 atoms, +or 1/3 of the fix virial due to SHAKE constraints applied to atoms in +a water molecule via the :doc:`fix shake ` command. In case of compute *centroid/stress/atom*, the virial contribution is: @@ -99,71 +95,69 @@ In case of compute *centroid/stress/atom*, the virial contribution is: W_{ab} & = \sum_{n = 1}^{N_p} r_{I0_a} F_{I_b} + \sum_{n = 1}^{N_b} r_{I0_a} F_{I_b} + \sum_{n = 1}^{N_a} r_{I0_a} F_{I_b} + \sum_{n = 1}^{N_d} r_{I0_a} F_{I_b} + \sum_{n = 1}^{N_i} r_{I0_a} F_{I_b} \\ & + {\rm Kspace}(r_{i_a},F_{i_b}) + \sum_{n = 1}^{N_f} r_{i_a} F_{i_b} -As with compute *stress/atom*, the first, second, third, fourth and fifth terms -are pairwise, bond, angle, dihedral and improper contributions, -but instead of assigning the virial contribution equally to each atom, -only the force :math:`\mathbf{F}_I` acting on atom :math:`I` -due to the interaction and the relative -position :math:`\mathbf{r}_{I0}` of the atom :math:`I` to the geometric center -of the interacting atoms, i.e. centroid, is used. -As the geometric center is different -for each interaction, the :math:`\mathbf{r}_{I0}` also differs. -The sixth and seventh terms, Kspace and :doc:`fix ` contribution -respectively, are computed identical to compute *stress/atom*. -Although the total system virial is the same as compute *stress/atom*, -compute *centroid/stress/atom* is know to result in more consistent -heat flux values for angle, dihedrals and improper contributions -when computed via :doc:`compute heat/flux `. +As with compute *stress/atom*, the first, second, third, fourth and +fifth terms are pairwise, bond, angle, dihedral and improper +contributions, but instead of assigning the virial contribution +equally to each atom, only the force :math:`\mathbf{F}_I` acting on +atom :math:`I` due to the interaction and the relative position +:math:`\mathbf{r}_{I0}` of the atom :math:`I` to the geometric center +of the interacting atoms, i.e. centroid, is used. As the geometric +center is different for each interaction, the :math:`\mathbf{r}_{I0}` +also differs. The sixth and seventh terms, Kspace and :doc:`fix +` contribution respectively, are computed identical to compute +*stress/atom*. Although the total system virial is the same as +compute *stress/atom*, compute *centroid/stress/atom* is know to +result in more consistent heat flux values for angle, dihedrals and +improper contributions when computed via :doc:`compute heat/flux +`. -If no extra keywords are listed, the kinetic contribution -all of the virial contribution terms are -included in the per-atom stress tensor. If any extra keywords are -listed, only those terms are summed to compute the tensor. The -*virial* keyword means include all terms except the kinetic energy -*ke*\ . +If no extra keywords are listed, the kinetic contribution all of the +virial contribution terms are included in the per-atom stress tensor. +If any extra keywords are listed, only those terms are summed to +compute the tensor. The *virial* keyword means include all terms +except the kinetic energy *ke*\ . Note that the stress for each atom is due to its interaction with all other atoms in the simulation, not just with other atoms in the group. -Details of how compute *stress/atom* obtains the virial for individual atoms for -either pairwise or many-body potentials, and including the effects of -periodic boundary conditions is discussed in :ref:`(Thompson) `. -The basic idea for many-body potentials is to treat each component of -the force computation between a small cluster of atoms in the same -manner as in the formula above for bond, angle, dihedral, etc -interactions. Namely the quantity :math:`\mathbf{r} \cdot \mathbf{F}` -is summed over the atoms in -the interaction, with the :math:`r` vectors unwrapped by periodic boundaries -so that the cluster of atoms is close together. The total +Details of how compute *stress/atom* obtains the virial for individual +atoms for either pairwise or many-body potentials, and including the +effects of periodic boundary conditions is discussed in +:ref:`(Thompson) `. The basic idea for many-body +potentials is to treat each component of the force computation between +a small cluster of atoms in the same manner as in the formula above +for bond, angle, dihedral, etc interactions. Namely the quantity +:math:`\mathbf{r} \cdot \mathbf{F}` is summed over the atoms in the +interaction, with the :math:`r` vectors unwrapped by periodic +boundaries so that the cluster of atoms is close together. The total contribution for the cluster interaction is divided evenly among those -atoms. Details of how compute *centroid/stress/atom* obtains -the virial for individual atoms -is given in :ref:`(Surblys) `, -where the idea is that the virial of the atom :math:`I` -is the result of only the force :math:`\mathbf{F}_I` on the atom due -to the interaction -and its positional vector :math:`\mathbf{r}_{I0}`, -relative to the geometric center of the -interacting atoms, regardless of the number of participating atoms. -The periodic boundary treatment is identical to +atoms. + +Details of how compute *centroid/stress/atom* obtains the virial for +individual atoms is given in :ref:`(Surblys) `, where the +idea is that the virial of the atom :math:`I` is the result of only +the force :math:`\mathbf{F}_I` on the atom due to the interaction and +its positional vector :math:`\mathbf{r}_{I0}`, relative to the +geometric center of the interacting atoms, regardless of the number of +participating atoms. The periodic boundary treatment is identical to that of compute *stress/atom*, and both of them reduce to identical -expressions for two-body interactions, -i.e. computed values for contributions from bonds and two-body pair styles, -such as :doc:`Lennard-Jones `, will be the same, -while contributions from angles, dihedrals and impropers will be different. +expressions for two-body interactions, i.e. computed values for +contributions from bonds and two-body pair styles, such as +:doc:`Lennard-Jones `, will be the same, while contributions +from angles, dihedrals and impropers will be different. The :doc:`dihedral_style charmm ` style calculates pairwise interactions between 1-4 atoms. The virial contribution of these terms is included in the pair virial, not the dihedral virial. The KSpace contribution is calculated using the method in -:ref:`(Heyes) ` for the Ewald method and by the methodology described -in :ref:`(Sirk) ` for PPPM. The choice of KSpace solver is specified -by the :doc:`kspace_style pppm ` command. Note that for -PPPM, the calculation requires 6 extra FFTs each timestep that -per-atom stress is calculated. Thus it can significantly increase the -cost of the PPPM calculation if it is needed on a large fraction of -the simulation timesteps. +:ref:`(Heyes) ` for the Ewald method and by the methodology +described in :ref:`(Sirk) ` for PPPM. The choice of KSpace +solver is specified by the :doc:`kspace_style pppm ` +command. Note that for PPPM, the calculation requires 6 extra FFTs +each timestep that per-atom stress is calculated. Thus it can +significantly increase the cost of the PPPM calculation if it is +needed on a large fraction of the simulation timesteps. The *temp-ID* argument can be used to affect the per-atom velocities used in the kinetic energy contribution to the total stress. If the @@ -189,10 +183,10 @@ See the :doc:`compute voronoi/atom ` command for one possible way to estimate a per-atom volume. Thus, if the diagonal components of the per-atom stress tensor are -summed for all atoms in the system and the sum is divided by :math:`dV`, where -:math:`d` = dimension and :math:`V` is the volume of the system, -the result should be :math:`-P`, where :math:`P` -is the total pressure of the system. +summed for all atoms in the system and the sum is divided by +:math:`dV`, where :math:`d` = dimension and :math:`V` is the volume of +the system, the result should be :math:`-P`, where :math:`P` is the +total pressure of the system. These lines in an input script for a 3d system should yield that result. I.e. the last 2 columns of thermo output will be the same: @@ -207,33 +201,43 @@ result. I.e. the last 2 columns of thermo output will be the same: .. note:: The per-atom stress does not include any Lennard-Jones tail - corrections to the pressure added by the :doc:`pair_modify tail yes ` command, since those are contributions to the - global system pressure. + corrections to the pressure added by the :doc:`pair_modify tail yes + ` command, since those are contributions to the global + system pressure. Output info """"""""""" -This compute *stress/atom* calculates a per-atom array with 6 columns, which can be -accessed by indices 1-6 by any command that uses per-atom values from -a compute as input. -The compute *centroid/stress/atom* produces a per-atom array with 9 columns, -but otherwise can be used in an identical manner to compute *stress/atom*. -See the :doc:`Howto output ` doc page -for an overview of LAMMPS output options. +Compute *stress/atom* calculates a per-atom array with 6 columns, +which can be accessed by indices 1-6 by any command that uses per-atom +values from a compute as input. Compute *centroid/stress/atom* +produces a per-atom array with 9 columns, but otherwise can be used in +an identical manner to compute *stress/atom*. See the :doc:`Howto +output ` doc page for an overview of LAMMPS output +options. -The per-atom array values will be in pressure\*volume -:doc:`units ` as discussed above. +The per-atom array values will be in pressure\*volume :doc:`units +` as discussed above. Restrictions """""""""""" -Currently (Spring 2020), compute *centroid/stress/atom* does not support -pair styles with many-body interactions, such as :doc:`Tersoff -`, or pair styles with long-range Coulomb interactions. -LAMMPS will generate an error in such cases. In principal, equivalent -formulation to that of angle, dihedral and improper contributions in the -virial :math:`W_{ab}` formula can also be applied to the many-body pair -styles, and is planned in the future. +Currently, compute *centroid/stress/atom* does not support pair styles +with many-body interactions (:doc:`EAM ` is an exception, +since its computations are performed pairwise), nor granular pair +styles with pairwise forces which are not aligned with the vector +between the pair of particles. All bond styles are supported. All +angle, dihedral, improper styles are supported with the exception of +USER-INTEL and KOKKOS variants of specific styles. It also does not +support models with long-range Coulombic or dispersion forces, +i.e. the kspace_style command in LAMMPS. It also does not support the +following fixes which add rigid-body constraints: :doc:`fix shake +`, :doc:`fix rattle `, :doc:`fix rigid +`, :doc:`fix rigid/small `. + +LAMMPS will generate an error if one of these options is included in +your model. Extension of centroid stress calculations to these force +and fix styles is planned for the futre. Related commands """""""""""""""" diff --git a/src/CLASS2/pair_lj_class2.cpp b/src/CLASS2/pair_lj_class2.cpp index 79e0c02fd7..3727b17493 100644 --- a/src/CLASS2/pair_lj_class2.cpp +++ b/src/CLASS2/pair_lj_class2.cpp @@ -35,7 +35,7 @@ PairLJClass2::PairLJClass2(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; writedata = 1; - centroidstressflag = 1; + centroidstressflag = CENTROID_SAME; } /* ---------------------------------------------------------------------- */ diff --git a/src/CLASS2/pair_lj_class2_coul_cut.cpp b/src/CLASS2/pair_lj_class2_coul_cut.cpp index c27dd0aee3..a767b70e82 100644 --- a/src/CLASS2/pair_lj_class2_coul_cut.cpp +++ b/src/CLASS2/pair_lj_class2_coul_cut.cpp @@ -33,7 +33,7 @@ using namespace MathConst; PairLJClass2CoulCut::PairLJClass2CoulCut(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; + centroidstressflag = CENTROID_SAME; } /* ---------------------------------------------------------------------- */ diff --git a/src/GRANULAR/pair_gran_hooke_history.cpp b/src/GRANULAR/pair_gran_hooke_history.cpp index 9d775c7f5b..7b7586d355 100644 --- a/src/GRANULAR/pair_gran_hooke_history.cpp +++ b/src/GRANULAR/pair_gran_hooke_history.cpp @@ -34,7 +34,6 @@ #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ @@ -43,6 +42,7 @@ PairGranHookeHistory::PairGranHookeHistory(LAMMPS *lmp) : Pair(lmp) { single_enable = 1; no_virial_fdotr_compute = 1; + centroidstressflag = CENTROID_NOTAVAIL; history = 1; size_history = 3; diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 1d2c7c3627..fe88d0755f 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -66,6 +66,7 @@ PairGranular::PairGranular(LAMMPS *lmp) : Pair(lmp) { single_enable = 1; no_virial_fdotr_compute = 1; + centroidstressflag = CENTROID_NOTAVAIL; single_extra = 12; svector = new double[single_extra]; diff --git a/src/KOKKOS/angle_charmm_kokkos.cpp b/src/KOKKOS/angle_charmm_kokkos.cpp index 6451b8c393..b363062c20 100644 --- a/src/KOKKOS/angle_charmm_kokkos.cpp +++ b/src/KOKKOS/angle_charmm_kokkos.cpp @@ -42,6 +42,8 @@ AngleCharmmKokkos::AngleCharmmKokkos(LAMMPS *lmp) : AngleCharmm(lmp) execution_space = ExecutionSpaceFromDevice::space; datamask_read = X_MASK | F_MASK | ENERGY_MASK | VIRIAL_MASK; datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; + + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/angle_class2_kokkos.cpp b/src/KOKKOS/angle_class2_kokkos.cpp index d4866201e7..a640289c76 100644 --- a/src/KOKKOS/angle_class2_kokkos.cpp +++ b/src/KOKKOS/angle_class2_kokkos.cpp @@ -42,6 +42,8 @@ AngleClass2Kokkos::AngleClass2Kokkos(LAMMPS *lmp) : AngleClass2(lmp) execution_space = ExecutionSpaceFromDevice::space; datamask_read = X_MASK | F_MASK | ENERGY_MASK | VIRIAL_MASK; datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; + + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/angle_cosine_kokkos.cpp b/src/KOKKOS/angle_cosine_kokkos.cpp index 6dbba1f86e..9f70df3baf 100644 --- a/src/KOKKOS/angle_cosine_kokkos.cpp +++ b/src/KOKKOS/angle_cosine_kokkos.cpp @@ -42,6 +42,8 @@ AngleCosineKokkos::AngleCosineKokkos(LAMMPS *lmp) : AngleCosine(lmp) execution_space = ExecutionSpaceFromDevice::space; datamask_read = X_MASK | F_MASK | ENERGY_MASK | VIRIAL_MASK; datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; + + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/angle_harmonic_kokkos.cpp b/src/KOKKOS/angle_harmonic_kokkos.cpp index 42d42dd6ca..088f1072dc 100644 --- a/src/KOKKOS/angle_harmonic_kokkos.cpp +++ b/src/KOKKOS/angle_harmonic_kokkos.cpp @@ -42,6 +42,8 @@ AngleHarmonicKokkos::AngleHarmonicKokkos(LAMMPS *lmp) : AngleHarmoni execution_space = ExecutionSpaceFromDevice::space; datamask_read = X_MASK | F_MASK | ENERGY_MASK | VIRIAL_MASK; datamask_modify = F_MASK | ENERGY_MASK | VIRIAL_MASK; + + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/dihedral_charmm_kokkos.cpp b/src/KOKKOS/dihedral_charmm_kokkos.cpp index 86539bc090..92efd9d082 100644 --- a/src/KOKKOS/dihedral_charmm_kokkos.cpp +++ b/src/KOKKOS/dihedral_charmm_kokkos.cpp @@ -48,6 +48,8 @@ DihedralCharmmKokkos::DihedralCharmmKokkos(LAMMPS *lmp) : DihedralCh k_warning_flag = Kokkos::DualView("Dihedral:warning_flag"); d_warning_flag = k_warning_flag.template view(); h_warning_flag = k_warning_flag.h_view; + + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ @@ -76,7 +78,7 @@ void DihedralCharmmKokkos::compute(int eflag_in, int vflag_in) // insure pair->ev_tally() will use 1-4 virial contribution - if (weightflag && vflag_global == 2) + if (weightflag && vflag_global == VIRIAL_FDOTR) force->pair->vflag_either = force->pair->vflag_global = 1; // reallocate per-atom arrays if necessary diff --git a/src/KOKKOS/dihedral_class2_kokkos.cpp b/src/KOKKOS/dihedral_class2_kokkos.cpp index 4848b86756..59a20b26da 100644 --- a/src/KOKKOS/dihedral_class2_kokkos.cpp +++ b/src/KOKKOS/dihedral_class2_kokkos.cpp @@ -47,6 +47,8 @@ DihedralClass2Kokkos::DihedralClass2Kokkos(LAMMPS *lmp) : DihedralCl k_warning_flag = DAT::tdual_int_scalar("Dihedral:warning_flag"); d_warning_flag = k_warning_flag.view(); h_warning_flag = k_warning_flag.h_view; + + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/dihedral_harmonic_kokkos.cpp b/src/KOKKOS/dihedral_harmonic_kokkos.cpp index 2681e84e94..764981f231 100644 --- a/src/KOKKOS/dihedral_harmonic_kokkos.cpp +++ b/src/KOKKOS/dihedral_harmonic_kokkos.cpp @@ -47,6 +47,8 @@ DihedralHarmonicKokkos::DihedralHarmonicKokkos(LAMMPS *lmp) : Dihedr k_warning_flag = DAT::tdual_int_scalar("Dihedral:warning_flag"); d_warning_flag = k_warning_flag.view(); h_warning_flag = k_warning_flag.h_view; + + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/dihedral_opls_kokkos.cpp b/src/KOKKOS/dihedral_opls_kokkos.cpp index 01c7f89e0c..ae2b522f81 100644 --- a/src/KOKKOS/dihedral_opls_kokkos.cpp +++ b/src/KOKKOS/dihedral_opls_kokkos.cpp @@ -47,6 +47,8 @@ DihedralOPLSKokkos::DihedralOPLSKokkos(LAMMPS *lmp) : DihedralOPLS(l k_warning_flag = DAT::tdual_int_scalar("Dihedral:warning_flag"); d_warning_flag = k_warning_flag.view(); h_warning_flag = k_warning_flag.h_view; + + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/improper_class2_kokkos.cpp b/src/KOKKOS/improper_class2_kokkos.cpp index f0b27afb3f..68f7223c44 100644 --- a/src/KOKKOS/improper_class2_kokkos.cpp +++ b/src/KOKKOS/improper_class2_kokkos.cpp @@ -43,6 +43,8 @@ ImproperClass2Kokkos::ImproperClass2Kokkos(LAMMPS *lmp) : ImproperCl k_warning_flag = DAT::tdual_int_scalar("Dihedral:warning_flag"); d_warning_flag = k_warning_flag.view(); h_warning_flag = k_warning_flag.h_view; + + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/improper_harmonic_kokkos.cpp b/src/KOKKOS/improper_harmonic_kokkos.cpp index 444088c156..8a64e679bf 100644 --- a/src/KOKKOS/improper_harmonic_kokkos.cpp +++ b/src/KOKKOS/improper_harmonic_kokkos.cpp @@ -44,6 +44,8 @@ ImproperHarmonicKokkos::ImproperHarmonicKokkos(LAMMPS *lmp) : Improp k_warning_flag = Kokkos::DualView("Dihedral:warning_flag"); d_warning_flag = k_warning_flag.template view(); h_warning_flag = k_warning_flag.h_view; + + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ diff --git a/src/KOKKOS/pair_hybrid_kokkos.cpp b/src/KOKKOS/pair_hybrid_kokkos.cpp index d4888883b0..0d28e4a15d 100644 --- a/src/KOKKOS/pair_hybrid_kokkos.cpp +++ b/src/KOKKOS/pair_hybrid_kokkos.cpp @@ -64,23 +64,23 @@ void PairHybridKokkos::compute(int eflag, int vflag) int i,j,m,n; // if no_virial_fdotr_compute is set and global component of - // incoming vflag = 2, then - // reset vflag as if global component were 1 + // incoming vflag = VIRIAL_FDOTR, then + // reset vflag as if global component were VIRIAL_PAIR // necessary since one or more sub-styles cannot compute virial as F dot r int neighflag = lmp->kokkos->neighflag; if (neighflag == FULL) no_virial_fdotr_compute = 1; - if (no_virial_fdotr_compute && vflag % 4 == 2) vflag = 1 + vflag/4 * 4; + if (no_virial_fdotr_compute && vflag & VIRIAL_FDOTR) vflag = VIRIAL_PAIR | (vflag & ~VIRIAL_FDOTR); ev_init(eflag,vflag); - // check if global component of incoming vflag = 2 - // if so, reset vflag passed to substyle as if it were 0 + // check if global component of incoming vflag = VIRIAL_FDOTR + // if so, reset vflag passed to substyle as if it were VIRIAL_NONE // necessary so substyle will not invoke virial_fdotr_compute() int vflag_substyle; - if (vflag % 4 == 2) vflag_substyle = vflag/4 * 4; + if (vflag & VIRIAL_FDOTR) vflag_substyle = VIRIAL_NONE | (vflag & ~VIRIAL_FDOTR); else vflag_substyle = vflag; double *saved_special = save_special(); diff --git a/src/KOKKOS/pair_reaxc_kokkos.cpp b/src/KOKKOS/pair_reaxc_kokkos.cpp index 12d4c3924a..eea397dc69 100644 --- a/src/KOKKOS/pair_reaxc_kokkos.cpp +++ b/src/KOKKOS/pair_reaxc_kokkos.cpp @@ -3629,7 +3629,9 @@ void PairReaxCKokkos::v_tally3_atom(EV_FLOAT_REAX &ev, const int &i, /* ---------------------------------------------------------------------- setup for energy, virial computation - see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) + see integrate::ev_set() for values of eflag and vflag + see pair::ev_setup() for values of eflag_* and vflag_* + VIRIAL_CENTROID bitflag is not yet supported by ReaxFF ------------------------------------------------------------------------- */ template @@ -3640,12 +3642,12 @@ void PairReaxCKokkos::ev_setup(int eflag, int vflag, int) evflag = 1; eflag_either = eflag; - eflag_global = eflag % 2; - eflag_atom = eflag / 2; + eflag_global = eflag & ENERGY_GLOBAL; + eflag_atom = eflag & ENERGY_ATOM; vflag_either = vflag; - vflag_global = vflag % 4; - vflag_atom = vflag / 4; + vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); + vflag_atom = vflag & VIRIAL_ATOM; // reallocate per-atom arrays if necessary @@ -3673,11 +3675,11 @@ void PairReaxCKokkos::ev_setup(int eflag, int vflag, int) Kokkos::parallel_for(Kokkos::RangePolicy(0,maxvatom),*this); } - // if vflag_global = 2 and pair::compute() calls virial_fdotr_compute() + // if vflag_global = VIRIAL_FDOTR and pair::compute() calls virial_fdotr_compute() // compute global virial via (F dot r) instead of via pairwise summation // unset other flags as appropriate - if (vflag_global == 2 && no_virial_fdotr_compute == 0) { + if (vflag_global == VIRIAL_FDOTR && no_virial_fdotr_compute == 0) { vflag_fdotr = 1; vflag_global = 0; if (vflag_atom == 0) vflag_either = 0; diff --git a/src/MANYBODY/pair_adp.cpp b/src/MANYBODY/pair_adp.cpp index fffa53b3d8..401c1f35ad 100644 --- a/src/MANYBODY/pair_adp.cpp +++ b/src/MANYBODY/pair_adp.cpp @@ -69,6 +69,7 @@ PairADP::PairADP(LAMMPS *lmp) : Pair(lmp) single_enable = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- diff --git a/src/MANYBODY/pair_airebo.cpp b/src/MANYBODY/pair_airebo.cpp index 5955e8c29c..4a0c2a874d 100644 --- a/src/MANYBODY/pair_airebo.cpp +++ b/src/MANYBODY/pair_airebo.cpp @@ -69,6 +69,7 @@ PairAIREBO::PairAIREBO(LAMMPS *lmp) nC = nH = nullptr; map = nullptr; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; sigwid = 0.84; sigcut = 3.0; diff --git a/src/MANYBODY/pair_atm.cpp b/src/MANYBODY/pair_atm.cpp index 5b15909845..460e1f6828 100644 --- a/src/MANYBODY/pair_atm.cpp +++ b/src/MANYBODY/pair_atm.cpp @@ -52,6 +52,7 @@ PairATM::PairATM(LAMMPS *lmp) : Pair(lmp) restartinfo = 1; one_coeff = 0; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- diff --git a/src/MANYBODY/pair_bop.cpp b/src/MANYBODY/pair_bop.cpp index fb2b06f48b..29fd8b63c6 100644 --- a/src/MANYBODY/pair_bop.cpp +++ b/src/MANYBODY/pair_bop.cpp @@ -61,6 +61,7 @@ PairBOP::PairBOP(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; ghostneigh = 1; allocated = 0; diff --git a/src/MANYBODY/pair_comb.cpp b/src/MANYBODY/pair_comb.cpp index 916b00d8d4..b850e9d2c1 100644 --- a/src/MANYBODY/pair_comb.cpp +++ b/src/MANYBODY/pair_comb.cpp @@ -53,6 +53,7 @@ PairComb::PairComb(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; nmax = 0; NCo = nullptr; diff --git a/src/MANYBODY/pair_comb3.cpp b/src/MANYBODY/pair_comb3.cpp index 4465cfdcb0..ddf4695df7 100644 --- a/src/MANYBODY/pair_comb3.cpp +++ b/src/MANYBODY/pair_comb3.cpp @@ -51,6 +51,8 @@ PairComb3::PairComb3(LAMMPS *lmp) : Pair(lmp) single_enable = 0; restartinfo = 0; one_coeff = 1; + manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; ghostneigh = 1; nmax = 0; diff --git a/src/MANYBODY/pair_eim.cpp b/src/MANYBODY/pair_eim.cpp index 5606dc137b..dbb16125c7 100644 --- a/src/MANYBODY/pair_eim.cpp +++ b/src/MANYBODY/pair_eim.cpp @@ -39,6 +39,7 @@ PairEIM::PairEIM(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); setfl = nullptr; diff --git a/src/MANYBODY/pair_gw.cpp b/src/MANYBODY/pair_gw.cpp index 4a67e51b79..675ec1b3fe 100644 --- a/src/MANYBODY/pair_gw.cpp +++ b/src/MANYBODY/pair_gw.cpp @@ -48,6 +48,7 @@ PairGW::PairGW(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; diff --git a/src/MANYBODY/pair_lcbop.cpp b/src/MANYBODY/pair_lcbop.cpp index 6184ebfe29..df43aac9d2 100644 --- a/src/MANYBODY/pair_lcbop.cpp +++ b/src/MANYBODY/pair_lcbop.cpp @@ -45,6 +45,7 @@ PairLCBOP::PairLCBOP(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; ghostneigh = 1; maxlocal = 0; diff --git a/src/MANYBODY/pair_nb3b_harmonic.cpp b/src/MANYBODY/pair_nb3b_harmonic.cpp index 369d3bfc04..ffb059c067 100644 --- a/src/MANYBODY/pair_nb3b_harmonic.cpp +++ b/src/MANYBODY/pair_nb3b_harmonic.cpp @@ -46,7 +46,8 @@ PairNb3bHarmonic::PairNb3bHarmonic(LAMMPS *lmp) : Pair(lmp) single_enable = 0; restartinfo = 0; one_coeff = 1; - manybody_flag = 1; + manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; diff --git a/src/MANYBODY/pair_polymorphic.cpp b/src/MANYBODY/pair_polymorphic.cpp index b045ec66d1..e256063b73 100644 --- a/src/MANYBODY/pair_polymorphic.cpp +++ b/src/MANYBODY/pair_polymorphic.cpp @@ -44,6 +44,8 @@ PairPolymorphic::PairPolymorphic(LAMMPS *lmp) : Pair(lmp) single_enable = 0; restartinfo = 0; one_coeff = 1; + manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; nelements = 0; elements = nullptr; diff --git a/src/MANYBODY/pair_sw.cpp b/src/MANYBODY/pair_sw.cpp index a300c11e1a..708027913c 100644 --- a/src/MANYBODY/pair_sw.cpp +++ b/src/MANYBODY/pair_sw.cpp @@ -43,6 +43,7 @@ PairSW::PairSW(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; diff --git a/src/MANYBODY/pair_tersoff.cpp b/src/MANYBODY/pair_tersoff.cpp index be285d7268..dbb7cd2f48 100644 --- a/src/MANYBODY/pair_tersoff.cpp +++ b/src/MANYBODY/pair_tersoff.cpp @@ -47,6 +47,7 @@ PairTersoff::PairTersoff(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; diff --git a/src/MANYBODY/pair_vashishta.cpp b/src/MANYBODY/pair_vashishta.cpp index 550fd661d6..d5a1374a80 100644 --- a/src/MANYBODY/pair_vashishta.cpp +++ b/src/MANYBODY/pair_vashishta.cpp @@ -45,6 +45,7 @@ PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; diff --git a/src/MLIAP/pair_mliap.cpp b/src/MLIAP/pair_mliap.cpp index f5895d6e5d..a422271f19 100644 --- a/src/MLIAP/pair_mliap.cpp +++ b/src/MLIAP/pair_mliap.cpp @@ -38,7 +38,7 @@ PairMLIAP::PairMLIAP(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; - + centroidstressflag = CENTROID_NOTAVAIL; } /* ---------------------------------------------------------------------- */ diff --git a/src/MOLECULE/bond_quartic.cpp b/src/MOLECULE/bond_quartic.cpp index 20269a0966..c7151c3aca 100644 --- a/src/MOLECULE/bond_quartic.cpp +++ b/src/MOLECULE/bond_quartic.cpp @@ -63,7 +63,7 @@ void BondQuartic::compute(int eflag, int vflag) // insure pair->ev_tally() will use 1-4 virial contribution - if (vflag_global == 2) + if (vflag_global == VIRIAL_FDOTR) force->pair->vflag_either = force->pair->vflag_global = 1; double **cutsq = force->pair->cutsq; diff --git a/src/MOLECULE/dihedral_charmm.cpp b/src/MOLECULE/dihedral_charmm.cpp index c3d94976ac..cc60dbec2b 100644 --- a/src/MOLECULE/dihedral_charmm.cpp +++ b/src/MOLECULE/dihedral_charmm.cpp @@ -79,7 +79,7 @@ void DihedralCharmm::compute(int eflag, int vflag) // insure pair->ev_tally() will use 1-4 virial contribution - if (weightflag && vflag_global == 2) + if (weightflag && vflag_global == VIRIAL_FDOTR) force->pair->vflag_either = force->pair->vflag_global = 1; double **x = atom->x; diff --git a/src/MOLECULE/dihedral_charmmfsw.cpp b/src/MOLECULE/dihedral_charmmfsw.cpp index 1f86014933..1164e93f18 100644 --- a/src/MOLECULE/dihedral_charmmfsw.cpp +++ b/src/MOLECULE/dihedral_charmmfsw.cpp @@ -82,7 +82,7 @@ void DihedralCharmmfsw::compute(int eflag, int vflag) // insure pair->ev_tally() will use 1-4 virial contribution - if (weightflag && vflag_global == 2) + if (weightflag && vflag_global == VIRIAL_FDOTR) force->pair->vflag_either = force->pair->vflag_global = 1; double **x = atom->x; diff --git a/src/MOLECULE/fix_cmap.cpp b/src/MOLECULE/fix_cmap.cpp index ab5c0a7323..f48d5254cf 100644 --- a/src/MOLECULE/fix_cmap.cpp +++ b/src/MOLECULE/fix_cmap.cpp @@ -42,8 +42,6 @@ #include "memory.h" #include "error.h" - - using namespace LAMMPS_NS; using namespace FixConst; using namespace MathConst; @@ -74,6 +72,7 @@ FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) : restart_peratom = 1; peatom_flag = 1; virial_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; thermo_virial = 1; peratom_freq = 1; scalar_flag = 1; diff --git a/src/POEMS/fix_poems.cpp b/src/POEMS/fix_poems.cpp index fb39c21f2a..a220216fdd 100644 --- a/src/POEMS/fix_poems.cpp +++ b/src/POEMS/fix_poems.cpp @@ -74,6 +74,7 @@ FixPOEMS::FixPOEMS(LAMMPS *lmp, int narg, char **arg) : time_integrate = 1; rigid_flag = 1; virial_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; thermo_virial = 1; dof_flag = 1; diff --git a/src/PYTHON/pair_python.cpp b/src/PYTHON/pair_python.cpp index ab726c8a28..eb54d1d4e5 100644 --- a/src/PYTHON/pair_python.cpp +++ b/src/PYTHON/pair_python.cpp @@ -41,7 +41,7 @@ PairPython::PairPython(LAMMPS *lmp) : Pair(lmp) { one_coeff = 1; reinitflag = 0; cut_global = 0.0; - centroidstressflag = 1; + centroidstressflag = CENTROID_SAME; py_potential = nullptr; skip_types = nullptr; diff --git a/src/SNAP/pair_snap.cpp b/src/SNAP/pair_snap.cpp index e6ff85b4b6..a7ea6de591 100644 --- a/src/SNAP/pair_snap.cpp +++ b/src/SNAP/pair_snap.cpp @@ -40,6 +40,7 @@ PairSNAP::PairSNAP(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; nelements = 0; elements = nullptr; diff --git a/src/USER-FEP/pair_coul_cut_soft.cpp b/src/USER-FEP/pair_coul_cut_soft.cpp index 3540fb0504..0928e18d6a 100644 --- a/src/USER-FEP/pair_coul_cut_soft.cpp +++ b/src/USER-FEP/pair_coul_cut_soft.cpp @@ -33,7 +33,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ PairCoulCutSoft::PairCoulCutSoft(LAMMPS *lmp) : Pair(lmp) { - centroidstressflag = 1; + centroidstressflag = CENTROID_SAME; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-FEP/pair_lj_class2_coul_cut_soft.cpp b/src/USER-FEP/pair_lj_class2_coul_cut_soft.cpp index 7e9f210905..8b462b0ac7 100644 --- a/src/USER-FEP/pair_lj_class2_coul_cut_soft.cpp +++ b/src/USER-FEP/pair_lj_class2_coul_cut_soft.cpp @@ -33,7 +33,7 @@ using namespace MathConst; PairLJClass2CoulCutSoft::PairLJClass2CoulCutSoft(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; + centroidstressflag = CENTROID_SAME; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-FEP/pair_lj_class2_soft.cpp b/src/USER-FEP/pair_lj_class2_soft.cpp index 716bf890cf..df8d1aaab6 100644 --- a/src/USER-FEP/pair_lj_class2_soft.cpp +++ b/src/USER-FEP/pair_lj_class2_soft.cpp @@ -32,7 +32,7 @@ using namespace MathConst; PairLJClass2Soft::PairLJClass2Soft(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; + centroidstressflag = CENTROID_SAME; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp b/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp index f776944ee3..f69ea9a5d9 100644 --- a/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp +++ b/src/USER-FEP/pair_lj_cut_coul_cut_soft.cpp @@ -37,7 +37,7 @@ using namespace MathConst; PairLJCutCoulCutSoft::PairLJCutCoulCutSoft(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; + centroidstressflag = CENTROID_SAME; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-FEP/pair_lj_cut_soft.cpp b/src/USER-FEP/pair_lj_cut_soft.cpp index 6d292acbb3..c3574746e2 100644 --- a/src/USER-FEP/pair_lj_cut_soft.cpp +++ b/src/USER-FEP/pair_lj_cut_soft.cpp @@ -43,7 +43,7 @@ PairLJCutSoft::PairLJCutSoft(LAMMPS *lmp) : Pair(lmp) respa_enable = 1; writedata = 1; allocated = 0; - centroidstressflag = 1; + centroidstressflag = CENTROID_SAME; } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-INTEL/dihedral_charmm_intel.cpp b/src/USER-INTEL/dihedral_charmm_intel.cpp index fb9e8781b1..d84658c797 100644 --- a/src/USER-INTEL/dihedral_charmm_intel.cpp +++ b/src/USER-INTEL/dihedral_charmm_intel.cpp @@ -91,7 +91,7 @@ void DihedralCharmmIntel::compute(int eflag, int vflag, // insure pair->ev_tally() will use 1-4 virial contribution - if (weightflag && vflag_global == 2) + if (weightflag && vflag_global == VIRIAL_FDOTR) force->pair->vflag_either = force->pair->vflag_global = 1; if (evflag) { diff --git a/src/USER-INTEL/intel_preprocess.h b/src/USER-INTEL/intel_preprocess.h index 3e547b58a0..3c285871cf 100644 --- a/src/USER-INTEL/intel_preprocess.h +++ b/src/USER-INTEL/intel_preprocess.h @@ -502,7 +502,7 @@ enum {TIME_PACK, TIME_HOST_NEIGHBOR, TIME_HOST_PAIR, TIME_OFFLOAD_NEIGHBOR, acc_t *f_scalar = &f_start[0].x; \ int f_stride4 = f_stride * 4; \ int t; \ - if (vflag == 2) t = 4; else t = 1; \ + if (vflag == VIRIAL_FDOTR) t = 4; else t = 1; \ acc_t *f_scalar2 = f_scalar + f_stride4 * t; \ for ( ; t < nthreads; t++) { \ _use_simd_pragma("vector aligned") \ @@ -512,7 +512,7 @@ enum {TIME_PACK, TIME_HOST_NEIGHBOR, TIME_HOST_PAIR, TIME_OFFLOAD_NEIGHBOR, f_scalar2 += f_stride4; \ } \ \ - if (vflag == 2) { \ + if (vflag == VIRIAL_FDOTR) { \ int nt_min = MIN(4,nthreads); \ IP_PRE_fdotr_acc_force_l5(iifrom, iito, minlocal, nt_min, f_start, \ f_stride, pos, ov0, ov1, ov2, ov3, ov4, \ diff --git a/src/USER-INTEL/pair_buck_coul_cut_intel.cpp b/src/USER-INTEL/pair_buck_coul_cut_intel.cpp index 80993f212e..3da131684c 100644 --- a/src/USER-INTEL/pair_buck_coul_cut_intel.cpp +++ b/src/USER-INTEL/pair_buck_coul_cut_intel.cpp @@ -245,7 +245,7 @@ void PairBuckCoulCutIntel::eval(const int offload, const int vflag, fxtmp = fytmp = fztmp = (acc_t)0; if (EFLAG) fwtmp = sevdwl = secoul = (acc_t)0; if (NEWTON_PAIR == 0) - if (vflag==1) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; + if (vflag == VIRIAL_PAIR) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; #if defined(LMP_SIMD_COMPILER) #pragma vector aligned diff --git a/src/USER-INTEL/pair_buck_intel.cpp b/src/USER-INTEL/pair_buck_intel.cpp index 95b37796ba..afa9b448b5 100644 --- a/src/USER-INTEL/pair_buck_intel.cpp +++ b/src/USER-INTEL/pair_buck_intel.cpp @@ -227,7 +227,7 @@ void PairBuckIntel::eval(const int offload, const int vflag, fxtmp = fytmp = fztmp = (acc_t)0; if (EFLAG) fwtmp = sevdwl = (acc_t)0; if (NEWTON_PAIR == 0) - if (vflag==1) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; + if (vflag == VIRIAL_PAIR) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; #if defined(LMP_SIMD_COMPILER) #pragma vector aligned diff --git a/src/USER-INTEL/pair_dpd_intel.cpp b/src/USER-INTEL/pair_dpd_intel.cpp index 5adcb4c6ea..a0a71e4a37 100644 --- a/src/USER-INTEL/pair_dpd_intel.cpp +++ b/src/USER-INTEL/pair_dpd_intel.cpp @@ -265,7 +265,7 @@ void PairDPDIntel::eval(const int offload, const int vflag, fxtmp = fytmp = fztmp = (acc_t)0; if (EFLAG) fwtmp = sevdwl = (acc_t)0; if (NEWTON_PAIR == 0) - if (vflag==1) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; + if (vflag == VIRIAL_PAIR) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; if (rngi + jnum > rng_size) { #ifdef LMP_USE_MKL_RNG diff --git a/src/USER-INTEL/pair_eam_intel.cpp b/src/USER-INTEL/pair_eam_intel.cpp index 6ebd7abe3d..a25d588ecf 100644 --- a/src/USER-INTEL/pair_eam_intel.cpp +++ b/src/USER-INTEL/pair_eam_intel.cpp @@ -481,7 +481,7 @@ void PairEAMIntel::eval(const int offload, const int vflag, fxtmp = fytmp = fztmp = (acc_t)0; if (EFLAG) fwtmp = sevdwl = (acc_t)0; if (NEWTON_PAIR == 0) - if (vflag==1) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; + if (vflag == VIRIAL_PAIR) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; int ej = 0; #if defined(LMP_SIMD_COMPILER) diff --git a/src/USER-INTEL/pair_gayberne_intel.cpp b/src/USER-INTEL/pair_gayberne_intel.cpp index a9aa00340e..c34c6965de 100644 --- a/src/USER-INTEL/pair_gayberne_intel.cpp +++ b/src/USER-INTEL/pair_gayberne_intel.cpp @@ -399,7 +399,7 @@ void PairGayBerneIntel::eval(const int offload, const int vflag, if (EFLAG) fwtmp = sevdwl = (acc_t)0.0; if (NEWTON_PAIR == 0) - if (vflag==1) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0.0; + if (vflag == VIRIAL_PAIR) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0.0; bool multiple_forms = false; int packed_j = 0; diff --git a/src/USER-INTEL/pair_lj_charmm_coul_charmm_intel.cpp b/src/USER-INTEL/pair_lj_charmm_coul_charmm_intel.cpp index bf55a53981..8f4dba4a24 100644 --- a/src/USER-INTEL/pair_lj_charmm_coul_charmm_intel.cpp +++ b/src/USER-INTEL/pair_lj_charmm_coul_charmm_intel.cpp @@ -265,7 +265,7 @@ void PairLJCharmmCoulCharmmIntel::eval(const int offload, const int vflag, fxtmp = fytmp = fztmp = (acc_t)0; if (EFLAG) fwtmp = sevdwl = secoul = (acc_t)0; if (NEWTON_PAIR == 0) - if (vflag==1) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; + if (vflag == VIRIAL_PAIR) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; int ej = 0; #if defined(LMP_SIMD_COMPILER) diff --git a/src/USER-INTEL/pair_lj_charmm_coul_long_intel.cpp b/src/USER-INTEL/pair_lj_charmm_coul_long_intel.cpp index 2f375fd08c..ffed4853ba 100644 --- a/src/USER-INTEL/pair_lj_charmm_coul_long_intel.cpp +++ b/src/USER-INTEL/pair_lj_charmm_coul_long_intel.cpp @@ -282,7 +282,7 @@ void PairLJCharmmCoulLongIntel::eval(const int offload, const int vflag, fxtmp = fytmp = fztmp = (acc_t)0; if (EFLAG) fwtmp = sevdwl = secoul = (acc_t)0; if (NEWTON_PAIR == 0) - if (vflag==1) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; + if (vflag == VIRIAL_PAIR) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; int ej = 0; #if defined(LMP_SIMD_COMPILER) diff --git a/src/USER-INTEL/pair_lj_cut_coul_long_intel.cpp b/src/USER-INTEL/pair_lj_cut_coul_long_intel.cpp index fa3cfd7bc3..c1d85ad37f 100644 --- a/src/USER-INTEL/pair_lj_cut_coul_long_intel.cpp +++ b/src/USER-INTEL/pair_lj_cut_coul_long_intel.cpp @@ -275,7 +275,7 @@ void PairLJCutCoulLongIntel::eval(const int offload, const int vflag, fxtmp = fytmp = fztmp = (acc_t)0; if (EFLAG) fwtmp = sevdwl = secoul = (acc_t)0; if (NEWTON_PAIR == 0) - if (vflag==1) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; + if (vflag == VIRIAL_PAIR) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; int ej = 0; #if defined(LMP_SIMD_COMPILER) diff --git a/src/USER-INTEL/pair_lj_cut_intel.cpp b/src/USER-INTEL/pair_lj_cut_intel.cpp index 426abb7660..8697a4f548 100644 --- a/src/USER-INTEL/pair_lj_cut_intel.cpp +++ b/src/USER-INTEL/pair_lj_cut_intel.cpp @@ -233,7 +233,7 @@ void PairLJCutIntel::eval(const int offload, const int vflag, fxtmp = fytmp = fztmp = (acc_t)0; if (EFLAG) fwtmp = sevdwl = (acc_t)0; if (NEWTON_PAIR == 0) - if (vflag==1) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; + if (vflag == VIRIAL_PAIR) sv0 = sv1 = sv2 = sv3 = sv4 = sv5 = (acc_t)0; #if defined(LMP_SIMD_COMPILER) #pragma vector aligned nog2s diff --git a/src/USER-MEAMC/pair_meamc.cpp b/src/USER-MEAMC/pair_meamc.cpp index e99c17cb31..d0fb8a0463 100644 --- a/src/USER-MEAMC/pair_meamc.cpp +++ b/src/USER-MEAMC/pair_meamc.cpp @@ -17,7 +17,6 @@ #include "pair_meamc.h" - #include #include "meam.h" @@ -30,8 +29,6 @@ #include "memory.h" #include "error.h" - - using namespace LAMMPS_NS; #define MAXLINE 1024 @@ -52,6 +49,7 @@ PairMEAMC::PairMEAMC(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; allocated = 0; diff --git a/src/USER-MESONT/pair_mesocnt.cpp b/src/USER-MESONT/pair_mesocnt.cpp index b4e21888e0..39a97321a7 100644 --- a/src/USER-MESONT/pair_mesocnt.cpp +++ b/src/USER-MESONT/pair_mesocnt.cpp @@ -23,7 +23,6 @@ #include - #include "atom.h" #include "comm.h" #include "force.h" @@ -33,9 +32,6 @@ #include "memory.h" #include "error.h" #include "update.h" - - - #include "math_const.h" #include "math_extra.h" @@ -58,6 +54,7 @@ PairMesoCNT::PairMesoCNT(LAMMPS *lmp) : Pair(lmp) respa_enable = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; no_virial_fdotr_compute = 0; writedata = 0; ghostneigh = 0; diff --git a/src/USER-MISC/pair_agni.cpp b/src/USER-MISC/pair_agni.cpp index 52f1db923c..afb457d9e4 100644 --- a/src/USER-MISC/pair_agni.cpp +++ b/src/USER-MISC/pair_agni.cpp @@ -20,6 +20,7 @@ #include "atom.h" #include "citeme.h" #include "comm.h" +#include "force.h" #include "error.h" #include "math_const.h" #include "math_special.h" @@ -81,6 +82,7 @@ PairAGNI::PairAGNI(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; no_virial_fdotr_compute = 1; diff --git a/src/USER-MISC/pair_drip.cpp b/src/USER-MISC/pair_drip.cpp index a58b77782a..45b49ee293 100644 --- a/src/USER-MISC/pair_drip.cpp +++ b/src/USER-MISC/pair_drip.cpp @@ -34,7 +34,6 @@ #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; #define MAXLINE 1024 @@ -48,6 +47,7 @@ PairDRIP::PairDRIP(LAMMPS *lmp) : Pair(lmp) single_enable = 0; restartinfo = 0; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; params = nullptr; nearest3neigh = nullptr; diff --git a/src/USER-MISC/pair_edip.cpp b/src/USER-MISC/pair_edip.cpp index 52fcfa85fb..b29f499de8 100644 --- a/src/USER-MISC/pair_edip.cpp +++ b/src/USER-MISC/pair_edip.cpp @@ -36,7 +36,6 @@ #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; #define MAXLINE 1024 @@ -63,6 +62,7 @@ PairEDIP::PairEDIP(LAMMPS *lmp) : restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; nelements = 0; elements = nullptr; diff --git a/src/USER-MISC/pair_edip_multi.cpp b/src/USER-MISC/pair_edip_multi.cpp index 7213ca37ed..811f80abfc 100644 --- a/src/USER-MISC/pair_edip_multi.cpp +++ b/src/USER-MISC/pair_edip_multi.cpp @@ -68,6 +68,7 @@ PairEDIPMulti::PairEDIPMulti(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; nelements = 0; elements = nullptr; diff --git a/src/USER-MISC/pair_extep.cpp b/src/USER-MISC/pair_extep.cpp index 2bf5aee947..2d0001c680 100644 --- a/src/USER-MISC/pair_extep.cpp +++ b/src/USER-MISC/pair_extep.cpp @@ -49,6 +49,7 @@ PairExTeP::PairExTeP(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; ghostneigh = 1; nelements = 0; diff --git a/src/USER-MISC/pair_meam_sw_spline.cpp b/src/USER-MISC/pair_meam_sw_spline.cpp index c27cdcf610..17511e90b8 100644 --- a/src/USER-MISC/pair_meam_sw_spline.cpp +++ b/src/USER-MISC/pair_meam_sw_spline.cpp @@ -47,6 +47,7 @@ PairMEAMSWSpline::PairMEAMSWSpline(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; nelements = 0; elements = nullptr; diff --git a/src/USER-MISC/pair_tersoff_table.cpp b/src/USER-MISC/pair_tersoff_table.cpp index feaa789688..c4e2d2c993 100644 --- a/src/USER-MISC/pair_tersoff_table.cpp +++ b/src/USER-MISC/pair_tersoff_table.cpp @@ -61,6 +61,7 @@ PairTersoffTable::PairTersoffTable(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; unit_convert_flag = utils::get_supported_conversions(utils::ENERGY); nelements = 0; diff --git a/src/USER-MOFFF/angle_cosine_buck6d.cpp b/src/USER-MOFFF/angle_cosine_buck6d.cpp index 3ab7327325..277ee8830b 100644 --- a/src/USER-MOFFF/angle_cosine_buck6d.cpp +++ b/src/USER-MOFFF/angle_cosine_buck6d.cpp @@ -72,7 +72,7 @@ void AngleCosineBuck6d::compute(int eflag, int vflag) // insure pair->ev_tally() will use 1-3 virial contribution - if (vflag_global == 2) + if (vflag_global == VIRIAL_FDOTR) force->pair->vflag_either = force->pair->vflag_global = 1; double **x = atom->x; diff --git a/src/USER-OMP/bond_quartic_omp.cpp b/src/USER-OMP/bond_quartic_omp.cpp index 73fa0614d8..1f39a7cee1 100644 --- a/src/USER-OMP/bond_quartic_omp.cpp +++ b/src/USER-OMP/bond_quartic_omp.cpp @@ -45,7 +45,7 @@ void BondQuarticOMP::compute(int eflag, int vflag) // insure pair->ev_tally() will use 1-4 virial contribution - if (vflag_global == 2) + if (vflag_global == VIRIAL_FDOTR) force->pair->vflag_either = force->pair->vflag_global = 1; const int nall = atom->nlocal + atom->nghost; diff --git a/src/USER-OMP/dihedral_charmm_omp.cpp b/src/USER-OMP/dihedral_charmm_omp.cpp index 6b399be09d..c75378a8be 100644 --- a/src/USER-OMP/dihedral_charmm_omp.cpp +++ b/src/USER-OMP/dihedral_charmm_omp.cpp @@ -49,7 +49,7 @@ void DihedralCharmmOMP::compute(int eflag, int vflag) // insure pair->ev_tally() will use 1-4 virial contribution - if (weightflag && vflag_global == 2) + if (weightflag && vflag_global == VIRIAL_FDOTR) force->pair->vflag_either = force->pair->vflag_global = 1; const int nall = atom->nlocal + atom->nghost; diff --git a/src/USER-OMP/thr_omp.cpp b/src/USER-OMP/thr_omp.cpp index 0a8611707b..e590aedc87 100644 --- a/src/USER-OMP/thr_omp.cpp +++ b/src/USER-OMP/thr_omp.cpp @@ -70,20 +70,20 @@ void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom, if (tid == 0) thr_error = 0; if (thr_style & THR_PAIR) { - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { thr->eatom_pair = eatom + tid*nall; if (nall > 0) memset(&(thr->eatom_pair[0]),0,nall*sizeof(double)); } // per-atom virial and per-atom centroid virial are the same for two-body // many-body pair styles not yet implemented - if (vflag & 12) { + if (vflag & (VIRIAL_ATOM | VIRIAL_CENTROID)) { thr->vatom_pair = vatom + tid*nall; if (nall > 0) memset(&(thr->vatom_pair[0][0]),0,nall*6*sizeof(double)); } // check cvatom_pair, because can't access centroidstressflag - if ((vflag & 8) && cvatom) { + if ((vflag & VIRIAL_CENTROID) && cvatom) { thr->cvatom_pair = cvatom + tid*nall; if (nall > 0) memset(&(thr->cvatom_pair[0][0]),0,nall*9*sizeof(double)); @@ -94,13 +94,13 @@ void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom, } if (thr_style & THR_BOND) { - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { thr->eatom_bond = eatom + tid*nall; if (nall > 0) memset(&(thr->eatom_bond[0]),0,nall*sizeof(double)); } // per-atom virial and per-atom centroid virial are the same for bonds - if (vflag & 12) { + if (vflag & (VIRIAL_ATOM | VIRIAL_CENTROID)) { thr->vatom_bond = vatom + tid*nall; if (nall > 0) memset(&(thr->vatom_bond[0][0]),0,nall*6*sizeof(double)); @@ -108,17 +108,17 @@ void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom, } if (thr_style & THR_ANGLE) { - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { thr->eatom_angle = eatom + tid*nall; if (nall > 0) memset(&(thr->eatom_angle[0]),0,nall*sizeof(double)); } - if (vflag & 4) { + if (vflag & VIRIAL_ATOM) { thr->vatom_angle = vatom + tid*nall; if (nall > 0) memset(&(thr->vatom_angle[0][0]),0,nall*6*sizeof(double)); } - if (vflag & 8) { + if (vflag & VIRIAL_CENTROID) { thr->cvatom_angle = cvatom + tid*nall; if (nall > 0) memset(&(thr->cvatom_angle[0][0]),0,nall*9*sizeof(double)); @@ -126,17 +126,17 @@ void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom, } if (thr_style & THR_DIHEDRAL) { - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { thr->eatom_dihed = eatom + tid*nall; if (nall > 0) memset(&(thr->eatom_dihed[0]),0,nall*sizeof(double)); } - if (vflag & 4) { + if (vflag & VIRIAL_ATOM) { thr->vatom_dihed = vatom + tid*nall; if (nall > 0) memset(&(thr->vatom_dihed[0][0]),0,nall*6*sizeof(double)); } - if (vflag & 8) { + if (vflag & VIRIAL_CENTROID) { thr->cvatom_dihed = cvatom + tid*nall; if (nall > 0) memset(&(thr->cvatom_dihed[0][0]),0,nall*9*sizeof(double)); @@ -144,17 +144,17 @@ void ThrOMP::ev_setup_thr(int eflag, int vflag, int nall, double *eatom, } if (thr_style & THR_IMPROPER) { - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { thr->eatom_imprp = eatom + tid*nall; if (nall > 0) memset(&(thr->eatom_imprp[0]),0,nall*sizeof(double)); } - if (vflag & 4) { + if (vflag & VIRIAL_ATOM) { thr->vatom_imprp = vatom + tid*nall; if (nall > 0) memset(&(thr->vatom_imprp[0][0]),0,nall*6*sizeof(double)); } - if (vflag & 8) { + if (vflag & VIRIAL_CENTROID) { thr->cvatom_imprp = cvatom + tid*nall; if (nall > 0) memset(&(thr->cvatom_imprp[0][0]),0,nall*9*sizeof(double)); @@ -222,29 +222,29 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, #pragma omp critical #endif { - if (eflag & 1) { + if (eflag & ENERGY_GLOBAL) { pair->eng_vdwl += thr->eng_vdwl; pair->eng_coul += thr->eng_coul; thr->eng_vdwl = 0.0; thr->eng_coul = 0.0; } - if (vflag & 3) + if (vflag & (VIRIAL_PAIR | VIRIAL_FDOTR)) for (int i=0; i < 6; ++i) { pair->virial[i] += thr->virial_pair[i]; thr->virial_pair[i] = 0.0; } } - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { data_reduce_thr(&(pair->eatom[0]), nall, nthreads, 1, tid); } // per-atom virial and per-atom centroid virial are the same for two-body // many-body pair styles not yet implemented - if (vflag & 12) { + if (vflag & (VIRIAL_ATOM | VIRIAL_CENTROID)) { data_reduce_thr(&(pair->vatom[0][0]), nall, nthreads, 6, tid); } // check cvatom_pair, because can't access centroidstressflag - if ((vflag & 8) && thr->cvatom_pair) { + if ((vflag & VIRIAL_CENTROID) && thr->cvatom_pair) { data_reduce_thr(&(pair->cvatom[0][0]), nall, nthreads, 9, tid); } } @@ -259,12 +259,12 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, #pragma omp critical #endif { - if (eflag & 1) { + if (eflag & ENERGY_GLOBAL) { bond->energy += thr->eng_bond; thr->eng_bond = 0.0; } - if (vflag & 3) { + if (vflag & (VIRIAL_PAIR | VIRIAL_FDOTR)) { for (int i=0; i < 6; ++i) { bond->virial[i] += thr->virial_bond[i]; thr->virial_bond[i] = 0.0; @@ -272,11 +272,11 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, } } - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { data_reduce_thr(&(bond->eatom[0]), nall, nthreads, 1, tid); } // per-atom virial and per-atom centroid virial are the same for bonds - if (vflag & 12) { + if (vflag & (VIRIAL_ATOM | VIRIAL_CENTROID)) { data_reduce_thr(&(bond->vatom[0][0]), nall, nthreads, 6, tid); } @@ -291,12 +291,12 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, #pragma omp critical #endif { - if (eflag & 1) { + if (eflag & ENERGY_GLOBAL) { angle->energy += thr->eng_angle; thr->eng_angle = 0.0; } - if (vflag & 3) { + if (vflag & (VIRIAL_PAIR | VIRIAL_FDOTR)) { for (int i=0; i < 6; ++i) { angle->virial[i] += thr->virial_angle[i]; thr->virial_angle[i] = 0.0; @@ -304,13 +304,13 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, } } - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { data_reduce_thr(&(angle->eatom[0]), nall, nthreads, 1, tid); } - if (vflag & 4) { + if (vflag & VIRIAL_ATOM) { data_reduce_thr(&(angle->vatom[0][0]), nall, nthreads, 6, tid); } - if (vflag & 8) { + if (vflag & VIRIAL_CENTROID) { data_reduce_thr(&(angle->cvatom[0][0]), nall, nthreads, 9, tid); } @@ -325,12 +325,12 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, #pragma omp critical #endif { - if (eflag & 1) { + if (eflag & ENERGY_GLOBAL) { dihedral->energy += thr->eng_dihed; thr->eng_dihed = 0.0; } - if (vflag & 3) { + if (vflag & (VIRIAL_PAIR | VIRIAL_FDOTR)) { for (int i=0; i < 6; ++i) { dihedral->virial[i] += thr->virial_dihed[i]; thr->virial_dihed[i] = 0.0; @@ -338,13 +338,13 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, } } - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { data_reduce_thr(&(dihedral->eatom[0]), nall, nthreads, 1, tid); } - if (vflag & 4) { + if (vflag & VIRIAL_ATOM) { data_reduce_thr(&(dihedral->vatom[0][0]), nall, nthreads, 6, tid); } - if (vflag & 8) { + if (vflag & VIRIAL_CENTROID) { data_reduce_thr(&(dihedral->cvatom[0][0]), nall, nthreads, 9, tid); } @@ -360,7 +360,7 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, #pragma omp critical #endif { - if (eflag & 1) { + if (eflag & ENERGY_GLOBAL) { dihedral->energy += thr->eng_dihed; pair->eng_vdwl += thr->eng_vdwl; pair->eng_coul += thr->eng_coul; @@ -369,7 +369,7 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, thr->eng_coul = 0.0; } - if (vflag & 3) { + if (vflag & (VIRIAL_PAIR | VIRIAL_FDOTR)) { for (int i=0; i < 6; ++i) { dihedral->virial[i] += thr->virial_dihed[i]; pair->virial[i] += thr->virial_pair[i]; @@ -379,23 +379,23 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, } } - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { data_reduce_thr(&(dihedral->eatom[0]), nall, nthreads, 1, tid); data_reduce_thr(&(pair->eatom[0]), nall, nthreads, 1, tid); } - if (vflag & 4) { + if (vflag & VIRIAL_ATOM) { data_reduce_thr(&(dihedral->vatom[0][0]), nall, nthreads, 6, tid); } - if (vflag & 8) { + if (vflag & VIRIAL_CENTROID) { data_reduce_thr(&(dihedral->cvatom[0][0]), nall, nthreads, 9, tid); } // per-atom virial and per-atom centroid virial are the same for two-body // many-body pair styles not yet implemented - if (vflag & 12) { + if (vflag & (VIRIAL_ATOM | VIRIAL_CENTROID)) { data_reduce_thr(&(pair->vatom[0][0]), nall, nthreads, 6, tid); } // check cvatom_pair, because can't access centroidstressflag - if ((vflag & 8) && thr->cvatom_pair) { + if ((vflag & VIRIAL_CENTROID) && thr->cvatom_pair) { data_reduce_thr(&(pair->cvatom[0][0]), nall, nthreads, 9, tid); } } @@ -409,12 +409,12 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, #pragma omp critical #endif { - if (eflag & 1) { + if (eflag & ENERGY_GLOBAL) { improper->energy += thr->eng_imprp; thr->eng_imprp = 0.0; } - if (vflag & 3) { + if (vflag & (VIRIAL_PAIR | VIRIAL_FDOTR)) { for (int i=0; i < 6; ++i) { improper->virial[i] += thr->virial_imprp[i]; thr->virial_imprp[i] = 0.0; @@ -422,13 +422,13 @@ void ThrOMP::reduce_thr(void *style, const int eflag, const int vflag, } } - if (eflag & 2) { + if (eflag & ENERGY_ATOM) { data_reduce_thr(&(improper->eatom[0]), nall, nthreads, 1, tid); } - if (vflag & 4) { + if (vflag & VIRIAL_ATOM) { data_reduce_thr(&(improper->vatom[0][0]), nall, nthreads, 6, tid); } - if (vflag & 8) { + if (vflag & VIRIAL_CENTROID) { data_reduce_thr(&(improper->cvatom[0][0]), nall, nthreads, 9, tid); } diff --git a/src/USER-QUIP/pair_quip.cpp b/src/USER-QUIP/pair_quip.cpp index 6b4b4ff576..0cbc71e564 100644 --- a/src/USER-QUIP/pair_quip.cpp +++ b/src/USER-QUIP/pair_quip.cpp @@ -16,7 +16,6 @@ Aidan Thompson (Sandia, athomps@sandia.gov) ------------------------------------------------------------------------- */ - #include #include @@ -44,6 +43,7 @@ PairQUIP::PairQUIP(LAMMPS *lmp) : Pair(lmp) one_coeff = 1; no_virial_fdotr_compute = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; map = nullptr; quip_potential = nullptr; diff --git a/src/USER-REAXC/pair_reaxc.cpp b/src/USER-REAXC/pair_reaxc.cpp index d13a48cb6d..d3507bdf85 100644 --- a/src/USER-REAXC/pair_reaxc.cpp +++ b/src/USER-REAXC/pair_reaxc.cpp @@ -40,7 +40,6 @@ #include "memory.h" #include "error.h" - #include "reaxc_defs.h" #include "reaxc_types.h" #include "reaxc_allocate.h" @@ -77,6 +76,7 @@ PairReaxC::PairReaxC(LAMMPS *lmp) : Pair(lmp) restartinfo = 0; one_coeff = 1; manybody_flag = 1; + centroidstressflag = CENTROID_NOTAVAIL; ghostneigh = 1; fix_id = new char[24]; diff --git a/src/angle.cpp b/src/angle.cpp index 97b2c4c3e1..110719d432 100644 --- a/src/angle.cpp +++ b/src/angle.cpp @@ -40,6 +40,7 @@ Angle::Angle(LAMMPS *lmp) : Pointers(lmp) vatom = nullptr; cvatom = nullptr; setflag = nullptr; + centroidstressflag = CENTROID_AVAIL; execution_space = Host; datamask_read = ALL_MASK; @@ -75,7 +76,19 @@ void Angle::init() /* ---------------------------------------------------------------------- setup for energy, virial computation - see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) + see integrate::ev_set() for bitwise settings of eflag/vflag + set the following flags, values are otherwise set to 0: + evflag != 0 if any bits of eflag or vflag are set + eflag_global != 0 if ENERGY_GLOBAL bit of eflag set + eflag_atom != 0 if ENERGY_ATOM bit of eflag set + eflag_either != 0 if eflag_global or eflag_atom is set + vflag_global != 0 if VIRIAL_PAIR or VIRIAL_FDOTR bit of vflag set + vflag_atom != 0 if VIRIAL_ATOM bit of vflag set + vflag_atom != 0 if VIRIAL_CENTROID bit of vflag set + and centroidstressflag != CENTROID_AVAIL + cvflag_atom != 0 if VIRIAL_CENTROID bit of vflag set + and centroidstressflag = CENTROID_AVAIL + vflag_either != 0 if any of vflag_global, vflag_atom, cvflag_atom is set ------------------------------------------------------------------------- */ void Angle::ev_setup(int eflag, int vflag, int alloc) @@ -85,13 +98,17 @@ void Angle::ev_setup(int eflag, int vflag, int alloc) evflag = 1; eflag_either = eflag; - eflag_global = eflag % 2; - eflag_atom = eflag / 2; + eflag_global = eflag & ENERGY_GLOBAL; + eflag_atom = eflag & ENERGY_ATOM; - vflag_global = vflag % 4; - vflag_atom = vflag & 4; - cvflag_atom = vflag & 8; - vflag_either = vflag_global || vflag_atom; + vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); + vflag_atom = vflag & VIRIAL_ATOM; + if (vflag & VIRIAL_CENTROID && centroidstressflag != CENTROID_AVAIL) + vflag_atom = 1; + cvflag_atom = 0; + if (vflag & VIRIAL_CENTROID && centroidstressflag == CENTROID_AVAIL) + cvflag_atom = 1; + vflag_either = vflag_global || vflag_atom || cvflag_atom; // reallocate per-atom arrays if necessary diff --git a/src/angle.h b/src/angle.h index 139380ff39..ffed437743 100644 --- a/src/angle.h +++ b/src/angle.h @@ -30,6 +30,11 @@ class Angle : protected Pointers { double *eatom,**vatom; // accumulated per-atom energy/virial double **cvatom; // accumulated per-atom centroid virial + int centroidstressflag; // centroid stress compared to two-body stress + // CENTROID_SAME = same as two-body stress + // CENTROID_AVAIL = different and implemented + // CENTROID_NOTAVAIL = different, not yet implemented + // KOKKOS host/device flag and data masks ExecutionSpace execution_space; diff --git a/src/bond.cpp b/src/bond.cpp index 1912d15aa7..16d3c43dea 100644 --- a/src/bond.cpp +++ b/src/bond.cpp @@ -80,7 +80,16 @@ void Bond::init() /* ---------------------------------------------------------------------- setup for energy, virial computation - see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) + see integrate::ev_set() for bitwise settings of eflag/vflag + set the following flags, values are otherwise set to 0: + evflag != 0 if any bits of eflag or vflag are set + eflag_global != 0 if ENERGY_GLOBAL bit of eflag set + eflag_atom != 0 if ENERGY_ATOM bit of eflag set + eflag_either != 0 if eflag_global or eflag_atom is set + vflag_global != 0 if VIRIAL_PAIR or VIRIAL_FDOTR bit of vflag set + vflag_atom != 0 if VIRIAL_ATOM or VIRIAL_CENTROID bit of vflag set + two-body and centroid stress are identical for bonds + vflag_either != 0 if vflag_global or vflag_atom is set ------------------------------------------------------------------------- */ void Bond::ev_setup(int eflag, int vflag, int alloc) @@ -90,13 +99,12 @@ void Bond::ev_setup(int eflag, int vflag, int alloc) evflag = 1; eflag_either = eflag; - eflag_global = eflag % 2; - eflag_atom = eflag / 2; + eflag_global = eflag & ENERGY_GLOBAL; + eflag_atom = eflag & ENERGY_ATOM; vflag_either = vflag; - vflag_global = vflag % 4; - // per-atom virial and per-atom centroid virial are the same for bonds - vflag_atom = vflag / 4; + vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); + vflag_atom = vflag & (VIRIAL_ATOM | VIRIAL_CENTROID); // reallocate per-atom arrays if necessary diff --git a/src/compute_centroid_stress_atom.cpp b/src/compute_centroid_stress_atom.cpp index d17d402a3d..a6755b5b53 100644 --- a/src/compute_centroid_stress_atom.cpp +++ b/src/compute_centroid_stress_atom.cpp @@ -124,10 +124,35 @@ void ComputeCentroidStressAtom::init() else biasflag = NOBIAS; } else biasflag = NOBIAS; - // check if pair styles support centroid atom stress + // check if force components and fixes support centroid atom stress + // all bond styles support it as CENTROID_SAME + if (pairflag && force->pair) - if (force->pair->centroidstressflag & 4) + if (force->pair->centroidstressflag == CENTROID_NOTAVAIL) error->all(FLERR, "Pair style does not support compute centroid/stress/atom"); + + if (angleflag && force->angle) + if (force->angle->centroidstressflag == CENTROID_NOTAVAIL) + error->all(FLERR, "Angle style does not support compute centroid/stress/atom"); + + if (dihedralflag && force->dihedral) + if (force->dihedral->centroidstressflag == CENTROID_NOTAVAIL) + error->all(FLERR, "Dihedral style does not support compute centroid/stress/atom"); + + if (improperflag && force->improper) + if (force->improper->centroidstressflag == CENTROID_NOTAVAIL) + error->all(FLERR, "Improper style does not support compute centroid/stress/atom"); + + if (kspaceflag && force->kspace) + if (force->kspace->centroidstressflag == CENTROID_NOTAVAIL) + error->all(FLERR, "KSpace style does not support compute centroid/stress/atom"); + + if (fixflag) { + for (int ifix = 0; ifix < modify->nfix; ifix++) + if (modify->fix[ifix]->virial_flag && + modify->fix[ifix]->centroidstressflag == CENTROID_NOTAVAIL) + error->all(FLERR, "Fix style does not support compute centroid/stress/atom"); + } } /* ---------------------------------------------------------------------- */ @@ -173,12 +198,12 @@ void ComputeCentroidStressAtom::compute_peratom() for (j = 0; j < 9; j++) stress[i][j] = 0.0; - // add in per-atom contributions from each force - - // per-atom virial and per-atom centroid virial are the same for two-body - // many-body pair styles not yet implemented + // add in per-atom contributions from all force components and fixes + + // pair styles are either CENTROID_SAME or CENTROID_AVAIL or CENTROID_NOTAVAIL + if (pairflag && force->pair && force->pair->compute_flag) { - if (force->pair->centroidstressflag & 2) { + if (force->pair->centroidstressflag == CENTROID_AVAIL) { double **cvatom = force->pair->cvatom; for (i = 0; i < npair; i++) for (j = 0; j < 9; j++) @@ -195,6 +220,10 @@ void ComputeCentroidStressAtom::compute_peratom() } // per-atom virial and per-atom centroid virial are the same for bonds + // bond styles are all CENTROID_SAME + // angle, dihedral, improper styles are CENTROID_AVAIL or CENTROID_NOTAVAIL + // KSpace styles are all CENTROID_NOTAVAIL, placeholder CENTROID_SAME below + if (bondflag && force->bond) { double **vatom = force->bond->vatom; for (i = 0; i < nbond; i++) { @@ -241,7 +270,8 @@ void ComputeCentroidStressAtom::compute_peratom() // possible during setup phase if fix has not initialized its vatom yet // e.g. fix ave/spatial defined before fix shake, // and fix ave/spatial uses a per-atom stress from this compute as input - + // fix styles are CENTROID_SAME or CENTROID_NOTAVAIL + if (fixflag) { for (int ifix = 0; ifix < modify->nfix; ifix++) if (modify->fix[ifix]->virial_flag) { @@ -258,7 +288,8 @@ void ComputeCentroidStressAtom::compute_peratom() // communicate ghost virials between neighbor procs - if (force->newton || (force->kspace && force->kspace->tip4pflag && force->kspace->compute_flag)) + if (force->newton || + (force->kspace && force->kspace->tip4pflag && force->kspace->compute_flag)) comm->reverse_comm_compute(this); // zero virial of atoms not in group diff --git a/src/dihedral.cpp b/src/dihedral.cpp index f778341bfa..4680b6f67d 100644 --- a/src/dihedral.cpp +++ b/src/dihedral.cpp @@ -40,6 +40,7 @@ Dihedral::Dihedral(LAMMPS *lmp) : Pointers(lmp) vatom = nullptr; cvatom = nullptr; setflag = nullptr; + centroidstressflag = CENTROID_AVAIL; execution_space = Host; datamask_read = ALL_MASK; @@ -74,7 +75,19 @@ void Dihedral::init() /* ---------------------------------------------------------------------- setup for energy, virial computation - see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) + see integrate::ev_set() for bitwise settings of eflag/vflag + set the following flags, values are otherwise set to 0: + evflag != 0 if any bits of eflag or vflag are set + eflag_global != 0 if ENERGY_GLOBAL bit of eflag set + eflag_atom != 0 if ENERGY_ATOM bit of eflag set + eflag_either != 0 if eflag_global or eflag_atom is set + vflag_global != 0 if VIRIAL_PAIR or VIRIAL_FDOTR bit of vflag set + vflag_atom != 0 if VIRIAL_ATOM bit of vflag set + vflag_atom != 0 if VIRIAL_CENTROID bit of vflag set + and centroidstressflag != CENTROID_AVAIL + cvflag_atom != 0 if VIRIAL_CENTROID bit of vflag set + and centroidstressflag = CENTROID_AVAIL + vflag_either != 0 if any of vflag_global, vflag_atom, cvflag_atom is set ------------------------------------------------------------------------- */ void Dihedral::ev_setup(int eflag, int vflag, int alloc) @@ -84,13 +97,17 @@ void Dihedral::ev_setup(int eflag, int vflag, int alloc) evflag = 1; eflag_either = eflag; - eflag_global = eflag % 2; - eflag_atom = eflag / 2; + eflag_global = eflag & ENERGY_GLOBAL; + eflag_atom = eflag & ENERGY_ATOM; - vflag_global = vflag % 4; - vflag_atom = vflag & 4; - cvflag_atom = vflag & 8; - vflag_either = vflag_global || vflag_atom; + vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); + vflag_atom = vflag & VIRIAL_ATOM; + if (vflag & VIRIAL_CENTROID && centroidstressflag != CENTROID_AVAIL) + vflag_atom = 1; + cvflag_atom = 0; + if (vflag & VIRIAL_CENTROID && centroidstressflag == CENTROID_AVAIL) + cvflag_atom = 1; + vflag_either = vflag_global || vflag_atom || cvflag_atom; // reallocate per-atom arrays if necessary diff --git a/src/dihedral.h b/src/dihedral.h index c96891b1c4..c571a74dd4 100644 --- a/src/dihedral.h +++ b/src/dihedral.h @@ -30,6 +30,11 @@ class Dihedral : protected Pointers { double *eatom,**vatom; // accumulated per-atom energy/virial double **cvatom; // accumulated per-atom centroid virial + int centroidstressflag; // centroid stress compared to two-body stress + // CENTROID_SAME = same as two-body stress + // CENTROID_AVAIL = different and implemented + // CENTROID_NOTAVAIL = different, not yet implemented + // KOKKOS host/device flag and data masks ExecutionSpace execution_space; diff --git a/src/fix.cpp b/src/fix.cpp index 905374ee05..ae8f8c7369 100644 --- a/src/fix.cpp +++ b/src/fix.cpp @@ -16,6 +16,7 @@ #include "atom.h" #include "atom_masks.h" #include "error.h" +#include "force.h" #include "group.h" #include "memory.h" @@ -103,7 +104,8 @@ Fix::Fix(LAMMPS *lmp, int /*narg*/, char **arg) : maxeatom = maxvatom = 0; vflag_atom = 0; - + centroidstressflag = CENTROID_SAME; + // KOKKOS per-fix data masks execution_space = Host; @@ -190,12 +192,12 @@ void Fix::ev_setup(int eflag, int vflag) evflag = 1; eflag_either = eflag; - eflag_global = eflag % 2; - eflag_atom = eflag / 2; + eflag_global = eflag & ENERGY_GLOBAL; + eflag_atom = eflag & ENERGY_ATOM; vflag_either = vflag; - vflag_global = vflag % 4; - vflag_atom = vflag / 4; + vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); + vflag_atom = vflag & (VIRIAL_ATOM | VIRIAL_CENTROID); // reallocate per-atom arrays if necessary @@ -235,7 +237,7 @@ void Fix::ev_setup(int eflag, int vflag) /* ---------------------------------------------------------------------- if thermo_virial is on: setup for virial computation - see integrate::ev_set() for values of vflag (0-6) + see integrate::ev_set() for values of vflag fixes call this if use v_tally() else: set evflag=0 ------------------------------------------------------------------------- */ @@ -251,8 +253,8 @@ void Fix::v_setup(int vflag) evflag = 1; - vflag_global = vflag % 4; - vflag_atom = vflag / 4; + vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); + vflag_atom = vflag & (VIRIAL_ATOM | VIRIAL_CENTROID); // reallocate per-atom array if necessary diff --git a/src/fix.h b/src/fix.h index 1a28e5c924..75c644e87f 100644 --- a/src/fix.h +++ b/src/fix.h @@ -102,6 +102,11 @@ class Fix : protected Pointers { double virial[6]; // accumulated virial double *eatom,**vatom; // accumulated per-atom energy/virial + int centroidstressflag; // centroid stress compared to two-body stress + // CENTROID_SAME = same as two-body stress + // CENTROID_AVAIL = different and implemented + // CENTROID_NOTAVAIL = different, not yet implemented + int restart_reset; // 1 if restart just re-initialized fix // KOKKOS host/device flag and data masks diff --git a/src/force.h b/src/force.h index a2cc37f0b8..b087b8e187 100644 --- a/src/force.h +++ b/src/force.h @@ -26,6 +26,26 @@ namespace LAMMPS_NS { class KSpace; class Pair; +enum { + ENERGY_NONE = 0x00, + ENERGY_GLOBAL = 0x01, + ENERGY_ATOM = 0x02 +}; + +enum { + VIRIAL_NONE = 0x00, + VIRIAL_PAIR = 0x01, + VIRIAL_FDOTR = 0x02, + VIRIAL_ATOM = 0x04, + VIRIAL_CENTROID = 0x08 +}; + +enum { + CENTROID_SAME = 0, + CENTROID_AVAIL = 1, + CENTROID_NOTAVAIL = 2 +}; + class Force : protected Pointers { public: double boltz; // Boltzmann constant (eng/degree-K) diff --git a/src/improper.cpp b/src/improper.cpp index faf9e2fa52..6e0f829733 100644 --- a/src/improper.cpp +++ b/src/improper.cpp @@ -37,6 +37,7 @@ Improper::Improper(LAMMPS *lmp) : Pointers(lmp) vatom = nullptr; cvatom = nullptr; setflag = nullptr; + centroidstressflag = CENTROID_AVAIL; execution_space = Host; datamask_read = ALL_MASK; @@ -72,7 +73,19 @@ void Improper::init() /* ---------------------------------------------------------------------- setup for energy, virial computation - see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) + see integrate::ev_set() for bitwise settings of eflag/vflag + set the following flags, values are otherwise set to 0: + evflag != 0 if any bits of eflag or vflag are set + eflag_global != 0 if ENERGY_GLOBAL bit of eflag set + eflag_atom != 0 if ENERGY_ATOM bit of eflag set + eflag_either != 0 if eflag_global or eflag_atom is set + vflag_global != 0 if VIRIAL_PAIR or VIRIAL_FDOTR bit of vflag set + vflag_atom != 0 if VIRIAL_ATOM bit of vflag set + vflag_atom != 0 if VIRIAL_CENTROID bit of vflag set + and centroidstressflag != CENTROID_AVAIL + cvflag_atom != 0 if VIRIAL_CENTROID bit of vflag set + and centroidstressflag = CENTROID_AVAIL + vflag_either != 0 if any of vflag_global, vflag_atom, cvflag_atom is set ------------------------------------------------------------------------- */ void Improper::ev_setup(int eflag, int vflag, int alloc) @@ -82,13 +95,17 @@ void Improper::ev_setup(int eflag, int vflag, int alloc) evflag = 1; eflag_either = eflag; - eflag_global = eflag % 2; - eflag_atom = eflag / 2; + eflag_global = eflag & ENERGY_GLOBAL; + eflag_atom = eflag & ENERGY_ATOM; - vflag_global = vflag % 4; - vflag_atom = vflag & 4; - cvflag_atom = vflag & 8; - vflag_either = vflag_global || vflag_atom; + vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); + vflag_atom = vflag & VIRIAL_ATOM; + if (vflag & VIRIAL_CENTROID && centroidstressflag != CENTROID_AVAIL) + vflag_atom = 1; + cvflag_atom = 0; + if (vflag & VIRIAL_CENTROID && centroidstressflag == CENTROID_AVAIL) + cvflag_atom = 1; + vflag_either = vflag_global || vflag_atom || cvflag_atom; // reallocate per-atom arrays if necessary diff --git a/src/improper.h b/src/improper.h index f3bed2607c..9e73be931c 100644 --- a/src/improper.h +++ b/src/improper.h @@ -30,6 +30,11 @@ class Improper : protected Pointers { double *eatom,**vatom; // accumulated per-atom energy/virial double **cvatom; // accumulated per-atom centroid virial + int centroidstressflag; // centroid stress compared to two-body stress + // CENTROID_SAME = same as two-body stress + // CENTROID_AVAIL = different and implemented + // CENTROID_NOTAVAIL = different, not yet implemented + // KOKKOS host/device flag and data masks ExecutionSpace execution_space; diff --git a/src/integrate.cpp b/src/integrate.cpp index 7d4bf36929..5c700527a1 100644 --- a/src/integrate.cpp +++ b/src/integrate.cpp @@ -111,21 +111,18 @@ void Integrate::ev_setup() /* ---------------------------------------------------------------------- set eflag,vflag for current iteration + based on computes that need energy/virial info on this timestep invoke matchstep() on all timestep-dependent computes to clear their arrays - eflag/vflag based on computes that need info on this ntimestep - eflag = 0 = no energy computation - eflag = 1 = global energy only - eflag = 2 = per-atom energy only - eflag = 3 = both global and per-atom energy - vflag = 0 = no virial computation (pressure) - vflag = 1 = global virial with pair portion via sum of pairwise interactions - vflag = 2 = global virial with pair portion via F dot r including ghosts - vflag = 4 = per-atom virial only - vflag = 5 or 6 = both global and per-atom virial - vflag = 8 = per-atom centroid virial only - vflag = 9 or 10 = both global and per-atom centroid virial - vflag = 12 = both per-atom virial and per-atom centroid virial - vflag = 13 or 15 = global, per-atom virial and per-atom centroid virial + eflag: set any or no bits + ENERGY_GLOBAL bit for global energy + ENERGY_ATOM bit for per-atom energy + vflag: set any or no bits, but PAIR/FDOTR bits cannot both be set + VIRIAL_PAIR bit for global virial as sum of pairwise terms + VIRIAL_FDOTR bit for global virial via F dot r + VIRIAL_ATOM bit for per-atom virial + VIRIAL_CENTROID bit for per-atom centroid virial + all force components (pair,bond,angle,...,kspace) use eflag/vflag + in their ev_setup() method to set local energy/virial flags ------------------------------------------------------------------------- */ void Integrate::ev_set(bigint ntimestep) @@ -136,13 +133,13 @@ void Integrate::ev_set(bigint ntimestep) int eflag_global = 0; for (i = 0; i < nelist_global; i++) if (elist_global[i]->matchstep(ntimestep)) flag = 1; - if (flag) eflag_global = 1; + if (flag) eflag_global = ENERGY_GLOBAL; flag = 0; int eflag_atom = 0; for (i = 0; i < nelist_atom; i++) if (elist_atom[i]->matchstep(ntimestep)) flag = 1; - if (flag) eflag_atom = 2; + if (flag) eflag_atom = ENERGY_ATOM; if (eflag_global) update->eflag_global = ntimestep; if (eflag_atom) update->eflag_atom = ntimestep; @@ -158,13 +155,13 @@ void Integrate::ev_set(bigint ntimestep) int vflag_atom = 0; for (i = 0; i < nvlist_atom; i++) if (vlist_atom[i]->matchstep(ntimestep)) flag = 1; - if (flag) vflag_atom = 4; + if (flag) vflag_atom = VIRIAL_ATOM; flag = 0; int cvflag_atom = 0; for (i = 0; i < ncvlist_atom; i++) if (cvlist_atom[i]->matchstep(ntimestep)) flag = 1; - if (flag) cvflag_atom = 8; + if (flag) cvflag_atom = VIRIAL_CENTROID; if (vflag_global) update->vflag_global = ntimestep; if (vflag_atom || cvflag_atom) update->vflag_atom = ntimestep; diff --git a/src/kspace.cpp b/src/kspace.cpp index 16668c137a..594fd1ba8d 100644 --- a/src/kspace.cpp +++ b/src/kspace.cpp @@ -88,7 +88,8 @@ KSpace::KSpace(LAMMPS *lmp) : Pointers(lmp) maxeatom = maxvatom = 0; eatom = nullptr; vatom = nullptr; - + centroidstressflag = CENTROID_NOTAVAIL; + execution_space = Host; datamask_read = ALL_MASK; datamask_modify = ALL_MASK; @@ -214,7 +215,17 @@ void KSpace::pair_check() /* ---------------------------------------------------------------------- setup for energy, virial computation - see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) + see integrate::ev_set() for bitwise settings of eflag/vflag + set the following flags, values are otherwise set to 0: + evflag != 0 if any bits of eflag or vflag are set + eflag_global != 0 if ENERGY_GLOBAL bit of eflag set + eflag_atom != 0 if ENERGY_ATOM bit of eflag set + eflag_either != 0 if eflag_global or eflag_atom is set + vflag_global != 0 if VIRIAL_PAIR or VIRIAL_FDOTR bit of vflag set + vflag_atom != 0 if VIRIAL_ATOM bit of vflag set + no current support for centroid stress + vflag_either != 0 if vflag_global or vflag_atom is set + evflag_atom != 0 if eflag_atom or vflag_atom is set ------------------------------------------------------------------------- */ void KSpace::ev_setup(int eflag, int vflag, int alloc) @@ -224,12 +235,12 @@ void KSpace::ev_setup(int eflag, int vflag, int alloc) evflag = 1; eflag_either = eflag; - eflag_global = eflag % 2; - eflag_atom = eflag / 2; + eflag_global = eflag & ENERGY_GLOBAL; + eflag_atom = eflag & ENERGY_ATOM; vflag_either = vflag; - vflag_global = vflag % 4; - vflag_atom = vflag / 4; + vflag_global = vflag & (VIRIAL_PAIR | VIRIAL_FDOTR); + vflag_atom = vflag & VIRIAL_ATOM; if (eflag_atom || vflag_atom) evflag_atom = 1; else evflag_atom = 0; diff --git a/src/kspace.h b/src/kspace.h index 110a790dfe..978daeace1 100644 --- a/src/kspace.h +++ b/src/kspace.h @@ -80,6 +80,11 @@ class KSpace : protected Pointers { int group_group_enable; // 1 if style supports group/group calculation + int centroidstressflag; // centroid stress compared to two-body stress + // CENTROID_SAME = same as two-body stress + // CENTROID_AVAIL = different and implemented + // CENTROID_NOTAVAIL = different, not yet implemented + // KOKKOS host/device flag and data masks ExecutionSpace execution_space; diff --git a/src/min.cpp b/src/min.cpp index ac711d8900..ca0d508e95 100644 --- a/src/min.cpp +++ b/src/min.cpp @@ -144,11 +144,12 @@ void Min::init() requestor = nullptr; // virial_style: - // 1 if computed explicitly by pair->compute via sum over pair interactions - // 2 if computed implicitly by pair->virial_compute via sum over ghost atoms + // VIRIAL_PAIR if computed explicitly in pair via sum over pair interactions + // VIRIAL_FDOTR if computed implicitly in pair by + // virial_fdotr_compute() via sum over ghosts - if (force->newton_pair) virial_style = 2; - else virial_style = 1; + if (force->newton_pair) virial_style = VIRIAL_FDOTR; + else virial_style = VIRIAL_PAIR; // setup lists of computes for global and per-atom PE and pressure @@ -801,19 +802,16 @@ void Min::ev_setup() invoke matchstep() on all timestep-dependent computes to clear their arrays eflag/vflag based on computes that need info on this ntimestep always set eflag_global = 1, since need energy every iteration - eflag = 0 = no energy computation - eflag = 1 = global energy only - eflag = 2 = per-atom energy only - eflag = 3 = both global and per-atom energy - vflag = 0 = no virial computation (pressure) - vflag = 1 = global virial with pair portion via sum of pairwise interactions - vflag = 2 = global virial with pair portion via F dot r including ghosts - vflag = 4 = per-atom virial only - vflag = 5 or 6 = both global and per-atom virial - vflag = 8 = per-atom centroid virial only - vflag = 9 or 10 = both global and per-atom centroid virial - vflag = 12 = both per-atom virial and per-atom centroid virial - vflag = 13 or 15 = global, per-atom virial and per-atom centroid virial + eflag: set any or no bits + ENERGY_GLOBAL bit for global energy + ENERGY_ATOM bit for per-atom energy + vflag: set any or no bits, but GLOBAL/FDOTR bit cannot both be set + VIRIAL_PAIR bit for global virial as sum of pairwise terms + VIRIAL_FDOTR bit for global virial via F dot r + VIRIAL_ATOM bit for per-atom virial + VIRIAL_CENTROID bit for per-atom centroid virial + all force components (pair,bond,angle,...,kspace) use eflag/vflag + in their ev_setup() method to set local energy/virial flags ------------------------------------------------------------------------- */ void Min::ev_set(bigint ntimestep) @@ -828,7 +826,7 @@ void Min::ev_set(bigint ntimestep) int eflag_atom = 0; for (i = 0; i < nelist_atom; i++) if (elist_atom[i]->matchstep(ntimestep)) flag = 1; - if (flag) eflag_atom = 2; + if (flag) eflag_atom = ENERGY_ATOM; if (eflag_global) update->eflag_global = update->ntimestep; if (eflag_atom) update->eflag_atom = update->ntimestep; @@ -844,13 +842,13 @@ void Min::ev_set(bigint ntimestep) int vflag_atom = 0; for (i = 0; i < nvlist_atom; i++) if (vlist_atom[i]->matchstep(ntimestep)) flag = 1; - if (flag) vflag_atom = 4; + if (flag) vflag_atom = VIRIAL_ATOM; flag = 0; int cvflag_atom = 0; for (i = 0; i < ncvlist_atom; i++) if (cvlist_atom[i]->matchstep(ntimestep)) flag = 1; - if (flag) cvflag_atom = 8; + if (flag) cvflag_atom = VIRIAL_CENTROID; if (vflag_global) update->vflag_global = update->ntimestep; if (vflag_atom || cvflag_atom) update->vflag_atom = update->ntimestep; diff --git a/src/pair.cpp b/src/pair.cpp index 1ff2322b48..74a03fbed9 100644 --- a/src/pair.cpp +++ b/src/pair.cpp @@ -74,9 +74,10 @@ Pair::Pair(LAMMPS *lmp) : Pointers(lmp) setflag = nullptr; cutsq = nullptr; - ewaldflag = pppmflag = msmflag = dispersionflag = tip4pflag = dipoleflag = spinflag = 0; + ewaldflag = pppmflag = msmflag = dispersionflag = + tip4pflag = dipoleflag = spinflag = 0; reinitflag = 1; - centroidstressflag = 4; + centroidstressflag = CENTROID_SAME; // pair_modify settings @@ -776,35 +777,47 @@ void Pair::del_tally_callback(Compute *ptr) /* ---------------------------------------------------------------------- setup for energy, virial computation - see integrate::ev_set() for values of eflag (0-3) and vflag (0-6) + see integrate::ev_set() for bitwise settings of eflag/vflag + set the following flags, values are otherwise set to 0: + eflag_global != 0 if ENERGY_GLOBAL bit of eflag set + eflag_atom != 0 if ENERGY_ATOM bit of eflag set + eflag_either != 0 if eflag_global or eflag_atom is set + vflag_global != 0 if VIRIAL_PAIR bit of vflag set + vflag_global != 0 if VIRIAL_FDOTR bit of vflag is set but no_virial_fdotr = 1 + vflag_fdotr != 0 if VIRIAL_FDOTR bit of vflag set and no_virial_fdotr = 0 + vflag_atom != 0 if VIRIAL_ATOM bit of vflag set + vflag_atom != 0 if VIRIAL_CENTROID bit of vflag set + and centroidstressflag != CENTROID_AVAIL + cvflag_atom != 0 if VIRIAL_CENTROID bit of vflag set + and centroidstressflag = CENTROID_AVAIL + vflag_either != 0 if any of vflag_global, vflag_atom, cvflag_atom is set + evflag != 0 if eflag_either or vflag_either is set + centroidstressflag is set by the pair style to one of these values: + CENTROID_SAME = same as two-body stress + CENTROID_AVAIL = different and implemented + CENTROID_NOTAVAIL = different but not yet implemented ------------------------------------------------------------------------- */ void Pair::ev_setup(int eflag, int vflag, int alloc) { int i,n; - evflag = 1; - eflag_either = eflag; - eflag_global = eflag % 2; - eflag_atom = eflag / 2; - - vflag_global = vflag % 4; - vflag_atom = vflag & 4; + eflag_global = eflag & ENERGY_GLOBAL; + eflag_atom = eflag & ENERGY_ATOM; + + vflag_global = vflag & VIRIAL_PAIR; + if (vflag & VIRIAL_FDOTR && no_virial_fdotr_compute == 1) vflag_global = 1; + if (vflag & VIRIAL_FDOTR && no_virial_fdotr_compute == 0) vflag_fdotr = 1; + vflag_atom = vflag & VIRIAL_ATOM; + if (vflag & VIRIAL_CENTROID && centroidstressflag != CENTROID_AVAIL) vflag_atom = 1; cvflag_atom = 0; + if (vflag & VIRIAL_CENTROID && centroidstressflag == CENTROID_AVAIL) cvflag_atom = 1; + vflag_either = vflag_global || vflag_atom || cvflag_atom; - if (vflag & 8) { - if (centroidstressflag & 2) { - cvflag_atom = 1; - } else { - vflag_atom = 1; - } - // extra check, because both bits might be set - if (centroidstressflag & 1) vflag_atom = 1; - } - - vflag_either = vflag_global || vflag_atom; - + evflag = eflag_either || vflag_either; + if (!evflag) return; + // reallocate per-atom arrays if necessary if (eflag_atom && atom->nmax > maxeatom) { @@ -869,18 +882,6 @@ void Pair::ev_setup(int eflag, int vflag, int alloc) } } - // if vflag_global = 2 and pair::compute() calls virial_fdotr_compute() - // compute global virial via (F dot r) instead of via pairwise summation - // unset other flags as appropriate - - if (vflag_global == 2 && no_virial_fdotr_compute == 0) { - vflag_fdotr = 1; - vflag_global = 0; - if (vflag_atom == 0 && cvflag_atom == 0) vflag_either = 0; - if (vflag_either == 0 && eflag_either == 0) evflag = 0; - } else vflag_fdotr = 0; - - // run ev_setup option for USER-TALLY computes if (num_tally_compute > 0) { @@ -1563,6 +1564,8 @@ void Pair::virial_fdotr_compute() double **x = atom->x; double **f = atom->f; + for (int i = 0; i < 6; i++) virial[i] = 0.0; + // sum over force on all particles including ghosts if (neighbor->includegroup == 0) { diff --git a/src/pair.h b/src/pair.h index ace233539f..ff0ce82cc8 100644 --- a/src/pair.h +++ b/src/pair.h @@ -68,10 +68,10 @@ class Pair : protected Pointers { int spinflag; // 1 if compatible with spin solver int reinitflag; // 1 if compatible with fix adapt and alike - int centroidstressflag; // compatibility with centroid atomic stress - // 1 if same as two-body atomic stress - // 2 if implemented and different from two-body - // 4 if not compatible/implemented + int centroidstressflag; // centroid stress compared to two-body stress + // CENTROID_SAME = same as two-body stress + // CENTROID_AVAIL = different and implemented + // CENTROID_NOTAVAIL = different, not yet implemented int tail_flag; // pair_modify flag for LJ tail correction double etail,ptail; // energy/pressure tail corrections diff --git a/src/pair_beck.cpp b/src/pair_beck.cpp index bfd39a56e5..a6e63fc4ec 100644 --- a/src/pair_beck.cpp +++ b/src/pair_beck.cpp @@ -26,15 +26,12 @@ #include "error.h" #include "math_special.h" - using namespace LAMMPS_NS; using namespace MathSpecial; /* ---------------------------------------------------------------------- */ -PairBeck::PairBeck(LAMMPS *lmp) : Pair(lmp) { - centroidstressflag = 1; -} +PairBeck::PairBeck(LAMMPS *lmp) : Pair(lmp) {} /* ---------------------------------------------------------------------- */ diff --git a/src/pair_buck.cpp b/src/pair_buck.cpp index 1954ae2891..989c341add 100644 --- a/src/pair_buck.cpp +++ b/src/pair_buck.cpp @@ -32,7 +32,6 @@ using namespace MathConst; PairBuck::PairBuck(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_buck_coul_cut.cpp b/src/pair_buck_coul_cut.cpp index 0940a3232a..05b36d3620 100644 --- a/src/pair_buck_coul_cut.cpp +++ b/src/pair_buck_coul_cut.cpp @@ -36,7 +36,6 @@ using namespace MathConst; PairBuckCoulCut::PairBuckCoulCut(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_coul_cut.cpp b/src/pair_coul_cut.cpp index c50cc75c0c..44ddd424d9 100644 --- a/src/pair_coul_cut.cpp +++ b/src/pair_coul_cut.cpp @@ -28,9 +28,7 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -PairCoulCut::PairCoulCut(LAMMPS *lmp) : Pair(lmp) { - centroidstressflag = 1; -} +PairCoulCut::PairCoulCut(LAMMPS *lmp) : Pair(lmp) {} /* ---------------------------------------------------------------------- */ diff --git a/src/pair_coul_dsf.cpp b/src/pair_coul_dsf.cpp index 889f2757aa..20bae21594 100644 --- a/src/pair_coul_dsf.cpp +++ b/src/pair_coul_dsf.cpp @@ -29,7 +29,6 @@ #include "math_const.h" #include "error.h" - using namespace LAMMPS_NS; using namespace MathConst; @@ -43,9 +42,7 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -PairCoulDSF::PairCoulDSF(LAMMPS *lmp) : Pair(lmp) { - centroidstressflag = 1; -} +PairCoulDSF::PairCoulDSF(LAMMPS *lmp) : Pair(lmp) {} /* ---------------------------------------------------------------------- */ diff --git a/src/pair_coul_wolf.cpp b/src/pair_coul_wolf.cpp index df22b5a67d..4466651237 100644 --- a/src/pair_coul_wolf.cpp +++ b/src/pair_coul_wolf.cpp @@ -36,7 +36,6 @@ using namespace MathConst; PairCoulWolf::PairCoulWolf(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; // NOTE: single() method below is not yet correct - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_hybrid.cpp b/src/pair_hybrid.cpp index 8e69deb21a..07206dc518 100644 --- a/src/pair_hybrid.cpp +++ b/src/pair_hybrid.cpp @@ -29,7 +29,6 @@ #include "suffix.h" - using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ @@ -42,10 +41,6 @@ PairHybrid::PairHybrid(LAMMPS *lmp) : Pair(lmp), outerflag = 0; respaflag = 0; - - // assume pair hybrid always supports centroid atomic stress, - // so that cflag_atom gets set when needed - centroidstressflag = 2; } /* ---------------------------------------------------------------------- */ @@ -82,10 +77,10 @@ PairHybrid::~PairHybrid() /* ---------------------------------------------------------------------- call each sub-style's compute() or compute_outer() function accumulate sub-style global/peratom energy/virial in hybrid - for global vflag = 1: + for global vflag = VIRIAL_PAIR: each sub-style computes own virial[6] sum sub-style virial[6] to hybrid's virial[6] - for global vflag = 2: + for global vflag = VIRIAL_FDOTR: call sub-style with adjusted vflag to prevent it calling virial_fdotr_compute() hybrid calls virial_fdotr_compute() on final accumulated f @@ -95,21 +90,22 @@ void PairHybrid::compute(int eflag, int vflag) { int i,j,m,n; - // if no_virial_fdotr_compute is set and global component of - // incoming vflag = 2, then - // reset vflag as if global component were 1 + // check if no_virial_fdotr_compute is set and global component of + // incoming vflag = VIRIAL_FDOTR + // if so, reset vflag as if global component were VIRIAL_PAIR // necessary since one or more sub-styles cannot compute virial as F dot r - if (no_virial_fdotr_compute && vflag % 4 == 2) vflag = 1 + vflag/4 * 4; + if (no_virial_fdotr_compute && (vflag & VIRIAL_FDOTR)) + vflag = VIRIAL_PAIR | (vflag & ~VIRIAL_FDOTR); ev_init(eflag,vflag); - // check if global component of incoming vflag = 2 - // if so, reset vflag passed to substyle as if it were 0 + // check if global component of incoming vflag = VIRIAL_FDOTR + // if so, reset vflag passed to substyle so VIRIAL_FDOTR is turned off // necessary so substyle will not invoke virial_fdotr_compute() int vflag_substyle; - if (vflag % 4 == 2) vflag_substyle = vflag/4 * 4; + if (vflag & VIRIAL_FDOTR) vflag_substyle = vflag & ~VIRIAL_FDOTR; else vflag_substyle = vflag; double *saved_special = save_special(); @@ -165,10 +161,13 @@ void PairHybrid::compute(int eflag, int vflag) for (j = 0; j < 6; j++) vatom[i][j] += vatom_substyle[i][j]; } + + // substyles may be CENTROID_SAME or CENTROID_AVAIL + if (cvflag_atom) { n = atom->nlocal; if (force->newton_pair) n += atom->nghost; - if (styles[m]->centroidstressflag & 2) { + if (styles[m]->centroidstressflag == CENTROID_AVAIL) { double **cvatom_substyle = styles[m]->cvatom; for (i = 0; i < n; i++) for (j = 0; j < 9; j++) @@ -386,6 +385,7 @@ void PairHybrid::flags() compute_flag = 0; respa_enable = 0; restartinfo = 0; + for (m = 0; m < nstyles; m++) { if (styles[m]->single_enable) ++single_enable; if (styles[m]->respa_enable) ++respa_enable; @@ -401,12 +401,26 @@ void PairHybrid::flags() if (styles[m]->dispersionflag) dispersionflag = 1; if (styles[m]->tip4pflag) tip4pflag = 1; if (styles[m]->compute_flag) compute_flag = 1; - if (styles[m]->centroidstressflag & 4) centroidstressflag |= 4; } single_enable = (single_enable == nstyles) ? 1 : 0; respa_enable = (respa_enable == nstyles) ? 1 : 0; restartinfo = (restartinfo == nstyles) ? 1 : 0; init_svector(); + + // set centroidstressflag for pair hybrid + // set to CENTROID_NOTAVAIL if any substyle is NOTAVAIL + // else set to CENTROID_AVAIL if any substyle is AVAIL + // else set to CENTROID_SAME if all substyles are SAME + + centroidstressflag = CENTROID_SAME; + + for (m = 0; m < nstyles; m++) { + if (styles[m]->centroidstressflag == CENTROID_NOTAVAIL) + centroidstressflag = CENTROID_NOTAVAIL; + if (centroidstressflag == CENTROID_SAME && + styles[m]->centroidstressflag == CENTROID_AVAIL) + centroidstressflag = CENTROID_AVAIL; + } } /* ---------------------------------------------------------------------- diff --git a/src/pair_lj_cubic.cpp b/src/pair_lj_cubic.cpp index 73e99d5ff8..7fe4757358 100644 --- a/src/pair_lj_cubic.cpp +++ b/src/pair_lj_cubic.cpp @@ -32,9 +32,7 @@ using namespace PairLJCubicConstants; /* ---------------------------------------------------------------------- */ -PairLJCubic::PairLJCubic(LAMMPS *lmp) : Pair(lmp) { - centroidstressflag = 1; -} +PairLJCubic::PairLJCubic(LAMMPS *lmp) : Pair(lmp) {} /* ---------------------------------------------------------------------- */ diff --git a/src/pair_lj_cut.cpp b/src/pair_lj_cut.cpp index cdd9167b65..32ac15da79 100644 --- a/src/pair_lj_cut.cpp +++ b/src/pair_lj_cut.cpp @@ -31,7 +31,6 @@ #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; using namespace MathConst; @@ -41,7 +40,6 @@ PairLJCut::PairLJCut(LAMMPS *lmp) : Pair(lmp) { respa_enable = 1; writedata = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_lj_cut_coul_cut.cpp b/src/pair_lj_cut_coul_cut.cpp index 2b164ede1d..9ed84bd1a4 100644 --- a/src/pair_lj_cut_coul_cut.cpp +++ b/src/pair_lj_cut_coul_cut.cpp @@ -33,7 +33,6 @@ using namespace MathConst; PairLJCutCoulCut::PairLJCutCoulCut(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_lj_cut_coul_dsf.cpp b/src/pair_lj_cut_coul_dsf.cpp index d1989ba020..3382ca1986 100644 --- a/src/pair_lj_cut_coul_dsf.cpp +++ b/src/pair_lj_cut_coul_dsf.cpp @@ -46,7 +46,6 @@ using namespace MathConst; PairLJCutCoulDSF::PairLJCutCoulDSF(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_lj_cut_coul_wolf.cpp b/src/pair_lj_cut_coul_wolf.cpp index 540cb5c959..f75422b684 100644 --- a/src/pair_lj_cut_coul_wolf.cpp +++ b/src/pair_lj_cut_coul_wolf.cpp @@ -37,7 +37,6 @@ PairLJCutCoulWolf::PairLJCutCoulWolf(LAMMPS *lmp) : Pair(lmp) { single_enable = 0; writedata = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_lj_expand.cpp b/src/pair_lj_expand.cpp index 6cb13f7390..bc6ef98785 100644 --- a/src/pair_lj_expand.cpp +++ b/src/pair_lj_expand.cpp @@ -32,7 +32,6 @@ using namespace MathConst; PairLJExpand::PairLJExpand(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_lj_gromacs_coul_gromacs.cpp b/src/pair_lj_gromacs_coul_gromacs.cpp index 30dc5b67fd..3f047d2334 100644 --- a/src/pair_lj_gromacs_coul_gromacs.cpp +++ b/src/pair_lj_gromacs_coul_gromacs.cpp @@ -26,7 +26,6 @@ #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ @@ -34,7 +33,6 @@ using namespace LAMMPS_NS; PairLJGromacsCoulGromacs::PairLJGromacsCoulGromacs(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_lj_smooth.cpp b/src/pair_lj_smooth.cpp index b6d4526f77..f45dbaa6bc 100644 --- a/src/pair_lj_smooth.cpp +++ b/src/pair_lj_smooth.cpp @@ -33,7 +33,6 @@ using namespace LAMMPS_NS; PairLJSmooth::PairLJSmooth(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_lj_smooth_linear.cpp b/src/pair_lj_smooth_linear.cpp index 2b494b14e7..54c0f2612d 100644 --- a/src/pair_lj_smooth_linear.cpp +++ b/src/pair_lj_smooth_linear.cpp @@ -25,14 +25,12 @@ #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ PairLJSmoothLinear::PairLJSmoothLinear(LAMMPS *lmp) : Pair(lmp) { single_hessian_enable = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_morse.cpp b/src/pair_morse.cpp index 3ec3db1e2b..6c67a8ef55 100644 --- a/src/pair_morse.cpp +++ b/src/pair_morse.cpp @@ -22,7 +22,6 @@ #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ @@ -30,7 +29,6 @@ using namespace LAMMPS_NS; PairMorse::PairMorse(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/pair_ufm.cpp b/src/pair_ufm.cpp index 45dcf3fbe2..8decaa42a7 100644 --- a/src/pair_ufm.cpp +++ b/src/pair_ufm.cpp @@ -28,7 +28,6 @@ #include "memory.h" #include "error.h" - using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ @@ -36,7 +35,6 @@ using namespace LAMMPS_NS; PairUFM::PairUFM(LAMMPS *lmp) : Pair(lmp) { writedata = 1; - centroidstressflag = 1; } /* ---------------------------------------------------------------------- */ diff --git a/src/respa.cpp b/src/respa.cpp index dcabcb11ef..aebee0f6eb 100644 --- a/src/respa.cpp +++ b/src/respa.cpp @@ -311,9 +311,10 @@ void Respa::init() if (force->pair && force->pair->respa_enable == 0) error->all(FLERR,"Pair style does not support rRESPA inner/middle/outer"); - // virial_style = 1 (explicit) since never computed implicitly like Verlet + // virial_style = VIRIAL_PAIR (explicit) + // since never computed implicitly with virial_fdotr_compute() like Verlet - virial_style = 1; + virial_style = VIRIAL_PAIR; // setup lists of computes for global and per-atom PE and pressure diff --git a/src/verlet.cpp b/src/verlet.cpp index 3d20a771f4..5e53daafbb 100644 --- a/src/verlet.cpp +++ b/src/verlet.cpp @@ -54,11 +54,12 @@ void Verlet::init() error->warning(FLERR,"No fixes defined, atoms won't move"); // virial_style: - // 1 if computed explicitly by pair->compute via sum over pair interactions - // 2 if computed implicitly by pair->virial_fdotr_compute via sum over ghosts + // VIRIAL_PAIR if computed explicitly in pair via sum over pair interactions + // VIRIAL_FDOTR if computed implicitly in pair by + // virial_fdotr_compute() via sum over ghosts - if (force->newton_pair) virial_style = 2; - else virial_style = 1; + if (force->newton_pair) virial_style = VIRIAL_FDOTR; + else virial_style = VIRIAL_PAIR; // setup lists of computes for global and per-atom PE and pressure